Observational and genetic analyses clarify the relationship between type 2 diabetes mellitus and gallstone disease

Background The relationship between type 2 diabetes mellitus (T2DM) and gallstone disease (GSD) have been incompletely understood. We aimed to investigate their phenotypic and genetic associations and evaluate the biological mechanisms underlying these associations. Methods We first evaluated the phenotypic association between T2DM and GSD using data from the UK Biobank (n>450,000) using a prospective observational design. We then conducted genetic analyses using summary statistics from a meta-analysis of genome-wide association studies of T2DM, with and without adjusting for body mass index (BMI) (Ncase=74,124, Ncontrol=824,006; T2DMadjBMI: Ncase=50,409, Ncontrol=523,897) and GSD (Ncase=43,639, Ncontrol=506,798). Results A unidirectional phenotypic association was observed, where individuals with T2DM exhibited a higher GSD risk (hazard ratio (HR)=1.39, P<0.001), but not in the reverse direction (GSD→T2DM: HR=1.00, P=0.912). The positive T2DM-GSD genetic correlation (rg =0.35, P=7.71×10-23) remained even after adjusting for BMI (T2DMadjBMI: rg =0.22, P=4.48×10-10). Mendelian randomization analyses provided evidence of a unidirectional causal relationship (T2DM→GSD: odds ratio (OR)=1.08, P=4.6×10-8; GSD→T2DM: OR=1.02, P=0.48), even after adjusting for important metabolic confounders (OR=1.02, P=0.02). This association was further corroborated through a comprehensive functional analysis reflected by 23 pleiotropic single nucleotide polymorphisms, as well as multiple neural and motor-enriched tissues. Conclusion Through comprehensive observational and genetic analyses, our study clarified the causal relationship between T2DM and GSD, but not in the reverse direction. These findings might provide new insights into prevention and treatment strategies for T2DM and GSD.


Introduction
Both type 2 diabetes mellitus (T2DM) and gallstone disease (GSD) are prevalent and expensive global public health issues (1,2).The coexistence of these two conditions, known as multimorbidity, poses complex challenges for clinical management (3).For example, while laparoscopic cholecystectomy is the gold standard treatment for GSD (4), individuals with diabetes are usually high-risk candidates for any surgery.Therefore, it is imperative to comprehend the association between T2DM and GSD, as well as the underlying biological mechanisms, to ensure effective prevention and management of this multimorbidity.
The association between T2DM and GSD has been studied, but results have been inconsistent (5,6).A Meta-analysis of prospective cohort studies found no association between diabetes and GSD risk (7), while two latest large-scale prospective cohorts reported a 31%-87% increased risk of GSD in individual with diabetes (5, 8).Conversely, three prospective cohort studies shown that GSD increased the risk of T2DM by 17%-42% (6,9,10).However, these observational studies are prone to bias, confounding, and reverse causality, making it difficult to establish a causal relationship.Mendelian randomization (MR) studies, which address confounding factors and reverse causation (11) have been used to evaluate the causal relationship between T2DM and GSD.Three MR studies consistently suggest that T2DM causally increases GSD risk (12)(13)(14), one study explored and refuted a reverse causal association (6).Despite the knowledge gained from exiting MR analyses advancing our understanding of causal relationships underlying T2DM and GSD, a few major gaps remain.Firstly, the only previous MR suggesting that GSD does not cause T2DM may be inaccurate due to its insufficient statistical power, using merely eight single nucleotide polymorphisms (SNPs) (15).Additionally, multivariable MR (MVMR) can be used to control for pleiotropic impact (16).However, previous MR studies either controlled for no potential confounders (6,14) or only controlled for body mass index (BMI) (12,13).
Despite the increasing number of studies reporting evidence between T2DM and GSD, the biological mechanisms linking these two conditions remain unclear.Genome-wide cross-trait analysis presents an opportunity to comprehensively characterize the shared genetic architectures across traits, shedding light on the underlying biological mechanisms of complex diseases (17).Moreover, both T2DM and GSD are moderately heritable, with SNP-heritability estimates of 25% for GSD (18) and 25%-69% for T2DM (19).Largescale genome-wide association studies (GWASs) have identified a number of disease-associated variants for GSD (20) (N=62) and T2DM (21) (N=386), of which, several risk loci, including PNPLA3, are shared by both conditions (22).However, to the best of our knowledge, no genome-wide cross-trait analysis has been performed to systematically investigate the underlying shared genetic architectures of T2DM and GSD.
Therefore, we conducted a comprehensive bidirectional analysis to investigate the phenotypic and genetic associations between T2DM and GSD.Furthermore, we sought to identify the shared genetic architecture between T2DM and GSD in order to elucidate the underlying biological mechanisms.The study design is outlined in Figure 1.

Data sources
Our observational analyses used data from UK Biobank (UKB).UK Biobank (UKB) is a large prospective cohort study that enrolled over 500,000 participants aged 40-69 years from England, Wales, and Scotland between 2006 and 2010.The study was approved by the National Health Service North West Multi-Centre Research Ethics Committee and all participants provided written informed consent.The diagnoses of T2DM and GSD were based on the International Classification of Diseases, Tenth Revision (ICD-10) and Ninth Revision (ICD-9), as well as UK-Biobank-specific codes, which are available in Supplementary Table S1.We included only 472,050 white participants and excluded those with a history of events at baseline.Participants who self-reported a history of T2DM/GSD at study enrolment (self-reported non-cancer illness, Data-Field 20002), but did not have an ICD-10 or ICD-9 code for T2DM/GSD were also excluded.The participant selection flow chart is presented in Supplementary Figure S1.

Survival analysis
To investigate the observational association between T2DM and GSD, we first assessed the proportional hazards (PH) assumption using stphtest in Stata.We then employed flexible parametric survival Overview of the study design, including the analysis methods in each phase.GSD, gallstone disease; GWAS, genome-wide association study; IVs, instrumental variables; PH, proportional hazards; T2DM, type 2 diabetes mellitus.models (FPSMs) and reported the results as hazard ratios (HRs) with 95% confidence intervals (CIs).This approach allows us to account for time-dependent effects (29).We utilized restricted cubic splines (with four degrees of freedom) in implementation of FPSMs using stpm2 in Stata.The determination of the optimal number of knots and model was based on minimizing the Akaike's information criterion and Bayesian information criterion (29).Covariates were selected based on existing literature and model selection.
The full model was ultimately implemented using potential confounders extracted from baseline questionnaires, including sex, age, assessment centre, the top 40 genetic principal components, BMI, Townsend deprivation index (TDI), education, smoking, alcohol intake, sleep duration, time spent watching television, physical activity (IPAQ), family history of T2DM, diastolic blood pressure (DBP), systolic blood pressure (SBP), TC, LDL and TG.
To ensure the reliability of our findings, we conducted two sensitivity analyses: 1) excluding participants who experienced events within the first two years after being diagnosed with exposure, to minimize potential reverse causation; and 2) using a flexible parametric competing risk regression model (30), with death as the competing event, to rectify overestimation of the probability for the event of interest occurring over time.All P values were two-sided, with statistical significance set at P<0.05.Survival analyses were performed using Stata version 13 (Stata Corp., College Station, TX).

Genome-wide genetic correlation analysis
Genetic correlation represents the average sharing of genetic effect between two traits that is not influenced by environmental confounders.We used linkage-disequilibrium score regression (LDSC) (31) to estimate the global genetic correlation (r g ) between T2DM and GSD using GWAS summary data.We utilized pre-calculated HapMap3 LD scores computed from ~1.2 million common SNPs in European ancestry, commonly acknowledged as well-imputed.A Bonferroni-corrected P-value (0.025 = 0.05/2) was used to define statistical significance.
To further validate the causal relationship identified in the univariable MR analysis, we conducted MVMR and considered potential confounders such as adult BMI, WHR, FI, TC, TG, HDL, LDL, LPA, ApoA, ApoB, smoking, and alcohol intake based on a review of existing literature.

Fine-mapping credible set and colocalization analysis
To investigate the causal SNP for these index SNPs identified from CPASSOC, we performed a fine-mapping (FM) analysis using Bayesian FM method (36) to estimate credible sets of SNPs that had a 99% likelihood of containing the causal SNPs.We estimated a posterior inclusion probability (PIP) for each index SNP using the steepest descent approximation.
We further conducted colocalization analysis to determine whether shared SNPs identified by CPASSOC were shared causal variant or distinct causal variants using a Bayesian method, Coloc (37).A locus was considered colocalized if the posterior probability for H4 (PPH4) exceeded 0.7.

Transcriptome-wide association study analysis
We conducted a transcriptome-wide association study (TWAS) analysis using FUSION (38) to investigate specific tissue-gene pairs shared by T2DM and GSD.Firstly, we performed single-trait TWAS using the expression weights from 49 post-mortem Genotype-Tissue Expression project tissues.Subsequently, we combined the single-trait TWAS result and determined the shared gene-tissue pairs across T2DM and GSD.Bonferroni correction was used to identify expression-trait associations.

Phenotypic association
The baseline characteristics of UKB participants included in the observational analysis are presented in Supplementary Tables S3 and S4.In terms of the impact of T2DM on GSD risk, 455,405 participants were followed for 5,306,344 person-years (11.7 ± 2.7 years), during which 1,553 participants with T2DM and 14,425 T2DM-free participants developed GSD.The PH assumption test indicated that T2DM and certain covariates, such as BMI, education, and alcohol intake, exhibited a timedependent effect on GSD risk.The FPSM analysis revealed that participants with T2DM had a higher risk of GSD (HR=1.71,95% CI=1.58-1.84,P<0.001), which decreased over time when adjusting for sex, age, assessment center, and the top 40 genetic principal components.Additional adjustment for confounders had minimal impact on the hazard of GSD in the final model (HR=1.39,95% CI=1.29-1.50,P<0.001).Furthermore, the association between T2DM and GSD remained consistent when considering death as a competing event (HR=1.23,95% CI=1.14-1.32,P<0.001), or when excluding participants who developed GSD within the first two years after T2DM diagnosis (HR=1.17,95% CI=1.08-1.28,P<0.001) (Table 1).

A B
Univariable and multivariable mendelian randomization analysis between type 2 diabetes mellitus and gallstone disease.CI, confidence intervals; GSD, gallstone disease; MR, mendelian randomization; No. of, the number of; SNPs, single nucleotide polymorphisms; T2DM, type 2 diabetes mellitus; *, removing the effect of body mass index on type 2 diabetes mellitus.
To investigate the causal SNP for these 23 pleiotropic SNPs, we estimated the 99% credible sets of causal SNPs.A total of 368 potential causal SNPs were obtained from the 23 pleiotropic SNPs.Note that the 99% credible set of rs1800961, rs1260326 and five other pleiotropic SNPs only included themselves.To further distinguish the shared causal SNPs from the distinct ones, we assessed statistical colocalization of these 23 pleiotropic SNPs.The results revealed that 56.2% (containing rs1800961, rs1260326 and 11 other pleiotropic SNPs) of shared loci colocalized at the same candidate causal SNPs.In summary, we found a good number of loci shared by T2DM and GSD, and in particular identified rs1800961 and rs1260326 as potential shared causal variant for T2DM and GSD.
After removing the effect of BMI on T2DM, we identified five additional pleiotropic SNPs, in addition to the ten pleiotropic SNPs also shared by T2DM and GSD.Notably, rs1800961 (P CPASSOC =1.44×10 -59 ) and rs1260326 (P CPASSOC =6.81×10 -40 ) were also identified as the most significant pleiotropic SNPs and shared causal SNPs for T2DM and GSD even after adjusting for BMI.

Discussion
To our knowledge, this study represents the most comprehensive research on the bidirectional relationships between T2DM and GSD, combining observational and genetic analysis.Our study identified a unidirectional causality running from T2DM to the risk of GSD.Furthermore, we identified specific genetic variants, such as rs1800961 (HNF4A) and rs1260326 (GCKR), that contribute to the biological links between T2DM and GSD.These findings advance our understanding of the complicated relationship underlying T2DM and GSD, providing important implications for preventing and treating these common diseases.
Expanding upon prior research, our study takes a comprehensive approach to investigate the association between T2DM and GSD risk through both observational and MR studies.Our findings, obtained through survival and MR analyses, consistently demonstrate a causal relationship between T2DM and GSD risk.These results align with the updated meta-analysis of prospective cohort studies (RR=1.49,95% CI=1.09-2.03,P=0.012, Supplementary Figure 2), which pooled data from previous meta-analysis (7) and two recent large-scale studies (5, 8).While in line with the findings of three existing MR studies (12)(13)(14), TABLE 2 The cross-trait meta-analysis, fine-mapping analysis and colocalization analysis for type 2 diabetes mellitus and gallstone disease ‡ .our MR expands previous results in two critical aspects.First, we minimize the influence of reverse causality by applying a bidirectional design.Second, we limit the impact of confounders by adjusting for potential confounding factors using the MVMR.Consequently, we assert that a causal relationship exists between T2DM and GSD risk.

SNPs
Regarding the association between GSD and T2DM risk, previous studies have reported conflicting results.Some prospective cohort studies (6, 9, 10) have shown a positive association, while one MR study fund no association (6).Our study used both survival and MR analyses and found no causal association between GSD and T2DM risk.The discrepancy between our study and previous cohort studies may be attributed to two aspects.Firstly, two of these previous cohort studies (9, 10) relied on self-reported data, which may have led to misclassification or underestimation of cases.Secondly, observational studies are prone to reverse causality, which the previous cohort study did not account for (6).However, after considering a 2-year lag time, our study found no association between GSD and T2DM.This implies that previous significant findings may have been distorted by reverse causality.MR studies are commonly employed to address reverse causality (11).Our MR analysis replicates the null findings of the previous MR study ( 6), yet with a much larger sample size (43,639 vs. 1,110 T2DM cases) and more IVs (62 vs. 8).Additionally, our study considered the significant contribution of BMI to T2DM development, which was not accounted for in the previous MR study.
Through cross-trait and colocalization analyses, we identified shared biological mechanisms.Specifically, the key SNPs associated with both conditions were rs1800961 mapped HNF4A and rs1260326 mapped GCKR.HNF4A regulates liver-specific gene expression involved in lipid transport, glucose, and bile metabolism (39).It is also essential for insulin secretion in pancreatic beta cells (40).Additionally, HNF4A has been linked to elevated levels of gamma-glutamyl transferase (GGT), a sensitive marker of cholestasis (41).The effects of HNF4A on insulin action and GGT contribute to the development of both T2DM and GSD.On the other hand, GCKR regulates glucose conversion to glucose-6-phosphate in the liver and pancreas (42, 43).It is associated with various metabolites involved in carbohydrate and lipid metabolism (42).The minor allele of GCKR has been associated with hepatospecific glucokinase activation, reduced plasma insulin levels, and protection against T2DM (42).Additionally, GCKR enhances hepatic cholesterol availability, leading to elevated bile cholesterol concentration and GSD development (43).
Our findings deliver important clinical and public health implications.First, our found that T2DM causally contributes to the development of GSD, but GSD does not increase T2DM risk.This suggests that preventing and treating GSD may not significantly reduce T2DM risk.However, our findings highlight the important of focusing on prevention and treatment of T2DM.Second, while new drugs for T2DM have been developed, further research is needed to explore novel treatments (44).It is noteworthy that certain medications commonly prescribed for T2DM have been associated with an increased risk of GSD in large meta-analyses of randomized controlled trials (45).Furthermore, our genetic analyses reveal shared genetic architectures between T2DM and  CPA, Cross-Phenotype Association; CumSum, cumulative sum of posterior inclusion probability; GSD, gallstone disease; SNPs, single nucleotide Polymorphisms; T2DM, type 2 diabetes mellitus; T2DM adj BMI, removing the effect of body mass index on type 2 diabetes mellitus; ‡, PCPASSOC < 5×10-8, single trait P-value < 1×10-5, clumping r2 = 0.2; *, a novel SNP (A significant index SNP met the following criteria was identified as novel index SNP: 1) did not reach genome-wide significance (5×10 -8 < P single-trait < 1×10 -5 ) in single-trait GWAS; 2) was not in LD (r 2 < 0.5) with any previously reported genome-wide significant SNPs in single-trait GWAS, and none of their neighboring SNPs (± 500 kb) reached P < 5×10 -8   in single-trait GWAS.).GSD, which can enhance our understanding of biological mechanisms and potentially identify therapeutic targets for T2DM.For example, our findings suggest that HNF4A and GCKR may be promising targets for T2DM therapies (42), but their effects on GSD and other comorbidities should be considered.Several limitations should be considered when interpreting the results and conclusions of this study.Firstly, fasting insulin may confound the association between T2DM and GSD risk, which we were unable to adjust for in our multivariable survival analysis due to a lack of information from prospective cohort.However, after adjusting for FI in our MVMR, we confirmed the causal relationship between T2DM and GSD risk.Additionally, dietary behaviors such as egg, vegetable, and fruit intake may act as other confounders for the association between T2DM and GSD.Unfortunately, due to a larger proportion (>50%) of missing data in the prospective cohort, we were unable to adjust for these factors in our multivariable survival analysis.Furthermore, it is worth noting that half of the genetic variants associated with dietary behaviors are a consequence of increased BMI (46).Considering the potential collinearity between BMI and dietary behaviors, we chose not to include the latter in the MVMR analysis.Secondly, our findings were restricted to European population to control for population stratification.However, this may limit the generalizability of our results.Additionally, despite over 90% of gallstones in Europe are cholesterol gallstones, there are also pigment gallstones with different etiologies, which were not specifically investigated in our study.

Conclusions
Our study demonstrates an association between T2DM and GSD risk at both the phenotypic and genetic levels using large-scale prospective cohort and GWAS data.This association is independent of BMI, WHR, and FI, suggesting an intrinsic and causal relationship.However, no causal association was found between GSD and T2DM risk.Furthermore, the shared genetic architecture between GSD and T2DM enhance our understanding of the underlying biological mechanisms.These findings might offer valuable insights for the identification of potential therapeutic targets for T2DM and novel perspectives on preventing GSD, ultimately contributing to a decrease in multimorbidity incidence.

TABLE 1
Phenotypic association between type 2 diabetes mellitus and gallstone disease.
CI, confidence interval; GSD, gallstone disease; HR, hazard ratio; T2DM, type 2 diabetes mellitus; Time interaction, time-varying effect of type 2 diabetes mellitus interacting with survival time; No. of, the number of; P interaction , the P value of interaction of T2DM and time; Model 1 was adjusted for sex, age, assessment center, and the top 40 genetic principal components; Model 2 was further adjusted for body mass index, Townsend deprivation index, education, smoking, alcohol intake, sleep duration, time spent watching TV, and physical activity; Model 3 was further adjusted for family history of Type 2 diabetes mellitus, cholesterol, low density lipoprotein and triglycerides; Model 4 was flexible parametric competing risk regression model, and death as the competing event; Model 5 was excluded the participant whose event occurring in the first two years after diagnosis of exposure.

TABLE 3
Shared transcriptome-wide association study significant genes between type 2 diabetes mellitus and gallstone disease.