Shared Genetic Links Between Digestive and Brain Degenerative Diseases

The GWAS data used in this study come from publicly available datasets from IEU and GWAS Catalog, including four GI diseases (IBD (Glanville et al. 2021), IBS (Eijsbouts et al. 2021), diverticular disease (DVD) (Schafmayer et al. 2019), and gastroesophageal reflux disease (GERD) (Ong et al. 2022)) and three ND (ALS (van Rheenen et al. 2021), AD excluding UK Biobank and 23andMe data (Wightman et al. 2021), and PD from the International Parkinson's Disease Genomics Consortium) (Supplementary Table S1). Blood eQTLs were sourced from the eQTLGen Consortium, which encompasses information on 10,317 single nucleotide polymorphisms (SNPs) associated with traits from 31,684 individuals (Võsa et al. 2021). The total sample size of the GWAS data exceeded 100,000, primarily from European populations to minimize ethnic differences. AD data SNPs were imputed using the hg19 genome version from https://hgdownload.soe.ucsc.edu/downloads.html↗. All GWAS analyses in this study were based on the human genome reference sequence GRCh37 (also known as hg19), chosen for its widespread support in genetic databases and analytical tools, ensuring accuracy and consistency in comparative research. In this study, the SNVs involved were excluded from all analyses except for LAVA in the MHC region (chromosome 6: 25–35 MB).

In this study, we applied Linkage Disequilibrium Score Regression (LDSC), High‐Definition Likelihood (HDL), and Local Analysis of Variant Annotation (LAVA) methods to assess the genome‐wide genetic correlations among 12 pairs of traits between four GI diseases and three ND (Werme et al. 2022). Quantifying the contribution of each element by examining the relationship between the test statistic and LDSC allowed determining whether the different traits share a common genetic basis (Bulik‐Sullivan et al. 2015). As a powerful method for estimating heritability and genetic correlations based on GWAS summary statistics, the LDSC can be used to distinguish between true polygenicity and mixed biases (e.g., population stratification and recessive associations) (Xu et al. 2022). The initial analysis phase involves applying cross‐trait LDSC (https://github.com/bulik/ldsc↗) and HDL (Ning et al. 2020) methods to estimate genetic correlations. Compared to LDSC, HDL comprehensively accounts for genome‐wide linkage disequilibrium (LD), thereby significantly improving the precision of the estimation. Local Analysis of Variant Annotation estimated local genetic correlations in 2495 genome‐wide regions. To adjust for multiple comparisons, the Benjamini–Hochberg method was applied across all pairwise comparisons.

Genetic correlations only reveal genome‐wide associations; hence, Genetic Analysis Incorporating Pleiotropy and Annotation (GPA) methods were further employed to explore comprehensive genetic overlaps among different traits (Gong et al. 2023). In GPA analysis, four models (M00, M10, M01, M11) categorized SNPs into four groups, assessing the overall genetic overlap's statistical significance through likelihood ratio tests (LRT). This method allowed a more thorough understanding of genetic interactions between traits and provided richer genetic background information for the study. Significance levels were set at p < 4.17 × 10⁻³ (0.05/12) after Bonferroni correction. p‐values between 0.05 and 4.17 × 10⁻³ were considered suggestive associations, indicating potential associations.

The core goal of GPA is to estimate the proportion of SNPs (PM) in these models and assess the statistical significance of the overall genetic overlap by LRT. This approach allows us to understand genetic interactions among different traits in a more comprehensive way and provides richer information about the genetic background of the study. The significance level was set at p < 4.17 × 10⁻³ (0.05/12) after Bonferroni correction. If the p‐value was between 0.05 and 4.17 × 10⁻³, we regarded it as the presence of a suggestive association.

In cases where trait pairs exhibited noteworthy genetic correlations or overlaps, this study implemented the PLACO approach to pinpoint potential pleiotropic SNVs (Ray and Chatterjee 2020, Ray et al. 2021). In the study, SNVs with a p‐value less than 5 × 10⁻⁸ in the PLACO analysis were identified as significant pleiotropic variants. To further characterize these loci, Functional Mapping and Annotation of Genetic Associations (FUMA) was employed (Watanabe et al. 2017). We used the default settings for FUMA annotation and used data from European pedigrees in 1000 Genomes Project Phase III as a reference panel. SNPs with p < 5 × 10⁻⁸ and independent of each other at r² < 0.6 within 1 Mb were defined as independently significant SNPs. Lead SNPs are a subset of independently significant SNPs that are independent of each other at r² < 0.1. Genomic risk loci are identified by combining LD blocks (< 250 kb) of independent significant SNPs that are closely linked to each other. Top lead SNPs were defined as those with the lowest p‐value in a given region. functional annotations, including ANNOVAR category (Wang et al. 2010). Functional annotations, including ANNOVAR category (Ray et al. 2021), Combinatorial Annotation Dependent Deletion (CADD) (Kircher et al. 2014) scoring, and RegulomeDB (Boyle et al. 2012) score, also obtained by FUMA. Bayesian colocalization analysis was then conducted to further identify shared causal variants within each pleiotropic locus. The analysis was based on the coloc method using default priors (p1 = 1 × 10⁻⁴, p2 = 1 × 10⁻⁴, and p12 = 1 × 10), with loci having a posterior probability for H4 (PP.H4) greater than 0.7 determined as colocalized loci (Giambartolomei et al. 2014).

Expanding on the PLACO findings, our study delved into the shared biological processes of these pleiotropic motifs. This involved conducting a GenoMic Annotation (MAGMA) analysis at the gene level, focusing on genes situated within or intersecting with pleiotropic motifs, guided by both PLACO outcomes and single‐trait GWAS data (de Leeuw et al. 2015). The research applied MAGMA analysis on PLACO results and single‐trait GWAS data to identify potential pleiotropic genes, with significance thresholds set differently for each. For the MAGMA analysis of PLACO data, a Bonferroni‐corrected p‐value of less than 0.05 was deemed significant, whereas for the original single‐trait GWAS in MAGMA analysis, a p‐value of less than 0.05 was used.

For gene set functional enrichment analysis of MAGMA outputs, we (https://www.gsea‐msigdb.org/gsea/msigdb/human/collections.jsp↗) obtained the latest KEGG and GO sets and used the clusterProfiler in RStudio to perform the analysis. Enrichment analysis (Yu et al. 2012) to obtain the results of gene set enrichment. Setting conditions of at least two enriched genes, as well as a p‐value of < 0.05 and an FDR of < 0.05, were considered statistically significant.

The Summary data‐based Mendelian Randomization (SMR) software tool was initially developed to implement the SMR & HEIDI method (Zhu et al. 2016). This method integrates data from GWAS and expression quantitative trait loci (eQTL) studies to assess the pleiotropic association between gene expression levels and complex traits of our interest. In terms of LD calculations, we adopted a method based on the 1000 Genomes European reference data (Auton et al. 2015). P‐values were calculated using the BH method, and results were considered significant with a P fdr < 0.05 and heterogeneity HEIDI p > 0.01 as the inclusion criteria (Zhu et al. 2016).

Finally, we performed two‐way Mendelian randomization analyses of 12 paired features between GI and ND to investigate potential causal features. We used an inverse variance weighting method as the primary analysis for our TSMR analyses using the "TwoSampleMR" package in RStudio (version 4.3.0, https://www.r‐project.org/↗). TSMR analyses (Emdin et al. 2017, Smith and Ebrahim 2003, Haycock et al. 2017) are based on three core assumptions: (1) there is a robust and significant association between Instrumental Variable (IV) and exposure factors; (2) IV is independent of confounders; and (3) IV affects the outcome only by influencing exposure factors and has no other pathways (Davey Smith and Hemani 2014).

To assess the heterogeneity among different Single Nucleotide Polymorphisms (SNP), we used Cochran's Q statistic methodology (Hemani et al. 2018). To verify the robustness of the results, we performed a sensitivity analysis. We also applied the MR‐Egger method to detect the presence of pleiotropy in SNPs, measuring directional pleiotropy by the intercept of MR‐Egger regression (Bowden et al. 2015). We chose a Benjamini–Hochberg p < 0.05 as the significance threshold.

Using summary statistics from publicly available genome‐wide association studies, significant genetic overlaps were identified in 12 paired traits. LDSC revealed suggestive associations in two trait pairs, IBD‐PD and GERD‐PD, with GERD‐PD also showing suggestive association in HDL. Notably, the other ten trait pairs were identified as having significant genetic overlap but not genetic correlation, resulting in a final combined set of 12 paired traits for subsequent analysis. The LDSC results were highly consistent with those from High‐definition Likelihood (Supplementary Tables S2–S5) (Figure 2).

Local genetic correlation analysis revealed 28 significant regions, of which 23 were unique significant regions (P_fdr < 0.05). The correlation coefficients (rho) of six loci were equal to −1 or 1, indicating strong genetic correlations. Notably, two regions between DVD and PD phenotypes exhibited positive and negative genetic associations, respectively, as detailed in Supplementary Tables S6–S7. Both regions between IBD and PD phenotypes also showed positive correlations, as specifically seen in Supplementary Tables S6–S7. Compared to global genetic correlation, local genetic correlation analysis can more accurately capture genetic links showing mixed effect directions. This means that even if a pair of phenotypes are not significant in global genetic correlation analysis, it is possibly due to a rough balance of positive and negative (opposite effect direction) local genetic correlations. Therefore, even in cases where significant correlations were not detected in LDSC and HDL analyses (P_fdr < 0.05), LAVA still identified some significant regions (P_fdr < 0.05). Moreover, although only suggestive associations were shown between IBD and PD and GERD and PD in global genetic correlation analysis, local genetic correlation analysis revealed potential associations between these phenotype pairs in multiple regions (p < 0.05), and in the IBD and PD phenotype pair, some regions showed significant genetic correlations (P_fdr < 0.05), as detailed in Supplementary Tables S6–S7 (Figure 3).

Through the PLACO method, the study identified 1457 SNVs (with 1213 being unique) as potential pleiotropic variants across a dozen trait pairs. Concurrently, FUMA successfully pinpointed 66 independent genomic loci with pleiotropic characteristics, covering 47 distinct chromosomal regions. (Supplementary Tables S8‐S9) (Supplementary Figure S1). Pleiotropic regions identified in multiple paired traits, such as 19q13.32 (mapped genes: APOE, CKM, PVRL2, and PVRL2:CTB‐129P6.4) and 1p31.3 (mapped genes: PDE4B, IL12RB2, IL23R, and C1orf141:IL23R), along with 2q14.3 (mapped gene: BIN1), indicate these sites' broad pleiotropy. Out of the 66 pleiotropic loci identified by FUMA, top SNVs in 58 pairs showed significance in both traits. Twenty‐seven of these pairs (about 47%) demonstrated consistent associations, suggesting that these variants might concurrently influence the risk factors for both GI and mental health conditions. Conversely, the remaining 31 pairs (around 53%) presented inconsistent associations, hinting at the possibility of diverse underlying biological processes (Supplementary Table S10).

The ANNOVAR category annotation revealed that within the 66 indexed SNVs, 33 (50%) were intronic, 15 (23%) intergenic, and only 2 (3%) exonic. Interestingly, these two exonic variants were located at distinct loci and corresponded to two different genes. Specifically, the index SNV rs429358 at 19q13.32 (PLACO p = 7.99 × 10⁻⁵⁷) in APOE, a gene associated with various physiological processes including cholesterol metabolism and Alzheimer's disease risk, was significantly associated with AD (original GWAS p = 0, recalculated p < 2.2 × 10⁻¹⁶) but not with DVD GWAS. The index SNV rs34311866 at 4p16.3 (PLACO p = 6.30 × 10⁻¹¹) in TMEM175 was identified with an exonic variant CADD score > 12.37, rs429358 (CADD score: 12.64; mapped gene: APOE), indicating potential deleteriousness. An intronic SNV, rs344816, identified in the GERD‐AD trait pair with a RegulomeDB (RDB) score of 1b, suggests this variant is likely to affect transcription factor binding and is closely related to the regulation of adjacent gene expression.

Subsequent colocalization analysis of the 66 identified potential pleiotropic loci showed that 9 loci (approximately 14%) had a PP.H4 value greater than 0.7 (Table 1). Within these 9 loci, 7 top SNVs were singled out as likely shared causal variants. Surprisingly, the 3q25.1 locus identified as a pleiotropic locus for two trait pairs, IBS and PD (PP.H4 = 0.9343) as well as DVD and PD (PP.H4 = 0.8808), co‐localized between IBS and PD (Top SNP rs9840232, mapped gene MED12L) and DVD and PD (Top SNP rs17204437, mapped gene MED12L:P2RY12) (Supplementary Table S11) (Supplementary Figure S2).

MAGMA analysis identified 82 significant pleiotropic genes (74 unique) from 185 potential pleiotropic genes located in or overlapping with 66 pleiotropic gene loci. Among these, 8 genes showed pleiotropy in two trait pairs, including CHRNB1, MED12L, P2RY14, GPR87, P2RY12, BCL7C, SNCA, and PVRL2 (Supplementary Table S12) (Figure 4).

The gene enrichment analysis revealed 17 significant pathways (p‐value < 0.05, P_adjust < 0.05). No significant enrichment was observed in KEGG; all pathways were from GO. These pathways mainly relate to signal transduction, amyloid fibril formation and regulation, and lipid metabolism and transport. Notably, four significant pathways, including G protein‐coupled purinergic nucleotide receptor signaling and activity, were identified in both IBS‐PD and DVD‐PD trait pairs. Multiple pathways related to amyloid fibril and protein regulation, as well as cholesterol and lipoprotein particle‐related pathways, were found in the GERD‐AD trait pair (Supplementary Table S13) (Supplementary Figure S3).

We incorporated 74 pleiotropic genes identified by MAGMA analysis into our study and mapped these genes using the eQTLgen gene set, successfully associating them with 57 genes. SMR analysis identified 26 genes with varying correlations to the seven diseases studied (P_fdr < 0.05, HEIDI p > 0.01). Among them, EP300 was related to four diseases. CHRNB1, KNOP1, P2RY14, and POLR2A were related to three diseases. C1orf106, FCF1, HSD3B7, RANGAP1, SCAPER, SETD1A, UBTF, VKORC1, and ZNF668 were related to two diseases (Supplementary Table S16) (Figure 5).

Additionally, this study retrieved drug‐gene interaction data from the Drug‐Gene Interaction database (DGIdb), identifying 8 genes with interactions with multiple drugs, most of which were inhibitory. Specifically, the EP300, PRKCB, and VKORC1 genes demonstrated inhibitory interactions with several drugs. Meanwhile, the CHRNB1 gene showed interactions with both antagonists and agonists (Supplementary Table S17).

Bidirectional TSMR was conducted for 12 pairs of GI and ND. Instrumental variable SNPs (p < 5×10⁻⁸) were selected, and linkage disequilibrium was calculated based on the 1000 Genomes European panel using the PLINK clump method, with an r²< 0.01 (clumping distance = 5000 kb) as the threshold for linkage equilibrium.

In each direction, one significant trait pair was identified. For GERD as the exposure and PD as the outcome, the Mendelian randomization results showed that an increased risk of GERD lowers the risk of PD (P_FDR = 0.0105, P_Q = 0.3542, P_{Egger intercept} = 0.7779), seemingly validating conclusions drawn from LDSC and HDL analyses (Supplementary Table S14). With ALS as the exposure and DVD as the outcome, increased ALS risk was associated with reduced DVD risk (P_FDR = 0.03540, P_Q = 0.8478, P_{Egger intercept} = 0.1356) (Supplementary Table S15) (Supplementary Figure S4).