Causal relationship between serum bilirubin, cholelithiasis and hepatobiliary and pancreatic malignancy: A bidirectional and 2-sample Mendelian randomization study.

1. Introduction
Bilirubin is a metabolite produced during red blood cell destruction and an important intermediate product of heme metabolism. It is metabolized primarily by the liver and excreted through bile. Although bilirubin has traditionally been regarded as a harmful waste product, recent studies have found it has significant antioxidant, anti-inflammatory, and cell-protective effects.[1] These characteristics make serum bilirubin levels the focus of research on a variety of diseases, including cardiovascular diseases, metabolic syndrome, and malignancies.[2]
Cholelithiasis (CHOL) is a common digestive disease caused by stones formed in the gallbladder or bile duct, which is usually divided into 2 categories: cholesterol stones and bilirubin stones.[3] The pathogenesis of cholelithiasis is complex, involving many metabolic abnormalities, especially the disorder of bilirubin metabolism. Studies have shown that patients with cholelithiasis have a significantly increased risk of gallbladder cancer (GDC), liver hepatocellular carcinoma (LIHC), and pancreatic adenocarcinoma (PADC).[4,5] However, most of these studies are based on observational data, and it is not easy to clearly distinguish causation from association, so there are potential confounding factors. Mendelian randomization (MR) is a method that uses genetic variation as an instrumental variable (IV) to study the causal relationship between exposure and outcome.[6] Due to the random distribution of genes, MR can essentially control confounding factors in observational studies and improve the accuracy of causal inference.
This study will use MR to investigate the causal relationship between serum bilirubin levels and cholelithiasis and gallbladder cancer, hepatocellular carcinoma, and pancreatic cancer. Specifically, this study will use genetic variants associated with exposure factors such as serum bilirubin levels and cholelithiasis, known as single nucleotide polymorphisms (SNPs), obtained from genome-wide association studies (GWAS), as IVs to assess their causal effects on 3 digestive malignancies. In this bidirectional 2-sample MR study, we found that genetically predicted serum direct bilirubin was a significant risk factor for gallbladder cancer, while no significant causal associations were observed for other exposure-outcome pairs or in the reverse direction.

2. Methods

2.1. Data sources
All data used in this study were obtained from publicly available online databases. Cholelithiasis, serum total bilirubin, direct bilirubin, and indirect bilirubin were selected as exposure factors, while gallbladder cancer, hepatocellular carcinoma, and pancreatic cancer were chosen as outcomes. The relevant SNPs were extracted from GWAS. This study utilized GWAS datasets through the OpenGWAS platform in R (https://gwas.mrcieu.ac.uk/).[7,8] Specifically, CHOL data were sourced from the FinnGen database (https://www.finngen.fi/en), with a study population of European ancestry comprising 19,883 cases and 1,95,144 controls. Direct bilirubin (DBIL) data were extracted from the GWAS catalog (https://www.ebi.ac.uk/gwas/studies/),[8] covering a European population of 3,72,420 individuals and containing 42,06,076 SNPs. Total bilirubin data (TBIL), also sourced from the GWAS catalog, included 3,42,829 participants of European ancestry. Indirect bilirubin (IBIL) data were retrieved from the UK Biobank (https://pan.ukbb.broadinstitute.org/), focusing on an Asian population with 6972 cases.[8] GDC GWAS data were obtained from the eLife database (https://elifesciences.org/articles/34408), involving an Asian population with 41 cases and 866 controls. Similarly, LIHC data were sourced from the eLife database comprising a European cohort with 168 cases and 3,72,016 controls.[9] Data on PADC were drawn from the GWAS Catalog, with the study population consisting of 499 Asian cases and 1,59,201 controls.[9,10]

2.2. Process
This study adheres to the STROBE-MR guidelines for the reporting of MRstudies, ensuring transparency and reproducibility in our methodology.

2.2.1. MR analysis is based on 3 assumptions
Independence assumption: The selected IV should be strongly correlated with exposure factors.

Independence assumption: The association between the selected IV and the outcome must be realized through exposure and cannot be associated with other confounding factors.

Exclusion restriction assumption: IV can only affect the outcome through exposure but not directly, meaning that the association between gene variation and the outcome cannot have a direct effect through a path other than exposure.[11]

All analysis and mapping in this study were performed using R 4.3.2 (https://www.r-project.org/) and related toolkits.

2.2.2. Select the IV
Set the significance threshold: Since multiple hypothesis tests in the same data set will increase the probability of type I errors, Bonferroni correction was carried out for the significance level (P-value) to maximize the statistical significance of the selected IVs. After correction, we selected P1 = 5e–8 or P1 = 1e–5 as the threshold value to filter the selection of genetic variation SNP.[11]

Removal of chained disequilibrium (LD): Based on the theory of LD, we removed SNPS within the ±1 Mb range of IVs and r² < 0.001.[12]

Calculate F-statistic and weak IV: F statistic is a statistic to measure the strength of correlation between IV and exposure factor. The calculation formula is as follows: or .

IVs with F-statistic < 10 are considered weak IV and removed.[13] Table S7, Supplemental Digital Content, https://links.lww.com/MD/Q510 presents the IVs for the 4 exposure factors.
[4]Allele harmonization: In order to ensure consistency in exposure and outcome analysis, we performed allele harmonization for each SNP, that is, to ensure that the allele direction of the same SNP was consistent in different databases or analyses.[14] Inconsistencies in alleles may lead to misjudgments in the direction of causal effects. That is, if the risk alleles of the same SNP are incorrectly assigned in different directions in the analysis of exposure and outcome, bias may occur (Fig. 1).

2.2.3. MR analysis
[1a].Forward analysis: First, the causal effects of exposure (bilirubin and cholelithiasis) on outcomes (GDC, LIHC, PADC) are analyzed.

[1b].Reverse analysis: Since direct bilirubin showed a causal relationship with gallbladder cancer, a reverse MR analysis was performed to rule out possible reverse causality.

Analysis method: MR-Egger regression, weighted median method, inverse variance weighted (IVW) method, simple mode method, weighted mode method.[15,16]

[2].Sensitivity analysis – through multiple IVs analysis: To evaluate the robustness of multiple IVs, different statistical models are used, including the weighted median method, MR-Egger regression, and model-based estimation.[17]
Heterogeneity tests: Sources of heterogeneity include pleiotropy, weak effects of IVs, measurement errors, gene-environment interactions, and different genetic backgrounds. Heterogeneity test was used to detect the consistency of the effects of IVs. In this study, Q statistics were used to assess heterogeneity among IVs to detect the presence of multiple independent effects that could affect outcomes.

Pleiotropic test: In MR analyses, results may be biased if the association between exposure and outcome variables is due to pleiotropy rather than causation. In this study, the MR-Egger regression method or MR-PRESSO (MR-Pleiotropy RESidual Sum and Outlier) was used to test the pleiotropy.[18]

3. Results

3.1. Causal effects of serum bilirubin and cholelithiasis on hepatobiliary-pancreatic cancers
Based on 2-sample and bidirectional MR, 5 of the most widely used statistical methods were used for each MR analysis in this study. Since the IVW method adopted a random effects model, it could deal with heterogeneity to a large extent, and the analytical results were more reliable than the other 4 methods. Therefore, in MR, researchers generally focus on the results of IVW analysis. The OR values between different exposures and outcomes were calculated, and P = .05 was set as the significance threshold. It was considered that when P < .05, the study results were statistically significant.
We first employed a 2-sample MR framework to assess the causal relationships between genetically predicted levels of serum bilirubin subtypes (direct, indirect, and total bilirubin), cholelithiasis, and the risks of 3 hepatobiliary-pancreatic cancers: GDC, LIHC, and pancreatic adenocarcinoma. Follow-up sensitivity analysis was performed for each analysis, and the sensitivity analysis results of different analyses were presented in the supplementary data (Figs. S1–S15, Tables S1–S6, Supplemental Digital Content, https://links.lww.com/MD/Q508, https://links.lww.com/MD/Q509, https://links.lww.com/MD/Q510). All analyses in this study passed a pleiotropic test, and sensitivity analyses supported the robustness of the results.
The overall results are summarized in Figure 2. A statistically significant causal effect was observed for DBIL on the risk of GDC (OR = 3.07, 95% confidence interval [CI]: 1.02–9.22, β = 1.121, P = .046). This indicates that genetically predicted higher levels of DBIL are associated with an increased risk of developing GDC. In this analysis, 59 SNPs were included as IV, all of which surpassed the threshold for strength (F-statistic > 10), minimizing the risk of weak instrument bias (Fig. 2).
In contrast, none of the other exposure-outcome pairs demonstrated statistically significant causal associations. The associations between TBIL, IBIL, CHOL, and GDC, as well as all exposure factors with LIHC and pancreatic adenocarcinoma, yielded P-values > .05, suggesting a lack of robust evidence for a causal relationship based on the available genetic instruments (Fig. 2).

3.2. Validation and sensitivity analyses for the DBIL-GDC association
Given the significant finding between direct bilirubin and gallbladder cancer, we further validated this association using an independent outcome dataset (GCST001404) and conducted a comprehensive set of sensitivity analyses to assess the robustness of the result.
The causal estimate was consistent across both outcome datasets. In the primary IEU-A-1057 dataset, the IVW estimate was OR = 3.07 (95% CI: 1.02–9.22, P = .046). The validation analysis in the GCST001404 dataset yielded a highly similar result (IVW OR = 3.04, 95% CI: 1.01–9.12, P = .048; Fig. 3). The scatter plots of SNP effects showed a consistent positive slope for the IVW method across both datasets, visually supporting the directional effect of DBIL on increasing GDC risk (Fig. 4B, D).
Results from supplementary MR methods largely supported the primary IVW finding. The weighted median method, which is robust to invalid instruments if <50% of the weight comes from pleiotropic variants, also indicated a significant positive effect (Fig. 3). The MR-Egger regression intercept test revealed no evidence of directional pleiotropy (Egger intercept: −0.019, P = .63 for IEU-A-1057; −0.020, P = .59 for GCST001404; Table 1), suggesting that the causal estimate is unlikely to be biased by unbalanced horizontal pleiotropy. Cochran’s Q statistic indicated no significant heterogeneity among the instrumental variants (Q = 55.63, P > .05 for IEU-A-1057; Q = 55.86, P > .05 for GCST001404; Table 1).
Leave-one-out sensitivity analysis further confirmed the robustness of the association. Sequentially removing each SNP from the analysis did not substantially alter the overall causal estimate, indicating that the finding was not driven by any single influential genetic variant (Fig. 5). The funnel plots for both datasets showed a symmetrical distribution of SNP estimates, providing additional visual evidence against significant pleiotropy and supporting the reliability of the results (Fig. 4A, C).

3.3. Assessment of reverse causality
To examine the possibility of reverse causation – that genetic predisposition to GDC might influence circulating DBIL levels – we performed a bidirectional MR analysis. Using 15 genetic instruments associated with GDC (IEU-A-1057 dataset), we found no evidence for a causal effect of GDC on DBIL (IVW OR ≈ 1, P = .9869; Fig. 6A). The scatter plot (Fig. 6B), funnel plot (Fig. 6C), and forest plot (Fig. 6D) of the reverse MR analysis collectively support the conclusion. This absence of a effect in the reverse direction strengthens the inference that the causal direction is from higher DBIL levels to increased GDC risk, rather than the converse; It should also be noted that insufficient eligible IV were identified during the reverse MR analysis of the GDC GWAS dataset (GCST001400), indicating that the reverse causal relationship in this MR analysis can be safely rejected.

4. Discussion
In the occurrence and development of malignant tumors in the hepatobiliary and pancreatic system, abnormalities in bilirubin metabolism and cholelithiasis are believed to play an important role. Bilirubin, as one of the main components of bile, is not only a key participant in the balance of soluble cholesterol but may also be associated with tumorigenesis through its antioxidant mechanisms. Studies have shown that elevated serum bilirubin levels are closely associated with the diagnosis of PADC and other hepatobiliary malignant tumors, and in some patients, they may even serve as one of the auxiliary diagnostic markers[19]; particularly in the context of obstructive jaundice, the different distributions of bilirubin levels (total bilirubin, direct, and indirect bilirubin) have been shown to have statistical significance in predicting the presence of malignant lesions.[20,21]
Gallstone disease is considered one of the primary risk factors for gallbladder cancer. Data indicates that over 80% of gallbladder cancer patients also have gallstones, particularly larger, long-standing gallstones, which are believed to carry a higher carcinogenic risk.[22] Recent studies have further revealed that the formation of pigment gallstones is closely associated with disorders in bilirubin metabolism. These stones are typically formed by the deposition of bilirubin calcium salts and are commonly found in cases of chronic cholangitis or hemolytic states. Prolonged stimulation of the bile duct epithelium may trigger the carcinogenic process.[23]
In this study, we systematically investigated the causal relationship between serum bilirubin levels, cholelithiasis, and 3 hepatobiliary pancreatic malignancies (GDC, LIHC, and PADC) using 2-sample and 2-way MR methods. The results showed a significant causal relationship between serum DBIL and GDC, while no significant association was found between other exposure factors and outcomes. These findings provide new clues to understanding the underlying mechanisms of bilirubin in developing hepatobiliary and pancreatic malignancies and provide a scientific basis for future prevention and treatment strategies for related diseases.

4.1. Advantages of MR and its application in this study
MR, as a widely used causal inference method in epidemiological studies, has many advantages. Firstly, the MR method uses genetic variation as an IV to effectively control confounding factors through random allocation of genes, avoiding bias in traditional observational studies.[6] Secondly, MR can eliminate the interference of reverse causality to more accurately reveal the causal link between exposure factors and outcomes. In this study, through rigorous MR design, we selected significance level (P < 5e–08) as the threshold to screen genetic variation as an IV. We conducted various sensitivity analyses to verify the robustness of the results.[17,18]

The causal relationship between direct bilirubin and gallbladder cancer.

One of the main findings of this study is the significant causal relationship between serum DBIL and GDC. The OR value calculated by inverse variance weighting (IVW) was 3.07, the 95% CI was 1.02 to 9.22, and the P-value was .04, indicating that elevated direct bilirubin level was a risk factor for gallbladder cancer. Further, reverse MR Analysis found no reverse causation, which supports a causal link between direct bilirubin and gallbladder cancer. Bilirubin, as an intermediate product of heme metabolism, has traditionally been considered a harmful waste product. However, recent studies have shown that bilirubin has antioxidant, anti-inflammatory, and cell-protective effects.[1,2] Our MR suggests that the link “direct bilirubin → gallbladder cancer” is biologically plausible: in the presence of biliary microbiota, direct bilirubin is readily deconjugated by bacterial β-glucuronidase, yielding unconjugated bilirubin that precipitates with Ca²⁺ as calcium bilirubinate; together with increased epithelial mucins (e.g., MUC5AC) and biofilm-driven nucleation, this promotes pigment stone formation and bile stasis.[24] Long-standing mechanical irritation/infection by stones induces chronic cholecystitis and progression along the metaplasia–dysplasia–carcinoma sequence.[25] In parallel, cholestasis enriches hydrophobic secondary bile acids (DCA/LCA), causing oxidative and mitochondrial injury and activating pro-inflammatory pathways (COX-2, NF-κB, IL-6/STAT3), thereby maintaining a tumor-promoting microenvironment.[26] Epidemiologically, chronic Salmonella Typhi carriage is strongly associated with gallbladder cancer, supporting synergy between infection and the stone - inflammation pathway.[27] Genetically, variants in UGT1A1, ABCC2, and SLCO1B1/1B3 determine the production and biliary excretion of conjugated bilirubin, reinforcing the causal chain from upstream exposure to intermediate phenotypes (stones/cholestasis) and downstream risk.[28]
The causal relationship between other exposure factors and outcomes.

In addition to the significant association between direct bilirubin and gallbladder cancer, no significant causal relationship was found between other exposure factors in this study (such as TBIL, IBIL, and CHOL) and hepatobiliary and pancreatic malignancies. For TBIL and IBIL, although their roles in other diseases (such as cardiovascular diseases) have been reported no significant causal link was shown between them and hepatobiliary and pancreatic malignancies in this study.[3,4] Possible reasons include the fact that the mechanisms of action of these exposures are more complex or that their effects on the hepatobiliary pancreatic system are insufficient to be detected by current analytical methods. However, previous research indicates that bilirubin exhibits organ-specific bidirectional effects in cancer development: A MR study demonstrates that in hepatocellular carcinoma and cholangiocarcinoma (CCA), bilirubin metabolites (e.g., “Bilirubin [E, Z or Z, E] levels”) serve as independent risk factors (OR > 1), whose elevation significantly increases tumor risk by inducing hepatocyte oxidative stress and disrupting bile acid homeostasis; whereas in other cancers, UGT1A1-driven bilirubin elevation shows protective effects – reducing risks of squamous cell lung cancer (OR = 0.80, 95% CI 0.65–0.99) and Hodgkin lymphoma.[29] This paradoxical effect stems from organ-specific metabolic microenvironments: In the liver, bilirubin accumulation directly promotes DNA oxidative damage and interferes with bile acid metabolic pathways, while in pulmonary and lymphatic systems, it suppresses inflammatory microenvironments by neutralizing free radicals (e.g., viral oxidative stress in EBV-associated lymphoma).[30] Concurrently, the diet-microbiome axis plays a key regulatory role – high-fat diets upregulate mollicutes RF9 microbiota (OR = 1.71) to exacerbate hepatic injury, whereas zinc intake enriches anaerotruncus, mediating 31% of the protective effect against biliary tract cancer.[31] CHOL, a common digestive disease, has previously been associated with GDC, LIHC, and PADC in observational studies.[5] The MR Analysis in this study did not find a significant causal relationship between CHOL and these malignancies. This may be because the etiology and pathological mechanism of CHOL are complex, involving a variety of metabolic abnormalities and inflammatory responses, and it is difficult to reveal the causal relationship between CHOL and tumorgenesis through a single MR Analysis. In addition, the formation of CHOL may be related to specific environmental factors or lifestyle factors that are not adequately reflected in the genetic data in this study.

4.2. Limitations and future directions of the study
Despite the methodological rigor of this study, there are still some limitations. First, although the MR method can effectively control confounding factors, its results still depend on the validity of the selected IV.[6] If the IV is associated with other unaccounted confounding factors, this may result in bias. Secondly, although we chose a variety of MR methods for analysis, the statistical assumptions and sensitivities of different methods are different, which may affect the interpretation of the results.[12,13] Finally, although this study found a causal association between DBIL and GDC, the limited sample size, especially in the outcome variables (such as GDC, LIHC, and PADC), may affect the statistical power and generalization of the results. Future studies should consider expanding the sample size, especially for rare outcome variables like gallbladder cancer.
Despite the robust causal inference offered by the MR design, our study has several limitations that warrant consideration. First, and most importantly, the sample sizes for the cancer outcomes, particularly for GDC (41 cases), were relatively small. Although we employed rigorous genetic instrument selection to maximize statistical power, the limited number of cases inevitably reduces the precision of our effect estimates and increases the risk of false negatives (type II errors). This means that while we identified a significant causal effect of DBIL on GDC, we might have lacked sufficient power to detect genuine associations for other exposure-outcome pairs, such as between CHOL and GDC. The findings for GDC, although significant, should be interpreted with caution and require validation in larger, well-powered GWAS cohorts in the future. Similarly, the nonsignificant results for hepatocellular carcinoma and pancreatic cancer do not definitively rule out a causal relationship but may reflect power limitations.
The GWAS summary statistics for the exposure factors (bilirubin and CHOL) were primarily derived from European ancestry populations, whereas the outcome data for gallbladder cancer (the IEU-A-1057 and GCST001404 datasets) were based on mixed or other ancestry populations. This cross-ancestry data usage introduces a potential risk of bias due to population stratification, as the genetic architecture and linkage disequilibrium patterns of the IV may differ across ethnic groups. Although we employed rigorous quality control and standard MR assumptions to minimize bias, and the consistency of our main findings across 2 independent outcome datasets provides some reassurance, the residual confounding from population structure cannot be entirely ruled out. Therefore, the generalizability of our causal estimates to non-European populations may be limited. Future studies utilizing large-scale, ancestry-matched GWAS data are warranted to validate and extend our findings.
Finally, given the limitations of MR Analysis, future studies should consider combining other causal inference methods, such as Bayesian network analysis and multivariate MR, to improve the reliability and explanatory power of research results.[12] In addition, the specific biological mechanism of bilirubin in hepatobiliary and pancreatic malignancy should be further explored, and the effects of bilirubin on the cellular level and tissue microenvironment should be discussed. Environmental and lifestyle factors: Given the correlation between cholelithiasis and gallbladder cancer, future studies should comprehensively analyze nongenetic factors such as environmental factors and lifestyle factors to explore how these factors interact with genetic factors to affect the risk of hepatobiliary and pancreatic malignancies.[5] Integration of multi-omics data: With the development of omics technology, combining multi-omics data (such as metabolomics, proteomics, etc) with genetic data may provide a new perspective for revealing the complex relationship between bilirubin and hepatobiliary pancreatic system cancer.

5. Conclusion
This study confirms that direct bilirubin is a potential risk factor for gallbladder cancer, with its mechanism potentially related to bile acid metabolism dysregulation and oxidative stress; however, no significant causal relationship was found between several exposure factors (TBIL, IBIL, and CHOL) and hepatobiliary-pancreatic malignancies, suggesting complex associations between these exposures and cancer risk. Future research should validate biomarker stability through multicenter prospective cohorts and employ multi-omics technologies to elucidate biological functional differences among bilirubin subtypes; for clinical translation, researchers could develop dynamic prediction models integrating direct bilirubin levels with genetic risk scores to enable early screening of high-risk populations, and explore preventive interventional strategies using targeted bile acid pathway modulators.

Acknowledgments
This study was successfully completed thanks to all data providers and study participants. Special thanks to OpenGWAS, FinnGen, UK Biobank, GWAS catalog, and eLIFE databases for providing valuable data resources.

Author contributions
Conceptualization: Tao Luo.
Data curation: Tao Luo, Shangru Yang, Jianqi Qin.
Formal analysis: Xuean Zhao.
Funding acquisition: Wence Zhou.
Investigation: Xuean Zhao, Ade Su, Jianqi Qin.
Methodology: Tao Luo, Senxin Wang, Shangru Yang, Ade Su, Jianqi Qin.
Project administration: Yating Liu, Wence Zhou.
Resources: Tao Luo, Ade Su.
Software: Tao Luo, Shangru Yang.
Supervision: Wence Zhou.
Validation: Ade Su.
Visualization: Tao Luo.
Writing – original draft: Tao Luo.
Writing – review & editing: Xuean Zhao, Jianqi Qin.

Supplementary Material