GWAS of peptic ulcer disease implicates Helicobacter pylori infection, other gastrointestinal disorders and depression

Abstract

Genetic factors are recognized to contribute to peptic ulcer disease (PUD) and other gastrointestinal diseases, such as gastro-oesophageal reflux disease (GORD), irritable bowel syndrome (IBS) and inflammatory bowel disease (IBD). Here, genome-wide association study (GWAS) analyses based on 456,327 UK Biobank (UKB) individuals identify 8 independent and significant loci for PUD at, or near, genes MUC1, MUC6, FUT2, PSCA, ABO, CDX2, GAST and CCKBR. There are previously established roles in susceptibility to Helicobacter pylori infection, response to counteract infection-related damage, gastric acid secretion or gastrointestinal motility for these genes. Only two associations have been previously reported for duodenal ulcer, here replicated trans-ancestrally. The results highlight the role of host genetic susceptibility to infection. Post-GWAS analyses for PUD, GORD, IBS and IBD add insights into relationships between these gastrointestinal diseases and their relationships with depression, a commonly comorbid disorder.

Introduction

Gastrointestinal (GI) diseases are highly prevalent in western countries. They use substantial health care resources, have a heavy societal economic burden^1,2, and impact the quality of life of those affected. GI disorders include peptic ulcer disease (PUD), gastro-esophageal reflux disease (GORD), irritable bowel syndrome (IBS), and inflammatory bowel disease (IBD), and some of these disorders are commonly reported as correlated with each other³. PUD, a common GI disorder, involves breaks (ulcers) in the inner lining of the digestive tract, usually located in the stomach or proximal duodenum. In GORD, the stomach contents leak back from the stomach into the esophagus⁴. IBS is a chronic functional disorder of the GI system. Patients with IBS often manifest abdominal pain and altered bowel habit, with either predominantly diarrhea, constipation, or both. IBD includes Crohn’s disease (CD) and ulcerative colitis (UC), which are chronic immune dysregulation disorders causing inflammation of the GI tract.

PUD is a complex disorder, for which Helicobacter (H.) pylori infection and the use of non-steroidal anti-inflammatory drugs (NSAIDs) are the main risk factors⁵. The development of infection-relevant PUD is recognized to be a multistep process, with contributions from both H. pylori infection and subsequent inflammation and damage of mucosa⁵. Eradicating H. pylori is effective for infection-relevant PUD treatment⁵. However, understanding the host factors influencing H. pylori infection and subsequent response could contribute to earlier risk identification and/or prevention, especially given the increasing antimicrobial resistance worldwide⁵. Moreover, clinical presentation of PUD that is not associated with H. pylori infection, nor with the use of NSAIDs, are now also imposing substantial diagnostic and therapeutic challenges^5,6. Lifetime prevalence of PUD in the general population has been estimated to be about 5–10%⁵. GORD is a multifactorial disorder and is more common in individuals with obesity and hiatal hernia⁷. Lifetime risk estimates of GORD have a wide range (9–26%), with a sample size-weighted mean of 15%⁸. An increase in the prevalence of GORD since 1995 has been reported⁸. IBS, a common disorder with a population lifetime risk of 11% globally⁹, is also likely a multifactorial disease, where hypervigilance of the central nervous system, immune activation of the intestinal mucosa, microbiome, prior infections, and diet are all suspected to play a role¹⁰. IBD is associated with many lifestyle risk factors, particularly smoking¹¹, and lifetime risk for IBD is around 0.3% in most countries of Europe¹². The genetic contributions to PUD, GORD, IBS and IBD have been well-recognized^13,14,15,16, and well-powered genome-wide association studies (GWASs) have identified >200 approximately independent susceptibility loci associated with IBD¹⁷. These loci implicate pathways such as autophagy and the IL-17/IL-23 axis and provide insights into IBD pathogenesis¹⁷. IBD has been extensively studied through the GWAS paradigm¹⁷ but to a lesser extent for PUD¹⁸, GORD^19,20, and IBS^21,22,23. Notably, the only GWAS to date for duodenal ulcer was in a Japanese ancestry cohort¹⁸.

Our primary focus was to identify genetic susceptibility factors for PUD by conducting a GWAS using data from the UK Biobank (UKB). Given the relationship between PUD and GORD (both acid-related disorders), the comorbidity of GORD and IBS³ and the much increased sample sizes afforded by the UKB, we also present GWAS results for GORD and IBS and investigate the shared genetic architecture between them. Inflammation of the gut is a key characteristic of IBD but not IBS²⁴, and hence differences between IBD and IBS from a genetic perspective are expected. Here, we take the opportunity to evaluate formally the genetic relationship between IBD and the more common disorders PUD, GORD, and IBS. Given increasing evidence for the importance of bidirectional signaling between the brain and the gut^25,26,27,28, possibly contributing to observational associations between depression and PUD²⁹, GORD³⁰, IBS³¹, and IBD³², we explore the potential causal relationships between major depression (MD) and the four disorders using Mendelian randomization (MR).

Here, we identify eight independent and significant single nucleotide polymorphisms (SNPs) for PUD and the results highlight the role of host genetic susceptibility to infection, acid secretion, and gastric motility. We also identify a total of 11 independent and significant SNPs previously unreported for GORD and IBS. We show genetic similarity across PUD, GORD, and IBS and between these GI diseases with psychiatric disorders. Post-GWAS analyses find a link between PUD, GORD, and IBS with the nervous system. We find a significant relationship between PUD, GORD, IBS, and major depression through observational and MR analyses. Taken together, our results expand our understanding of the role of genetics in gastrointestinal diseases and add insights into relationships between these gastrointestinal diseases and their relationships with MD.

Results

The workflow for our study is given in Supplementary Fig. 1.

Prevalence and comorbidity

Based on disease-diagnosis (self-reported, hospital admission, primary care, and death register records) in the UKB, four case–control digestion disorder datasets were identified (Table 1), with prevalences for PUD, GORD, IBS and IBD of 3.7%, 12.0%, 6.5%, and 1.5%, respectively (Table 1). We note that the UKB prevalence rates are higher for IBD than the 0.3% aforementioned lifetime risk¹². This may reflect ascertainment biases or misdiagnosis within the UKB. Disease statistics (Table 1) and comorbidity analyses (Fig. 1) help describe the UKB phenotypes. The male/female odds ratio for being a PUD case is 1.66 while for IBS it is 0.42 (Table 1), consistent with PUD being more common in men and IBS more common in women. To allow readers to interpret the phenotype data that underpin our analyses, we describe in detail the disorder co-occurrences present in UKB participant records (Fig. 1 and Supplementary Note 1), which may reflect, in part, the natural course of the symptom presentations and/or misdiagnosis. For example, it is well recognized that IBD patients may have symptoms that provide a PUD or GORD diagnosis prior to, or as part of, the endoscopies that generate an IBD diagnosis and/or monitor the disease progress. Moreover, GORD can occur as a consequence of treatment of H. pylori for PUD as H. pylori infection can reduce gastric acidity³³. Here, while rates of PUD and IBD in those with GORD were significantly lower than the rate of GORD cases in the UKB as a whole, those with PUD or IBD were significantly more likely to also have GORD. For each of the PUD, GORD, IBS, and IBD (defined as the index disease), competitive comorbidity analyses tested, among the other three diseases, which disease is more prone to be comorbid with the index disease. We found that PUD is more prone to be comorbid with GORD while IBS is more likely to be comorbid with IBD (Fig. 1c).

Table 1 Full-sibling relative risk and heritability estimation for PUD, GORD, IBS, and IBD.

Fig. 1: Comorbidity analyses in unrelated European individuals and heritability estimation based on full-sibling relative risk for PUD, GORD, IBS, and IBD.

a The number of unrelated individuals with each diagnosis (cyan boxes) and the number of overlapped individuals between each pair of PUD, GORD, IBS, and IBD cases (yellow boxes). These differ from Table 1 because only unrelated individuals are considered here. b Cells represent ratio of the odds of disease cases from each column in those with disease from each row and the odds of each row disease cases in unrelated European-ancestry individuals. The diagonal elements are the sample risk rates and differ from “N case/(N case + N control)” in Table 1 because only unrelated individuals are included here. c Competitive comorbidity analyses for each of PUD, GORD, IBS, and IBD disease (defined, in turn, as index disease) to test among the other three diseases which disease is more prone to be comorbid with the index disease. The disease on the left of the red dashed line is the index disease and the number shows the number of cases without any comorbidity. The corresponding Venn diagram shows the number of individuals with recorded or self-report diagnosis with at least one of the other three diseases. The numbers outside of the Venn diagrams are the numbers of individuals with diagnosis of both the index disease and each of the other three diseases in turn, after removing the overlapped individuals for these three diseases. At the bottom of each Venn diagram is the proportion of the index disease cases in the other three diseases respectively, compared in pairs using a two-proportion Z test to test which disease is more prone to be comorbid with the index disease. d Heritability estimates and 95% confidence interval for PUD, GORD, IBS, and IBD based on full-sibling relative risk. See the Supplementary Note 1 for an example explaining b, c.

Full-sibling risk and heritability estimation

Using SNP-based estimates of coefficients of genetic relationship between individuals in the UKB, we estimated the full-sibling relative risk for each of the PUD, GORD, IBS, and IBD and the heritability of liability, with the assumption that the increased risk in relatives only reflect shared genetic factors (Table 1 and Fig. 1d). The estimated heritabilities for PUD, GORD, IBS, and IBD were 0.28 (95% CI: 0.18–0.37), 0.28 (95% CI: 0.23–0.33), 0.12 (95% CI: 0.05–0.19), and 0.49 (95% CI: 0.36–0.64), respectively, all significantly different from zero.

GWAS

Genome-wide association analyses were conducted for five digestion phenotypes, the four disease-diagnosis traits (PUD, GORD, IBS, and IBD, see “Methods” section) and one trait that combined the disease-diagnosis and taking of corresponding medications and treatments (Supplementary Tables 1, 2). In clinical practice, medications for PUD also have a therapeutic effect on GORD, hence we generated the PG₊M phenotype—a combination of disease-diagnosis of GORD and/or PUD and/or corresponding medications and treatments. We tested for association between 8,546,065 DNA variants and each of the five digestion phenotypes (PUD, GORD, PG₊M, IBS, and IBD) in 456,327 UKB participants. A total of 66 within-trait (61 across-trait) independent variants were genome-wide significant (P < 5.0E−8) for the five digestion phenotypes analysed, of which 8 were associated with PUD, 6 with GORD, 19 with PG₊M, 2 with IBS, and 31 with IBD. Two PUD-associated SNPs have been previously linked to PUD but in a Japanese cohort¹⁸. For GORD, a recent GWAS study²⁰ which used an earlier release of UKB data plus meta-analysis with other studies reported 25 genome-wide significant SNPs. Conditional analyses of our GORD and PG₊M association results on SNPs previously reported as genome-wide significantly associated gastroesophageal reflux²⁰ using GCTA-COJO^34,35 found that three SNPs remain significant for GORD and six SNPs remain significant for PG₊M. Given the aim of our study, Tables 2–3 list the 19 previously unreported genome-wide significant SNPs for PUD, GORD, PG₊M, and IBS. SNPs associated with IBD are in Supplementary Table 3; we note that 28 of the 31 SNPs have been previously linked with inflammatory bowel diseases (despite concerns of the high prevalence of IBD in the UKB). The three GORD-associated SNPs and 13 PG₊M-associated SNPs that have been reported²⁰ are listed in Supplementary Table 4. The GERA cohort data^36,37 were available as a replication sample for PUD (1004 cases, 60,843 controls) and IBS (3359 cases and 58,488 controls), albeit with limited power (~20% for each SNP, Supplementary Table 5). Six of the eight genome-wide significant SNPs for UKB PUD exist and one for UKB IBS in GERA and all have very similar effect size estimates as those in UKB, but only rs681343 was formally significant (P = 5.0E−4 < 0.05/8, a replication rate consistent with the power calculations) (Supplementary Table 5).

Table 2 Genome-wide significant SNPs associated with PUD in the UK Biobank.

Table 3 Genome-wide significant SNPs associated with GORD, PG₊M, and IBS in the UK Biobank.

Many of the PUD SNPs associations are previously unreported even when checking published results for any digestive diseases (Table 2). There are eight PUD-associated SNPs, which by physical distance (<22 kb) could implicate MUC1, MUC6, FUT2, PSCA, ABO, CDX2, GAST, and CCKBR, as annotated in Fig. 2. Pathway analysis using the GENE2FUNC of the FUMA pipeline³⁸ showed that these genes are highly overexpressed in human stomach tissue. Detail on these PUD results are found in the discussion section. For IBS, a previous study of UKB data (earlier phenotype release) reported rs10512344 as the only SNP genome-wide significantly associated (P = 3.6E−8) with self-report IBS, and in females only²³. The P value for this SNP in our analyses is 5.0E−05 (P = 4.4E−07 in females, i.e., less associated than in the published study despite sample overlap). Given that they used a self-report phenotype while we used a combination of three resources, we regenerated three subgroup phenotypes (self-report only, hospital admission only, and primary care only) for IBS and repeated GWAS analyses for the three phenotypes and the results were used for further analyses (Supplementary Note 2). We also conducted similar subgroup phenotypes GWAS for GORD and PUD. Figure 2 shows Manhattan plots for PUD. Figure 3 shows Manhattan plots for GORD, PG₊M, and IBS and Supplementary Fig. 2 shows Manhattan plots for IBD. Quantile–Quantile (Q–Q) plots of all the variants analysed in UKB are provided in Supplementary Fig. 3 for the five phenotypes. Regional visualization plots of the 66 independent variants are in Supplementary Figs. 4–8. Detailed pleiotropy results derived from the GWAS Catalog³⁹ are provided in Supplementary Data 1. The 54 GTEx tissue specific enrichment results for eight genes in Fig. 2 are in Supplementary Data 2. In sensitivity analysis GWAS, phenotypes for each of PUD, GORD, IBS, and IBD were regenerated by excluding individuals with more than one GI diagnosis among these four disorders from the corresponding original phenotype cases (Supplementary Note 3).

Fig. 2: Manhattan plot of peptic ulcer disease (PUD) for SNPs associated P < 1.0E−4 from BOLT-LMM association test.

SNPs highlighted with green triangles are independent loci with P < 5.0E−8. rs2976388 and rs687621 are the only two loci associated with duodenal ulcer in a Japanese cohort¹⁸ and rs681343 showed statistically significant association in GERA PUD GWAS, as annotated in the blue box. SNPs on odd/even chromosomes are presented in mauve/pink. Schematic diagram on the right side represents the reported biological evidence supporting involvement in peptic ulcers of genes physically located near PUD-associated loci (noting that we do not have direct evidence to link the associated SNPs with the genes). MUC1⁷⁶, MUC6⁷⁷, and FUT2⁷³ have been linked to susceptibility to H. pylori infection and PSCA and ABO have been proposed to be associated with subsequent response after infection¹⁸. Induced/enhanced CDX2 expression as a result of H. pylori infection of gastric epithelial cell lines has been observed⁷⁹. GAST encodes gastrin, which is a hormone whose main function is to stimulate secretion of hydrochloric acid by the gastric mucosa. CCKBR encodes cholecystokinin receptor which mediates a therapeutic effect for peptic ulcer treatment by reducing acid secretion and inhibiting gastrointestinal motility. The cholecystokinin receptor is also an effect-mediating target of itriglumide on phase II clinical trial for anxiety and panic disorder⁸¹. Using the GENE2FUNC of FUMA pipeline³⁸, we found that the eight annotated genes are highly overexpressed in human stomach tissue from GTEx 8th version data, as shown in the bar plot on the right side with data presented in Supplementary Data 2. The dashed line represents Bonferroni corrected significance at −log₁₀(0.05/54).

Fig. 3: Manhattan plot of gastroesophageal reflux disease (GORD), peptic ulcer disease, GORD and corresponding medications (PG₊M) and irritable bowel syndrome (IBS) for SNPs associated P < 1.0E−4 from BOLT-LMM association test.

a GORD: SNPs highlighted with orange squares are genome-wide statistically significant (P < 5.0E−8) independent loci, which correspond to the orange squares in b. SNPs highlighted with yellow represent loci that have not previously been reported to be associated with GORD. b PG₊M: SNPs highlighted with red diamond are independent loci with P < 5.0E−8. Only two of eight SNPs associated with peptic ulcer disease (highlighted with green triangles) are with P < 1.0E−5 in PG₊M and all the SNPs associated with GORD (highlighted with orange squares) have P < 1.0E−5 in PG₊M. SNPs highlighted with yellow represent loci that have not previously been reported to be associated with GORD. c IBS: SNPs highlighted with red diamond are independent loci with P < 5.0E−8. The blue dot is for rs10512344 that has been reported associated with female IBS in UKB previously released data²³. The IBD Manhattan plot is found in Supplementary Fig. 2.

SNP-based heritabilities and genetic correlations

We estimated the proportion of variance in trait liability attributable to genome-wide common SNPs (i.e., SNP-based heritability, $h_{{\rm{SNP}}}^2$) using Linkage Disequilibrium Score regression (LDSC)⁴⁰. $h_{{\rm{SNP}}}^2$ estimates on the liability scale were: PUD 0.06 (SE = 0.007), GORD 0.08 (SE = 0.004), PG₊M 0.09 (SE = 0.003), IBS 0.06 (SE = 0.005), and IBD 0.11 (SE = 0.016) (Fig. 4a and Supplementary Table 6), all significantly different from zero. The SNP-based genetic correlation (r_g) between PUD and GORD is 0.65 (SE = 0.05, P = 4.9E−36, phenotypic correlation (r_p) = 0.11), similar to the r_g estimate for GORD and IBS (0.65, SE = 0.05, P = 1.1E−46, r_p = 0.10). The r_g between PUD and IBS is 0.49 (SE = 0.08, P = 2.0E-10, r_p = 0.03) (Supplementary Table 7), while the r_g between IBD and each of PUD, GORD, PG₊M, and IBS are not statistically significantly different from zero after Bonferroni correction (Fig. 4b). In sensitivity analyses, all individuals with more than one GI diagnosis were excluded (Supplementary Fig. 9a); the $h_{{\rm{SNP}}}^2$ estimates were lower but still significantly different from zero (Supplementary Table 6 and Supplementary Fig. 9b), while the r_g between PUD, GORD, and IBS remain significant (all r_g > 0.25) and none showed statistically significant r_g with IBD (Supplementary Table 7 and Supplementary Fig. 9c₁). Detailed results of sensitivity analyses are discussed in Supplementary Note 3 with the corresponding data presented in Supplementary Tables 6–7 and Supplementary Fig. 9. We also investigated $h_{{\rm{SNP}}}^2$ and r_g for the three subgroup phenotypes (mentioned above) of each of PUD, GORD, and IBS, and found them to show significant $h_{{\rm{SNP}}}^2$ and high genetic correlation (see Supplementary Fig. 10 and Supplementary Tables 8–9).

Fig. 4: SNP-based heritability and genetic correlation analyses for the five digestion phenotypes from LD score regression analyses.

a SNP-based heritability of the five digestion phenotypes both on the observed and liability scales. The transformation to the liability scale uses the UKB sample risk, i.e., the proportion cases in the UKB cohort, as the population lifetime risk; the sample risk percentage of whole UKB European ancestry is shown below the x-axis in parentheses. The error bars represent 95% confidence interval for the estimated SNP-based heritability. “*” represents that the SNP-based heritability P value remain significant after Bonferroni correction (P < 0.05/10). b, c Results for genetic correlation within-digestion phenotypes and between digestion phenotypes and the nine psychiatric and neurological disorder traits are in b, while genetic correlation results between digestion phenotypes and traits from LD Hub that had a significant correlation with at least one of the digestion phenotypes are provided in b, c. “*” represent that genetic correlation estimates are still significant after Bonferroni correction (P < 0.05/((5 + 9 + 258) * 5)).

The r_g between each of the five digestion phenotypes and the 258 traits from LD Hub⁴¹ plus nine psychiatric^{42,43,44,45,46,47,48} and neurological traits^49,50 (Supplementary Table 10 and Supplementary Data 3, respectively) included 26, 37, 53, 16, and 3 significant correlations for PUD, GORD, PG₊M, IBS, and IBD, respectively, after Bonferroni correction (P < 3.7E−5). Figure 4b, c show the r_g between each of the five phenotypes and statistically significant correlations from the 267 traits. We summarize these results briefly. First, we observed significant positive r_g between four digestion phenotypes (PUD, GORD, PG₊M, and IBS) and depressive symptoms⁵¹, major depression (MD)⁴⁸, attention deficit hyperactivity disorder (ADHD)⁴², neuroticism⁵¹, and insomnia^52,53. However, there was no statistically significant r_g between IBD and these traits. Second, PUD, GORD, and PG₊M have significant positive r_g with body mass index (BMI)⁵⁴, body fat-related traits⁵⁵ and coronary artery disease⁵⁶. GORD and PG₊M are also genetically correlated with type 2 diabetes⁵⁷ in a positive direction. Third, the common variant genetic architecture of PUD, GORD, PG₊M, and IBS are positively correlated with educational attainment-related traits⁵⁸. In sensitivity analyses, all individuals with more than one GI diagnosis were excluded and the results are similar as above. Detailed results of sensitivity analysis are discussed in Supplementary Note 3 with the corresponding data presented in Supplementary Table 11, Supplementary Fig. 9, and Supplementary Data 4.

Links to gene expression, eQTLs, and mQTLs

We used partitioned SNP-based heritability analyses to determine if any genomic annotations showed enrichment of $h_{{\rm{SNP}}}^2$ relative to the null hypothesis that SNP-based heritability is partitioned proportional to the number of SNPs in an annotated set (see “Methods” section). After Bonferroni correction, GORD, PG₊M, and IBS showed significant enrichment of $h_{{\rm{SNP}}}^2$ in conserved regions and $h_{{\rm{SNP}}}^2$ enrichment for IBD was in the super enhancer category (Supplementary Fig. 11 and Supplementary Table 12). In analyses based on SNP annotations derived from cell-type histone-mark data (Fig. 5a and Supplementary Table 13), IBD showed significant $h_{{\rm{SNP}}}^2$ enrichment in immune and gastrointestinal cell-type groups, while GORD, PG₊M, and IBS showed enrichment in the CNS cell-type. Based on cell-type specific SNP annotations⁵⁹ derived from gene expression data of 205 different tissues (53 from GTEx⁶⁰ and 152 from Franke lab⁶¹), PG₊M showed significantly enriched association with genes expressed in the hippocampus, frontal cortex (Brodmann Area, BA9) and anterior cingulate cortex (BA24) of the brain and IBD showed enriched associations in leukocytes (Fig. 5b and Supplementary Table 14). Given these significant results for GORD, PG₊M, and IBS we conducted the same analyses using the fine-scale GTEx brain gene expression data which includes data from 13 brain regions (albeit smaller sample size)⁶⁰. GWAS associations for PG₊M were consistently enriched in the frontal cortex (BA9) (Fig. 5c and Supplementary Table 15). We conducted sensitivity analyses to investigate if these results are mediated by educational attainment (EA)⁶², BMI⁵⁵, and smoking-related traits⁶³. The results remain significant after conditioning the PG₊M GWAS results on the GWAS results of EA, BMI and smoking using mtCOJO³⁷ (Supplementary Tables 13–15). We also investigated whether associations between SNPs and the five digestion phenotypes were consistent with mediation through gene expression using the Summary data-based Mendelian randomization method, SMR⁶⁴. SMR combines genome-wide significant SNP-gene expression associations (i.e., eQTLs) with the SNP-trait association results. Significant SMR associations provide the best statistical evidence given available data that the trait associated SNP could be a causal SNP through its action on gene expression. A total of five unique genes for which expression is significantly associated with three digestion phenotypes, including two genes for PUD (PSCA and FUT2), one gene for PG₊M (SUOX), and one gene for IBD (RP11-129J12.2 and RPS23P10), were identified (Supplementary Table 16). The statistical framework of SMR⁶⁴ can be applied to mQTL (genome-wide significant SNP-methylation association) data to identify putative methylation-trait association. Hence, we repeated SMR analyses for PUD using blood mQTL data from McRae et al.⁶⁵. Among the mQTL SMR results (Supplementary Table 17), three DNA methylation probes (cg01656853, cg08873673, and cg04660111), located in promoter region of FUT2, are associated with PUD (Supplementary Fig. 12). From Roadmap epigenomics annotation⁶⁶, there is a digestive tissue specific active enhancer for FUT2, as shown in Supplementary Fig. 12.

Fig. 5: Analyses of partitioning SNP-based heritability by different annotations.

a SNP-based heritability enrichment analysis for each digestion phenotype partitioned by cell type groups annotated by histone marks. The dashed line represents the Bonferroni corrected significance threshold (0.05/(53 * 5 + 5 * 5)). b SNP-based heritability enrichment analysis for PG₊M and IBD partitioned by cell types annotated using cell-type specific gene expression data given the results from a. The dotted lines represent the Bonferroni correction threshold (P < 0.05/(205*2)). c SNP-based heritability enrichment analysis for GORD, PG₊M, and IBS partitioned by cell types annotated using fine-scale GTEx brain gene expression data from 13 brain regions selected given the central nervous system enrichment results of these three phenotypes from a. The dashed line represents the Bonferroni correction threshold (P < 0.05/(13*3)).

Gene-based and gene-set enrichment analyses

While individual SNP associations can implicate relevance of specific genes, given gene-specific genetic architectures putative roles of genes can also be identified by gene-based tests that combine the SNP-associations into genic annotations. We used MAGMA⁶⁷ software to identify genes significantly (P < 2.7E−6) associated with each of the five digestion phenotypes, identifying: 29 for PUD, 112 for GORD, 157 for PG₊M, 14 for IBS, and 97 for IBD (Supplementary Data 5). For gene-set enrichment analysis applied to gene-based summary statistics identified 11 gene ontology (GO) sets as significant for PG₊M and six for IBD. The top enriched gene set for PG₊M is “GO: POSTSYNAPSE” while that for IBD is “GO MHC CLASS II PROTEIN COMPLEX” (Supplementary Table 18).

Comorbidity and MR with depression

Given observational associations between depression and PUD²⁹, GORD³⁰, IBS³¹, and IBD³², we used the eight UKB depression phenotypes identified by Cai et al.⁶⁸, together with our four digestion diagnoses, to test whether each of the 32 depression-digestion phenotypes pairs show significant lifetime comorbidity relationships. Cai et al.⁶⁸ argued that some of the clinically derived depression phenotypes were more specific to major depressive disorder than self-report depression, and that biological interpretation should focus on analyses using clinically derived phenotypes. We were interested to know if the relationship between depression and common GI disorders depended on these different depression definitions. All eight depression phenotypes showed statistically significant comorbidity relationship with each of PUD, GORD, and IBS (Supplementary Tables 19, 20 and Fig. 6a). For IBD, the highest statistically significant ORs were for the more severe depression definitions of electronic medical record depression (ICD10Dep), DSM-V clinical guideline-defined major depressive disorder and recurrence (abbreviation: LifetimeMDD and MDDRecur) (Fig. 6a). Given the comorbidities of the digestion phenotypes with depression, we tested for statistical evidence of a putative causal association between MD and each of the five digestion phenotypes using Generalized Summary-data-based MR (GSMR)³⁷; we also tested for reverse causality (Supplementary Table 21 and Supplementary Figs. 13, 14). For the relationship between MD and IBD, GSMR estimates were not statistically significant either in forward direction (the effect of MD on IBD) or the reverse direction (the effect of IBD on MD), despite being well-powered to identify associations. We found a unidirectional effect between MD and PUD, 1.20-fold risk for PUD per standard deviation (SD) in liability to MD. The reverse direction (PUD as exposure and MD as outcome) had a small point estimate and is not statistically significant, but we note that we relaxed the association significance threshold to obtain more PUD genetic instruments (8 SNPs to 13 SNPs). When we repeated the analyses (PUD as exposure and MD as outcome) using the eight genome-wide significant SNPs, the result remained statistically non-significant, suggesting that MD is putatively causal for PUD. This analysis should be revisited when GWAS sample size for PUD increases identifying more genome-wide significant SNP instruments. The effect of MD on GORD and IBS showed statistically significant estimates of 1.23-fold and 1.48-fold respectively increase per SD in liability to MD (Fig. 6b). The point estimates for the reverse causality analyses were smaller (but statistically significant), and again these analyses should be revisited in when more genome-wide significant SNPs are identified (Fig. 6b and Supplementary Table 21). We observed bidirectional statistically significant results between MD and PG₊M, i.e., 1.27-fold increased risk for PG₊M per SD) in liability to MD (P = 2.5E−15), and 1.26-fold increased risk for MD per SD in liability to PG₊M (P = 2.7E−09). No SNPs were identified as outliers by the HEIDI test. The pattern of results was the same when other MR methods were applied, which, as expected, showed less significant results (see Supplementary Note 4, Supplementary Table 22, and Supplementary Fig. 15). The bidirectional MR statistical significance between MD and PG₊M could be consistent with reverse causality or pleiotropy. To exclude potential known confounders we repeated analyses with GWAS summary statistics conditioned on GWAS summary statistics of on EA, BMI, and smoking-related traits from mtCOJO analysis. The results remain significant in both directions (Supplementary Table 23). We used the latent causal variable (LCV) method⁶⁹, which is designed to better separate pleiotropy from causality. As expected, the genetic causality proportion is not significant for PG₊M and MD because of the strong bidirectional significance (Supplementary Note 4 and Supplementary Table 24). As sensitivity analyses, we repeated the analyses after removing the depression cases from both cases and controls from the five GI disorder phenotypes (Supplementary Note 5 and Supplementary Tables 25–27). In summary, the statistically non-significant GSMR results in this study neither support a causal relationship from MD to IBD nor from IBD to MD. The results for PG₊M do not support a unidirectional causal relationship given the significant bidirectional GSMR results. The MR results between MD and PUD, support a putative causal role of MD in PUD, but should be revisited when more SNP instruments are available.

Fig. 6: Comorbidity, Mendelian randomisation and polygenic risk score analyses with depression.

a Comorbidity relationship between eight depression phenotypes⁶⁸ and each of the PUD, GORD, IBS, and IBD. Odds ratio and 95% CI were calculated based on the 2 by 2 contingency table for number of cases and controls of each of the 32 digestion-depression phenotype pairs (Supplementary Tables 19 and 20). The definitions of the eight depression phenotypes are in the Methods section. “*” represents the P value for each digestion and depression pair that remain significant after Bonferroni correction (P < 0.05/32). b Mendelian randomization (MR) results between major depression (MD) and five digestion phenotypes. The left y axis is for MD while the right y axis is for the five digestion phenotypes. The arrow for each horizontal line represents the direction from exposure trait to outcome trait relative to the y-axis labels. OR and 95% CI are represented as diamond and horizontal lines taking values from x axis. Each digestive phenotype corresponds to two horizontal lines. “R” on the right side of the horizontal line represents relaxation of significance threshold of SNP associations to obtain more SNP genetic instruments (Supplementary Table 21). The number of the SNP instruments used in analyses are shown above the diamond. The common pathological characteristics or symptoms for these phenotype-related diseases are shown on the right side (noting that pathological characteristics or symptoms are not limited to these locations). c Decile of major depression polygenic score predicts PG₊M. Odds ratio (OR) and 95% confidence intervals (CI, blue dots and orange bars) relative to decile 1 were estimated using logistic regression. The blue dashed lines shows the OR for the highest decile, OR of 1.33 for PG₊M related disorders. d OR and 95% CIs (green dots and blue bars) for top decile of PG₊M polygenic score predicting depression and psychiatric disorders in the independent GERA cohort³⁶.

Out-of-sample polygenic score prediction

To further investigate the relationship between MD and PG₊M, we used MD GWAS summary statistics (European ancestry, excluding the UKB cohort)⁴⁸ to generate MD polygenic scores and used these to predict PG₊M risk in the UKB. We found that participants in the UKB with a high polygenic score for MD had a higher risk for PG₊M-related disorders. The top decile of individuals ranked on polygenic risk prediction for MD had an OR of 1.33 (95% CI: 1.29–1.38) for PG₊M risk compared to the bottom decile (Fig. 6c and Supplementary Table 28). We also selected genome-wide significant SNPs associated with PG₊M in the UKB to calculate polygenic scores and predict depression and psychiatric disorder risk in GERA cohort³⁶, as shown in Fig. 6d respectively. Further, we selected genome-wide significant SNPs associated with PUD, IBS in the UKB to calculate polygenic scores and predict peptic ulcer and IBS risk in GERA cohort³⁶. The top decile of individuals ranked on polygenic score for PUD had an OR of 1.80 (95% CI: 1.38–2.36) for PUD risk compared to the bottom decile and that for IBS had an OR 1.42 (95% CI: 1.21–1.67) for IBS risk compared to bottom decile (Supplementary Fig. 16).

Discussion

This study describes an analysis of PUD and its relationship with other digestion disorders using a single study cohort. We used both the phenotypes and genotypes of up to 456,327 individuals to study the genetic contributions to PUD, GORD, IBS, and IBD and the connection between these disorders with major depression. Our results provide insights into the genetic basis of, and inter-relation between, these gastrointestinal disorders and also their relation to depression.

GWAS of PUD identified eight independent associated loci and 6 of 8 have potential links to H. pylori infection, highlighting the role of host genetic susceptibility (Fig. 2). Only two SNPs (rs2294008 and rs505922) have been previously reported as associated with peptic ulcers, both from a Japanese cohort study¹⁸ of duodenal ulcers. The authors provided evidence for their role in with duodenal ulcer development after H. pylori infection¹⁸. In our PUD-associated SNPs, rs2976388, located in PSCA gene, is in high LD (r² = 0.94 in Europeans) with rs2294008, while rs687621, an intronic SNP located in the blood-group ABO gene, is in high LD (r² = 0.98 in Europeans) with rs505922 (Fig. 2). These two loci are reported here in Europeans with the same direction of effect. From published gene expression data, we found that allele A of SNP rs2976388 is associated with increased PSCA expression (b_eQTL = 0.73, P_eQTL = 8.8E−41), and through SMR analysis, the increased expression of PSCA decreased risk for PUD (b_SMR = −0.12, P_SMR = 3.0E−11). Decreased PSCA expression has been reported following H. pylori infection⁷⁰, indicating negative regulation of PSCA expression by H. pylori infection. Other data sets recorded for H. pylori infection status are needed to explore this proposed relationship. A previous study⁷¹ has shown that individuals with blood group O have a higher risk of peptic ulcers compared to those with other blood groups which offers a possible explanation that this observation may result from different susceptibilities and immunologic responses to H. pylori infection. The A allele of rs687621 is associated with PUD (OR = 1.08, Table 2), and this allele is in high LD (r² = 0.86 in Europeans and r² = 0.89 in Japanese) with the rs8176719 deletion allele which generates a premature stop codon that leads to the O blood group⁷². We note that rs505922 has similar LD (r² = 0.87) with rs8176719 in Europeans, but is in slightly higher LD with rs8176719 in Japanese (r² = 0.92). This finding deserves further investigation of the association between blood group and peptic ulcer. Of the six previously unreported PUD-associated SNPs, it is notable that rs681343 is statistically significant in both UKB discovery and GERA replication GWAS (a replication rate consistent with power, see Supplementary Table 5). This SNP is located in the FUT2 gene which has been implicated in susceptibility to H. pylori infection in humans⁷³ and mice gastric tissue⁷⁴. The FUT2 fucosyltransferase allows expression of blood group antigens on the gastrointestinal mucosa and in bodily secretions. European individuals who are homozygous (AA) at rs601338 (non-secretor phenotype) are unable to secrete blood group antigens into bodily fluids, or express them on mucosal surfaces given the allele AA encodes a stop codon that inactivates the FUT2 enzyme⁷⁵. Moreover, rs601338 is also associated with different infections (either resistance or susceptibility)⁷⁵ and rs681343 associated with PUD in our analyses is in high LD with rs601338 in Europeans (r² = 0.996), with the T allele of rs681343 paired with the A allele of rs601338.

Another two of the PUD-associated SNPs are in located in genes encoding mucin peptides, located on different chromosomes: rs147048677 (P = 9.0E−12) is a synonymous variant in the MUC1 gene. A mouse model study⁷⁶ has shown that Muc1 limits H. pylori colonization of gastric mucosa. rs78459074 (P = 2.6E−10) is an intronic variant in the MUC6 gene. We note that, MUC6, MUC2, and MUC5AC are sequentially co-located on chromosome 11 (Supplementary Fig. 4). The expression of MUC6 and MUC5AC have been linked to H. pylori infection^77,78. For the other three PUD-associated SNPs, rs9581957 is 22 kb away from the CDX2 gene (Supplementary Fig. 4). In nonneoplastic cell lines H. pylori infection led to prompt and robust expression of Cdx2⁷⁹. rs10500661 (P = 4.1E−14) is located in ~7 kb upstream of CCKBR (cholecystokinin B receptor) and this gene encodes a G-protein coupled receptor for both gastrin and cholecystokinin, regulatory peptides of the brain and gastrointestinal tract⁸⁰. This gene, as shown in the GTEx⁶⁰ portal (https://gtexportal.org/home/gene/CCKBR), is highly expressed in the brain frontal cortex (Brodmann area (BA) 9) and stomach. Moreover, this gene is a therapeutic-effect target gene for a peptic ulcer drug proglumide (ATC code: A02BX06), which inhibits gastrointestinal motility and reduces gastric acid secretions. In addition to these findings, itriglumide, an antagonist for the CCKBR-encoded protein, has been investigated as a potential treatment for anxiety and panic disorders⁸¹. rs34074411 is located approximately 2 kb upstream of the GAST gene which encodes gastrin, a hormone whose main function is to stimulate secretion of hydrochloric acid by the gastric mucosa. It is unusual to be able to link 6 out of 8 GWAS associated loci for a trait to a putative mechanism, but here the results highlight the role of host genetic variability to bacterial infection. Given the high genetic correlation between PUD and GORD (0.65, SE = 0.05), we also conducted mtCOJO analyses for PUD conditioning on GORD GWAS summary statistics. All the SNPs mentioned above still showed statistical significance (Supplementary Table 29), suggesting that most of those PUD-associated SNPs are specific to PUD, and may provide important disorder specific profiles for use in clinical trials.

We identified 6, 19, and 2 independent significant associations for GORD, PG₊M, and IBS, respectively. Some genes around these loci have biological support for their mechanistic involvement and are worthy of note. For PG₊M, two of the significantly associated SNPs, 19:18793695 and rs13097265, have been previously linked to Barrett’s esophagus and esophageal adenocarcinoma⁸² (Supplementary Table 4). rs13097265 is located ~60 kb from the FOXP1 gene, and recently a mouse study⁸³ has shown that Foxp1 protein is expressed in all layers of the murine GI tract (including the myenteric plexus, which is part of the enteric nervous system and regulates gut peristalsis and transit). Altered motility and achalasia has been observed in Foxp1^+/− mice. Heartburn and regurgitation, the main GORD symptoms, commonly occur during the early stages of achalasia and are consequently poor indicators of esophageal motility disorder. It remains controversial whether GORD and achalasia coexist or whether one disease transforms into the other⁸³. Interestingly, rs10891491 associated with GORD (P = 4.1E−08) and rs7947502 associated with IBS (P = 2.5E−08) are in LD (r²: 0.15) (Table 3). rs10891491 and rs7947502 are intronic variants of NCAM1, suggesting this genic region may show pleiotropic effect for these two disorders. NCAM1 encoded protein is involved in development of the nervous system, including the enteric nervous system⁸⁴. These GWAS results deserve further investigation to understand the role of genetic variants in the etiology of PUD, GORD, and IBS, especially since many of the genes putatively implicated by SNP associations have biological mechanism support for these disorders.

We provide direct genetic evidence for the difference between IBD and the other digestion phenotypes; although not unexpected, our analyses quantify these differences as illustrated by high genetic correlations among PUD, GORD, and IBS, which all show low genetic correlations with IBD (Fig. 4b). Both PUD and GORD are acid-related diseases; their high genetic correlation (r_g = 0.65, s.e. = 0.052) motivated the combination of GORD and PUD with medication-taking cases. Additional evidence for genetic differences between IBD and the other digestion phenotypes was provided by the partitioned SNP-based heritability analyses, which showed enrichment of PG₊M-associated SNPs in genes expressed in the brain regions (hippocampus, BA9 and BA24) while those for IBD are enriched in blood and immune related tissues (Supplementary Table 14). SNP-based heritability for PG₊M is still enriched in BA9 region when we used fine-scaled GTEx brain gene expression data (Supplementary Table 15), and also after controlling for the effect of EA, BMI, and smoking (Supplementary Tables 13–15). Gene set enrichment analysis showed PG₊M-associated SNPs are enriched in neuron-related gene sets. We note that a limitation of our brain enrichment analysis is our conclusions are limited by the availability of tissue specific gene expression data. The GTEx database does not report gene expression data for multiple cortical regions, so specificity to BA9 or frontal cortex (over and above other cortical regions) is not established. Given the non-availability of gene expression data from other human tissues, such as sympathetic, parasympathetic (vagus nerve)²⁶, and enteric nervous system^25,27, we cannot conduct key hypothesis-based enrichment analyses. However, despite these limitations, our findings indicate that a genetic contribution to PG₊M may highlight the potential link between the nervous system and esophagus, stomach and duodenum^28,85, although there is likely not just one causal tissue or cell type⁵⁹. Historically, vagotomy was used commonly to manage peptic ulcer diseases, as vagal stimulation promotes acid secretion⁸⁶ (now successfully treated by H2 receptor agonists), indicating the clinical importance of the link between the nervous system and gastrointestinal tract.

We observed significant increased risk of the four digestion disorders with depression (under multiple definitions of clinical and self-report) (Fig. 6a). The association between mental health and GORD has been addressed through observational studies^30,87. For example⁸⁷, a bidirectional association between GORD and depression was reported with risk factor roles for depression on GORD and for GORD on depression. Here, we conducted MR to investigate if this statistical approach provides evidence of a causal relationship between major depression and GI disorder phenotypes. Since, the number of genome-wide significant SNP instruments for GORD and PUD is too few to draw strong conclusions, we focus on the combined PG₊M phenotype that identified a higher number genome-wide significant SNPs for use as the MR instruments. We found an OR of 1.27 (P = 2.5E−15) for PG₊M per SD in liability to MD, which has direction and effect size estimates consistent with those previously reported between MD and drugs for GORD and PUD (OR: 1.23, P = 4.0E−6)⁸⁸. However, the reverse MR analysis (the effect of PG₊M on MD) is also significant (OR: 1.26, P = 2.7E−9). The bidirectional statistically significant results from GSMR usually include two interpretations, one is that there is the bidirectional causality and the other one is horizontal pleiotropy, including an indirect relationship through an intermediate endophenotype or confounder^69,89. Hence, our conclusions must include these interpretations, although we note that the MR Egger intercept test suggests no horizontal pleiotropy (Supplementary Table 22), and the results were similar after mtCOJO conditional analyses, which take into account correlated traits such as educational attainment, BMI, and smoking. In terms of bidirectional causality, there are several possible explanations. Patients with psychological comorbidity often perceive low intensity esophageal stimulation as being painful due to hypervigilance to these intra-esophageal events⁸⁵. Psychological factors can decrease the pressure of the lower esophageal sphincter and change esophageal motility⁹⁰. The reflux symptom itself could result in depression through potentially disabling effects on occupational or social function, or if patients are constantly feeling upset about their condition⁹⁰. Use of medications could conceivably mediate bidirectional associations between depression and PG₊M. Tricyclic antidepressants can lead to a decrease in lower esophageal sphincter pressure, and thus and increase in the number of reflux episodes (anticholinergic effect)⁹¹. Recently, a study shows that use of proton-pump inhibitors (PPI) for acid-related disorders are associated with the subsequent risk of major depression disorder⁹². Further studies are needed to clarify this association. We note a recent population based study of 83.9 million person-years⁹³ found that a diagnosis of mood disorder was associated with increased risk (hazard ratio 1.57) of subsequent peptic ulcer disease/gastritis. While a causal relationship cannot be confirmed between major depression and digestion related disorders (or vice versa), consideration of clinical implications of a possible relationship is justified. When treating patients with MD, awareness of the digestion symptoms for PUD/GORD could help to decide if further interventions are needed. Also, these results may provide clues for screening psychological factors in PUD/GORD patients. A previous study⁹⁴ shows that GORD patients who are also comorbid with psychological distress have more severe symptoms at baseline and more residual symptoms after PPI treatment. In those patients, treatment for the underlying psychological distress might improve the PPI response. A notable negative result was the absence of evidence for any causal relationship between IBD and major depression (an OR of 1.08 (P = 0.4) when major depression was the exposure trait and IBD as outcome and 1.01 (P = 0.11) for the reverse direction). A causal relationship might be expected from observational studies and mechanistic theory linking inflammation with depression⁹⁵, suggesting that observed phenotypic associations may be driven by residual confounding, and having potential implications for future research into the role of inflammation in the pathogenesis of depression.

Despite these interesting findings, our study has several limitations. First, the phenotypes of PUD, GORD, and IBS were combinations of self-reported illness, primary care, hospital admission record, and/or medication use. There is a potential influence regarding accuracy of self-report and misdiagnosis, however, the correlation in effect sizes from case-ness defined by hospital admission, primary care, or self-report data support the use of self-report data, Supplementary Fig. 10e). Given the existent co-reporting of some diagnoses, we conducted sensitivity analyses in which individuals recorded with more than one diagnosis were excluded, but these analyses did not impact our conclusions (Supplementary Tables 6–7, and 11). Formal diagnosis of PUD requires an endoscopy or a upper gastrointestinal barium test, moreover non-ulcer dyspepsia (NUD), not considered here, can be misdiagnosed as PUD. This could influence the downstream MR analyses with depression. Future large datasets of formal-diagnosed PUD are needed to further explore the MR analyses with major depression. Importantly, we note that the GWAS results from IBD derived from the UKB were highly consistent with results from published GWAS, despite a higher prevalence in the UKB than expected in a random population sample. We specifically focus on the broad definition of our phenotypes, and hence heterogeneity in the phenotypes has not been explored. Second, our study is conducted in the UKB cohort, which although a large population study, has recognized volunteer bias⁹⁶. Third, we do not have replication data sets for GORD and PG₊M genome-wide significant SNPs. Fourth, we do not have the H. pylori infection status and microbiome data for all the UKB participants, thus additional analyses on these factors cannot be further explored.

In summary, we identified 35 independent SNPs associated with different digestion disorders, of which 19 SNPs are previously unreported, including eight SNPs at or near MUC1, MUC6, FUT2, PSCA, ABO, CDX2, GAST, and CCKBR genes associated with peptic ulcer disease. Previously established roles of these genes in H. pylori infection, response to counteract infection-related damage and gastric secretion support their involvement. Post-GWAS analyses highlighted the link between the nervous system and the gastrointestinal tract. In addition, MR analyses imply a bi-directional relationship (the risk of GORD, PUD, and taking corresponding medications in liability to major depression and vice versa), which may reflect bidirectional causality, or pleiotropic effects between them. Taken together, our findings demonstrate the role of genetic variants in the etiology of common digestion disorders and the link between depression and PUD/GORD.

Methods

UK Biobank genotyping and quality control

The United Kingdom Biobank (UKB) cohort is a population-based volunteer longitudinal cohort of ~500,000 individuals recruited at 22 centers across the United Kingdom⁹⁷. Genotype data from these individuals were imputed using the Haplotype Reference Consortium (HRC) and UK10K as the reference sample. A European ancestry subset (456,327 individuals, including 348,441 unrelated individuals) was identified by projecting the UKB participants onto the 1000 Genome Project principal components coordinates⁹⁸. Genotype probabilities were converted to hard-call genotypes using PLINK2⁹⁹ (hard-call 0.1) and single nucleotide polymorphisms (SNPs) with minor allele count < 5, Hardy–Weinberg equilibrium test P value <1.0E−5, missing genotype rate > 0.05, or imputation accuracy (Info) score < 0.3 were excluded.

Phenotype definition

The UKB phenotypes used in analyses were derived from two categories: one is disease-diagnoses phenotypes from either a death register, self-reported, hospital admission, or primary care record for the corresponding disease and the other is treatment phenotypes based on the operation and medication-taking code. Supplementary Table 1 lists the corresponding case definition and number of cases for each phenotype. PUD disease-diagnoses cases are a combination of gastric ulcer cases (UKB data field: 131591), duodenal ulcer cases (UKB data field: 131593), other site peptic ulcer cases (UKB data field: 131595) and gastro-jejunal ulcer cases (UKB data field: 131597). The remaining individuals were PUD controls. There were 16,666 cases and 439,661 controls for PUD. For the GORD disease-diagnoses phenotype (UKB data field: 131585), there are 54,854 cases and 401,473 controls. Inflammatory bowel disease (IBD)¹⁰⁰ disease-diagnoses phenotype cases are a combination of Crohn’s diseases (UKB data field: 131627), ulcerative colitis (UKB data field: 131629) diagnoses, giving a total of 7045 cases and 449,282 controls. There were 29,524 cases and 426,803 controls for irritable bowel syndrome (IBS, UKB data field: 131639). For the comorbidity, full-sibling risk and heritability estimation analyses, we used this IBS phenotype definition. For the IBS GWAS analyses, we first removed 1006 participants who also had IBD diagnoses from the original IBS cases, the remaining 28,518 participants were assigned to case status for the IBS disease-diagnoses phenotype while the other 426,803 participants were coded as IBS controls. In clinical practice, medications for PUD also have a therapeutic effect on GORD. Thus, we combined 55,865 individuals taking medications that are mainly considered medications for GORD/PUD (UKB data field: 6154 and 20003), 1162 participants with endorsement for anti-reflux operations (UKB data field: 41272) and diseases-diagnoses cases for GORD and PUD, leaving a total of 90,175 cases (41,945 male cases) and 366,152 controls (phenotype abbreviation: PG₊M—for GORD, PUD, and corresponding medications and treatment). Since the medications for IBD and IBS from UKB can also be used to treat other diseases or relieve symptoms (i.e., medications are not specific), we did not incorporate the medication data for the IBD and IBS phenotypes. The Supplementary Data 1 of Wu et al.⁸⁸ provides UKB medication classification based on Anatomical Therapeutic Chemical (ATC) Classification System¹⁰¹, and we extracted medications for GORD/PUD (the first two ATC level: A02, Supplementary Table 2). For the validation of UKB PUD and IBS genome-wide significant SNPs, we used Genetic Epidemiology Research on Aging (GERA)³⁶ PUD and IBS summary statistics from Zhu et al.³⁷. The definition of PUD and IBS phenotypes from the GERA cohort is described in The Supplementary Table 4 of Zhu et al.³⁷. Briefly, the total 61,847 individuals were divided into two groups: 1004 peptic ulcer cases according to ICD9 code (531–534) and 60,843 controls, respectively. There are 3359 cases and 58,488 controls for IBS phenotypes using the ICD9 code (564). Sensitivity phenotypes for each of PUD, GORD, IBS, and IBD were regenerated by excluding individuals with more than one GI diagnosis among these four disorders from the corresponding original phenotype cases (Supplementary Note 3). To understand the genetic similarity and difference among the self-report, primary care and hospital admission data for each of the PUD, GORD, and IBS, we further divided the cases of each of the PUD, GORD, and IBS into three groups according to the UKB coding (3rd column of Supplementary Tables 1 and 30: primary care only, 40: hospital admission data only, and 50: self-report only) and generated subgroup phenotypes using those subgroup cases, together with the controls from each of the original phenotypes.

Comorbidity analyses

Comorbidity analyses, including comorbidity amongst PUD, GORD, IBS, and IBD and comorbidity between each with depression phenotypes in the UKB, were conducted in 348,441 unrelated individuals. Among the four digestive diseases, for each two of the PUD, GORD, IBS, and IBD cases (six pairs in total), we first checked whether the number of overlapped individuals between case groups is statistically significantly larger than the overlap expected by chance. For each of the PUD, GORD, IBS, and IBD disease (defined as the index disease), we conducted competitive comorbidity analyses to test among the other three diseases which disease is more prone to be comorbid with the index disease. Briefly, we calculated the proportion of the index disease cases in the other three diseases respectively and compared them in pairs by using two-proportions Z-test. One prerequisite for the two-proportion Z-test is that the two samples are independent of each other. Given this, we removed the overlapped cases among these three diseases when calculating the proportion of the index disease cases. For comorbidity between each of the PUD, GORD, IBS, and IBD with depression phenotypes, we first derived eight depression phenotypes based on different data field (UKB data field: 20002, 20216, 2090, 20440, 20442, 2100, 41202, and 41204) and mental health online follow-up (data category: 138) from the UKB. The details for depression phenotype definitions were described Cai et al.⁶⁸. Briefly, depression phenotypes were defined according to help seeking behavior and symptoms, including seen general practice (GP) for nerves, anxiety, tension or depression (abbreviation: GPpsy), seen psychiatrist for nerves, anxiety, tension, or depression (abbreviation: Psypsy), probable recurrent major depression or single probable major depression episode (abbreviation: DepAll), self-reported depression (abbreviation: SelfRepDep), ICD10 defined depression (abbreviation: ICD10Dep), DSM V clinical guideline defined major depression (abbreviation: LifetimeMDD), major depression recurrence (abbreviation: MDDRecur) and seen GP for depression but no cardinal symptoms (abbreviation: GPNoDep) phenotypes. We then checked the overlap individuals from each of the GORD, PUD, IBS, and IBD with the eight depression phenotypes. For each of 32 digestion-depression phenotype pairs, we tested whether the number of individuals who are both digestion phenotype case and depression phenotype case is statistically significantly different from the expected number. Odds ratio, 95% confidence interval (CI) and corresponding P value for each digestion-depression phenotype pair were calculated using fmsb R package (https://cran.r-project.org/web/packages/fmsb/fmsb.pdf). Bonferroni correction was used to account for multiple testing.

Full-sibling risk and estimation of heritability

To demonstrate the genetic component for the UKB case–control phenotypes of PUD, GORD, IBS, and IBD, we estimated the increased risk of the disorders in full-siblings of those affected (we did not incorporate other relative pairs data given the limited sample size for disease cases), compared to the risk in all UKB individuals (disease risk). As described by Bycroft et al.⁹⁷, 22,665 full-sibling pairs were inferred from the kinship coefficients estimated using KING¹⁰². Risk ratio, 95% CI and corresponding P value were calculated using fmsb R package. After obtaining the full-sibling relative risks, we used liability distribution theory^103,104,105 to estimate the heritability of each trait, under the assumption that the increased risk only reflects shared genetic factors.

Genome-wide association study (GWAS) analyses

We performed case–control GWAS analyses using BOLT-LMM¹⁰⁶ with sex, age and 20 ancestry principal components (PCs) fitted as covariates. 543,919 SNPs generated by linkage disequilibrium (LD) pruning (r² < 0.9) from HapMap3 SNPs were used to control for population structure and polygenic effects, including genetic relatedness between individuals. The effect size (β) from BOLT-LMM on the observed 0–1 scale were transformed to odds ratio (OR) using the following equation¹⁰⁷: ${\mathrm{OR}} = \frac{{\left( {k + \beta \left( {1 - p} \right)} \right. \times \left( {1 - k + \beta p} \right)}}{{\left( {k - \beta p} \right) \times \left( {1 - k - \beta \left( {1 - p} \right)} \right)}}$, where k is the proportion of sample that are cases, and P is the allele frequency in the full UKB European cohort. The standard error (s.e.) for OR were then calculated based on the OR and P value from the initial GWAS using the formula ${\mathrm{s}}.{\mathrm{e}}. = \left| {\frac{{\ln ({\mathrm{OR}})}}{{{\mathrm{{\Phi}}}^{ - 1}\left( {\frac{{\mathrm{P}}}{2}} \right)}}} \right|$. A total of 8,546,065 SNPs with minor allele frequency (MAF) > 0.01 were analysed. Independent trait-associated SNPs were generated using GCTA-COJO^34,35 analyses. The genotype data (8,546,065 SNPs with MAF > 0.01) of 20,000 random sampled unrelated European individuals were used to provide a LD reference. Due to the complexity of major histocompatibility complex (MHC) region (25–34 Mb), only the most significant SNP across that region was reported. Regional visualization plots were produced using LocusZoom¹⁰⁸. The genomic inflation factor (λ_GC) was also reported for each phenotype. We used PUD and IBS GWAS summary statistics³⁷ from the GERA³⁶ cohort for validation look-up of the UKB PUD and IBS genome-wide significant SNPs, respectively. We applied RPower package¹⁰⁹ to calculate power for UKB PUD and IBS genome-wide significant SNPs give GERA sample size for each trait. We also conducted pleiotropy (SNP associated with multiple traits) analysis. Briefly, we downloaded published GWAS associations from the GWAS Catalog³⁹ on April 6th 2020. For each of PUD, GORD, PG₊M, and IBS associated SNPs in our study (index SNP), we first selected SNPs from the GWAS Catalog within a ±1000 kb window size of the index SNP. We then selected the GWAS Catalog SNPs significantly associated (P < 5.0E−8) with either mental health-related traits or digestive diseases. We reported a pleiotropic association if selected GWAS Catalog SNPs are in LD (r² > 0.1) with the index SNP. Similarly, for IBD-associated SNPs in our study, we checked whether significant association (P < 5.0E−8) have been reported for inflammatory bowel diseases (including the subtypes) using the downloaded GWAS Catalog data. GWASs for sensitivity phenotypes of each of PUD, GORD, IBS, and IBD were conducted. We also conducted GWASs for the three subgroup phenotypes of each of PUD, GORD, and IBS for further post-GWAS analyses, following the same procedures mentioned above. Given the sex specific effect of rs10512344 on IBS, we further conducted sex-specific GWAS on IBS. To explore the specificity of PUD-associated SNPs, we conducted a multitrait conditional and joint (mtCOJO)³⁷ analysis to condition the PUD GWAS results on GORD GWAS summary statistics and check whether the original PUD-associated SNPs remain significant. We used the GENE2FUNC of FUMA pipeline³⁸, together with GTEx 8th version gene expression data⁶⁰, to investigate whether the genes surrounding PUD-associated SNPs are overexpressed in any specific tissue.

SNP-based heritability and genetic correlations of the five digestion phenotypes

Linkage disequilibrium score regression (LDSC)⁴⁰ was used to estimate SNP-based heritability ($h_{{\rm{SNP}}}^2$) from the GWAS summary statistics. The $h_{{\rm{SNP}}}^2$ estimated on the observed scale were transformed to the liability scale taking the sample lifetime risk (proportion of sample that are cases) as the disease lifetime risk estimates. The summary statistics for each phenotype were filtered using the LDSC default file, w_hm3.snplist, with the default LD scores computed using 1000 Genomes European data (eur_w_ld_chr) as a reference. Genetic correlations (r_g) for each two of the five UKB digestion phenotypes or between each of the five phenotypes and the 258 GWAS traits in LD Hub⁴¹ were calculated using bivariate LDSC¹¹⁰. Since we had a particular interest in the relationship between the five digestion phenotypes and disorders of nervous system, we also used more recent GWAS summary statistics than those included in LDHub for seven psychiatric traits (attention deficit/hyperactivity disorder (ADHD)⁴², schizophrenia (SCZ)⁴³, anxiety disorder⁴⁴, posttraumatic stress disorder (PTSD)⁴⁵, bipolar disorder (BIP)⁴⁶ and autism spectrum disorder (ASD)⁴⁷, major depression (MD)⁴⁸) and two neurological traits (Alzheimer’s disease based on family history⁴⁹ and Parkinson’s disease⁵⁰). As sensitivity analyses, we repeated analyses using sensitivity phenotypes (excluding any individual recorded to have more than one disorder for each of PUD, GORD, IBS, and IBD). Although removal of these individuals could make estimates of SNP-based heritability difficult to interpret, genetic correlations are more robust to such ascertainment⁸⁹. We also used three subgroup phenotypes for each of PUD, GORD, and IBS, as mentioned above, to calculate SNP-based heritibilities and within-diseases genetic correlations. Note, as shown by theory and simulation, estimates of genetic correlations are unbiased even in the context of sample overlap¹¹⁰. In the absence of sample overlap, the genetic covariance intercept term in the bivariate LDSC regression equation is expected to be zero. If sample overlap is present then the intercept term is expected to be a function of the phenotypic correlation, and the proportion of the samples that overlap between the two data sets.

Linking GWAS findings to gene expression

Following the estimation of $h_{{\rm{SNP}}}^2$, we partitioned the SNP-based heritability by genomic features¹¹¹. Briefly, in this method¹¹¹, genetic variants are assigned into 53 functional categories using 24 publicly available annotation data sets, such as UCSC coding, UTRs, promoter and intronic regions¹¹², conserved regions¹¹³ and functional genomic annotations constructed using ENCODE¹¹⁴ and Roadmap Epigenomics Consortium⁶⁶ data. The method evaluates the contribution of each functional category to the overall $h_{{\rm{SNP}}}^2$ of a trait. A category is enriched for $h_{{\rm{SNP}}}^2$ if the variants with high LD to that category have elevated χ² statistics, compared to the expectation given the number of SNPs in the category. In another analysis, genetic variants were annotated to histone marks (H3K4me1, H3K4me3, H3K9ac, and H3K27ac) by cell type specific classes and these annotations were allocated to ten groups: adrenal and pancreas, central nervous system (CNS), cardiovascular, connective and bone, gastrointestinal, immune and hematopoietic, kidney, liver, skeletal muscle, and other. We tested the enrichment of $h_{{\rm{SNP}}}^2$ in tissues relevant to the five digestion phenotypes: the adrenal and pancreas, gastrointestinal, immune and hematopoietic, and liver cell types. We also considered the CNS given the high r_g between four of the six digestion phenotypes and depressive symptoms. We also used LDSC specific expressed genes (SEG)⁵⁹ analysis to test the enrichment of $h_{{\rm{SNP}}}^2$ through gene expression derived cell-type specific annotations. First, given the strong contribution to PG₊M $h_{{\rm{SNP}}}^2$ from the CNS and to IBD $h_{{\rm{SNP}}}^2$ from the immune cell group, LDSC-SEG⁵⁹ was applied to test the enrichment of $h_{{\rm{SNP}}}^2$ in 205 different tissues (53 from GTEx and 152 from Franke lab). Second, given the observation that GORD, PG₊M, and IBS $h_{{\rm{SNP}}}^2$ were enriched in the CNS, we also applied LDSC SEG to test the enrichment of $h_{{\rm{SNP}}}^2$ for GORD, PG₊M, and IBS in 13 brain regions using the multiple brain regions available in the GTEx study (https://gtexportal.org/home) data⁶⁰ to identify specific brain regions implicated by the GWAS results for the three phenotypes. Bonferroni correction was used to account for multiple testing. As sensitivity analyses, we conducted a multitrait conditional and joint (mtCOJO)³⁷ analysis to condition the PG₊M GWAS results on educational attainment (EA)⁶², body mass index (BMI)⁵⁴, and tobacco-use⁶³ GWAS summary statistics and further conducted LDSC cell type group and LDSC-SEG enrichment analyses using gene expression data from the tissues which are statistically significantly enriched above. We checked eQTL (expression quantitative trait loci, i.e., SNPs significantly (P < 5 × 10⁻⁸) associated with gene expression) status for each genome-wide significant SNP using GTEx⁶⁰ results for gastrointestinal tissue and brain tissues. Summary-data-based Mendelian randomization (SMR)⁶⁴ was used to provide evidence for likely causal relationship between the trait-associated SNPs and gene expression. We used eQTLGen¹¹⁵ whole blood eQTL data since this is the largest eQTL data set and many eQTLs are shared across tissues¹¹⁶. To capture more tissue specific eQTL, we used GTEx⁶⁰ eQTL data from six tissues: esophagus-gastroesophageal junction, esophagus mucosa, esophagus muscularis, stomach, small intestine terminal ileum, colon sigmoid, and colon transverse tissue. We also used GTEx eQTL data from hippocampus, brain anterior cingulate (Brodmann area 24, BA24), frontal cortex (Brodmann area 9, BA9) tissue for PG₊M given $h_{{\rm{SNP}}}^2$ enrichment results. The Bonferroni corrected significance threshold was 0.05/131,295, where 131,295 is the number of total genes tested in SMR analyses. Because of its complexity, we do not report results of the MHC region (25–34 Mb)⁴⁸. The statistical framework of SMR⁶⁴ can be applied to mQTL (genome-wide significant SNP-methylation association) data to identify putative methylation-trait association. Hence, we repeated SMR analyses for PUD using blood mQTL data from McRae et al.⁶⁵.

Gene-based and gene-set enrichment analyses

MAGMA (v1.06)⁶⁷ (Multi-marker Analysis of GenoMic Annotation) was used to test for gene-based association based on the SNP association results of the five digestion phenotypes. Gene length boundaries were defined as 35 kilobase (kb) upstream and 10 kb downstream from start and stop site, respectively, to include regulatory elements. The NCBI 37.3 build was used to assign the genetic variants to each gene. SNPs with MAF > 0.01 from 20,000 randomly sampled and unrelated UKB European-ancestry individuals were used to provide a LD reference. A total of 18,402 genes were assessed for an association with each of the six digestion phenotypes with Bonferroni correction used to determine significance (α = 0.05/18402, P < 2.7E−6). We used the results obtained from gene-based analysis, together with gene ontology sets (c5.cc) from MSigDB^117,118 to conduct gene-set enrichment analyses. Competitive test P value for each gene set, as implemented in MAGMA, were computed taking gene size, density, minor allele count and gene–gene correlation into consideration⁶⁷. False discovery rate (FDR)-adjusted P values for biological pathways for each of the five digestion phenotypes were generated using Benjamini and Hochberg’s method¹¹⁹ to account for multiple testing.

Mendelian randomization

We applied the generalized summary-data-based Mendelian randomization (GSMR)³⁷ method to explore the potentially causal effect of MD as an exposure on the five UKB digestion phenotypes as outcome traits (defined as forward direction). GSMR uses summary-level data to test for causal associations between a putative risk factor (exposure) and an outcome trait. Independent genome-wide significant SNPs from the MD GWAS (excluding UKB cohort)⁴⁸ were used as the MR exposure instrument variables. The HEIDI outlier test³⁷ was used to remove outlier pleiotropic genetic instruments associated with both exposure phenotype and outcome phenotype from the analysis. We also conducted reverse causation analysis (i.e., testing the opposite hypothesis that the five digestion phenotypes cause MD). However, GSMR guidelines advise the use of at least ten independent lead SNPs as genetic instruments to achieve robust results. In order to test the effect of PUD (eight SNPs with P < 5.0E−8), GORD (six SNPs), and IBS (two SNPs) on MD, we relaxed the significance threshold to allow for at least ten SNPs for each of the three phenotypes. For each trait as exposure, we only keep the most significant SNP within MHC region for MR analysis (if applicable). We used 0.01 as a clumping threshold r² in GSMR analyses. The other parameters were set to software defaults. For comparison, we also conducted IVW-MR¹²⁰, MR-Egger¹²¹, weighted median-MR¹²², and MR-PRESSO¹²³ analyses following the STROBE-MR guideline¹²⁴ for the effect between PG₊M and MD. We also used LCV⁶⁹ approach for PG₊M and MD given the fact that this approach could control pleiotropy well. To explore whether the bidirectionally statistically significant GSMR results between PG₊M and MD is mediated by EA⁶², BMI⁵⁵, and tobacco-use⁶³ effect, we further used mtCOJO results to conduct bidirectional GSMR analyses, i.e., taking MD conditional GWAS results (the MD original GWAS results conditioned on EA⁶², BMI⁵⁵, and tobacco-use⁶³ GWAS summary statistics separately or together) as exposure and PG₊M non-conditional GWAS results as outcome and vice versa. As a sensitivity analysis, we further removed the depression cases (the combined cases from eight depression phenotypes in UKB as aforementioned) from the five GI disorders phenotypes and conducted GWASs for the five depression removed sensitivity phenotypes (defined as GI-DepComRMV phenotypes). We further repeated bivariate LDSC analyses for genetic correlation and GSMR analyses for MR between major depression and GI-DepComRMV phenotypes. In addition to this, we removed the depression cases of each of eight depression phenotypes in UKB from PUD phenotypes and generated eight PUD sensitivity phenotypes. We further repeated GSMR analyses for MR between major depression and eight PUD sensitivity phenotypes.

Out-of-sample polygenic score prediction

We used MD GWAS summary statistics⁴⁸ (European ancestry, excluding UKB cohort) as discovery data to predict PG₊M risk (risk for PUD, GORD, and likelihood for taking PUD, GORD medications or treatments). The MD data SNPs were matched with the PG₊M SNPs, then LD pruned and “clumped”, discarding variants within 1000 kb of, and in r² ≥ 0.1 with, another (more significant) marker using SNPs with MAF > 0.01 from 20,000 random sampled unrelated UKB European-ancestry individuals as the LD reference. Polygenic score of PG₊M sample individuals were generated for a range of MD GWAS summary statistics data association P value thresholds (5.0E−8, 1.0E−5, 1.0E−4, 1.0E−3, 1.0E−2, 0.05, 0.1, 0.5). For each discovery-target pair, three outcome variables were calculated. (1) The P value of case-control polygenic score difference from logistic regression. (2) Area under the receiver operator characteristic curve using R package pROC¹²⁵, which can be interpreted as the probability of ranking a randomly chosen case higher than a randomly chosen control. (3) Odds ratio and 95% confidence interval for the 2nd–10th polygenic score deciles group compared with 1st decile. The polygenic score decile presented in Fig. 6c is based on association P value thresholds 1.0E−2. Given the availability of depression and psychiatric disorders (not GORD) phenotypes in GERA cohort, we used UKB PG₊M GWAS summary statistics to calculate polygenic score for GERA individuals based on P value threshold 0.05 and conducted out-of-sample polygenic score prediction for depression and psychiatric disorders risk in GERA individuals, respectively. We further conducted out-of-sample prediction for PUD and IBS risk in GERA cohort following the same steps but the P value threshold for polygenic score calculation are 5E−8 and 0.1, respectively.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.