We conducted a two-sample bidirectional MR analysis to explore the genetic underpinnings of the causal relationship between PBC and IBD. Initially, a forward MR analysis was executed, wherein PBC served as the exposure variable, while UC and CD were considered as outcome variables. Subsequently, a reverse MR analysis was undertaken, with UC and CD as exposure variables and PBC as the outcome variable. The validation analysis conducted a bidirectional MR analysis, treating PBC and IBD as the exposure and outcome variables respectively. Our genetic causal assessment relied upon three fundamental assumptions of MR analysis: 1) the IVs are strongly correlated with the exposure variables, 2) the IVs are not associated with the outcome variables or confounding factors, and 3) the IVs solely influence outcomes through the exposure [23]. The GWAS summary data for both exposures and outcomes utilized in this investigation were sourced from publicly accessible databases, thus obviating the necessity for ethical statements and informed consent. Detailed information regarding the data employed in this study can be found in Supplementary Table S2.

The GWAS summary data for PBC were retrieved from the IEU OpenGWAS database. This dataset comprised 11,375 individuals of European descent, consisting of 2,861 PBC cases and 8,514 controls, and encompassed 119,756 SNPs. Participants self-reported British or Irish ancestry. Diagnosis of PBC adhered to established criteria, including serological evidence of AMAs at a titer of 1:40 or higher, liver biopsy findings consistent with PBC pathology, and biochemical markers indicative of the condition [25]. Genotyping of participants utilized the Illumina iSelect HD custom genotyping array, followed by the extraction of normalized probe intensities adhering to predefined laboratory quality control metrics, and genotype assignments conducted using optiCall software [26]. Further detailed information regarding the dataset can be found in the associated publication [25]. Conversely, GWAS summary data for UC and CD were sourced from the Finnish consortium. The UC dataset included 214,620 European participants, with 4,320 cases and 210,300 controls, encompassing 16,380,459 SNPs. The CD dataset comprised 211,107 European participants, with 807 cases and 210,300 controls, and included 16,380,453 SNPs. The IBD GWAS dataset was derived from a European-ancestry cohort of 218,792 individuals, comprising 5,673 cases and 213,119 controls, and provided genotypic information for 16,380,466 SNPs. Case identification was facilitated through the application of the M13 code within the International Classification of Diseases-Tenth Revision (ICD-10) classification system. Genotyping utilized Illumina and Affymetrix chip arrays (Illumina Inc., San Diego, CA, USA, and Thermo Fisher Scientific, Santa Clara, CA, USA, respectively). Researchers interested in further elucidation regarding the datasets employed are encouraged to consult the FinnGen website for Supplementary Material.

Ensuring the robustness and credibility of MR analysis outcomes hinges upon rigorous selection criteria for IVs. This study implemented a comprehensive IV screening protocol. Initially, SNPs exhibiting significant correlations with the exposure of interest were identified from the summary data of GWAS, employing a stringent correlation threshold of P < 5 × 10^-8. Subsequently, a criterion of F statistic >10 was applied to ascertain strong correlation, calculated using the formula: F = R2(N-K-1)/K(1-R2). Should the prerequisite number of SNPs for MR analysis not be attainable under these conditions, adjustments to the correlation threshold were informed by pertinent literature [27–29]. In the subsequent phase, potential biases stemming from linkage disequilibrium (LD) among SNPs were addressed by eliminating those exhibiting strong LD, defined by an r2 < 0.001 and a clumping distance of 10,000 kb. Following this, SNPs meeting the criteria were cross-referenced with GWAS summary data related to the outcome variable, supplemented by proxy SNP identification via online LD link platform when necessary. Palindromic SNPs, susceptible to confounding due to intermediate allele frequencies, were then systematically removed. Furthermore, SNPs exhibiting correlations with the outcome variable were excised using thresholds consistent with those employed for the exposure variable. Finally, potential confounding factors were meticulously considered and adjusted for in both forward and reverse MR analyses. In the forward MR analyses investigating the causal relationship between PBC as exposure and UC as outcome, confounding variables among which were familial history of IBD, antecedent gastrointestinal infections, and the utilization of oral contraceptives [30]. And the forward MR analyses of CD as outcome, confounding variables includes smoking, family history of IBD, antecedent gastrointestinal infections, and the utilization of oral contraceptives [31]. In parallel, in the reverse MR analyses exploring UC or CD as the exposure and PBC as the outcome, confounding variables encompassed variables such as smoking, history of infections, and hormone replacement therapy [17, 32, 33].

We employed a diverse array of MR analysis to comprehensively evaluate the genetic causality between exposure and outcome variables. Our primary analytical approach was the random-effects inverse variance weighted (IVW) method, widely recognized as the most statistically robust technique in MR analysis. Additionally, we employed four supplementary MR analysis methods: MR Egger, weighted median, simple mode, and weighted mode. Furthermore, three MR analysis methods—maximum likelihood, penalized weighted median, and fixed effects IVW—were utilized for validation results. Consistent with established MR analysis literature, random-effects IVW represents the gold standard in MR analysis. In instances where the findings of alternative MR methods diverged from those of random-effects IVW, precedence was given to the conclusions drawn from random-effects IVW analysis [20–22]. The IVW method offers the most precise estimates under the assumption of instrument validity [23]. However, it is important to note that in the presence of horizontal pleiotropy, the IVW method may produce biased estimates of causal effects [34].

To ensure the robustness of our MR analysis findings, we conducted a comprehensive set of sensitivity analyses. Initially, we executed tests to evaluate heterogeneity and horizontal pleiotropy. Heterogeneity was assessed using Cochran’s Q statistic for MR-IVW and Rucker’s Q statistic for MR Egger. For horizontal pleiotropy assessment, we employed the intercept test of MR Egger and the global test of MR Pleiotropy Residual Sum and Outlier (MR-PRESSO). Subsequently, outlier detection was conducted, employing the radial variants of IVW and MR Egger, and integrate a more robust distortion test from MR-PRESSO. It is imperative to highlight that identified outliers necessitate exclusion followed by reassessment to ascertain genetic causality. Additionally, we performed Leave-one-out analyses to identify individual SNPs influencing MR analysis outcomes. Any single SNPs found to exert influence on MR results underwent exclusion and subsequent reevaluation for genetic causality. Lastly, we subjected the MR analysis to a normal distribution test using MR Robust Adjusted Profile Score (MR-RAPS).

The statistical analyses in this study were conducted using R version 4.1.2. Genetic causal assessment was performed utilizing the “TwoSampleMR” software package within the R environment. A significance threshold of P < 0.05 was applied to determine genetic causality. Specifically, positive genetic causality was inferred when the odds ratio (OR) > 1 at P < 0.05, whereas negative genetic causality was indicated by an OR < 1 at P < 0.05. Heterogeneity and horizontal pleiotropy were assessed using significance thresholds, whereby P > 0.05 suggested the absence of these phenomena. Additionally, the normal distribution assumption was evaluated, with a P > 0.05 indicating compliance with this assumption.

After accounting for LD, we identified 22 SNPs exhibiting a strong correlation with PBC. These SNPs are all present in the GWAS summary data for UC. Notably, there were no palindromic SNPs among them, and none of the 22 SNPs demonstrated an association with UC. It is also important to highlight that there was no confounding SNPs detected. Consequently, these 22 SNPs were utilized as IVs for assessing the genetic causality between PBC on UC (Supplementary Table S3).

Random-effects IVW results indicated a positive genetic causality between PBC on UC (P = 0.024, OR = 1.065, 95% confidence interval [CI] = 1.008–1.124). Additional supplementary methods, except MR Egger, yielded consistent findings with the random-effects IVW results (Supplementary Figures S1, S2A). Consistency was also observed across three verification methods, reinforcing the random-effects IVW analysis (Supplementary Figure S1). Cochran’s Q statistic for MR-IVW and Rucker’s Q statistic for MR Egger analysis demonstrated heterogeneity in the genetic causal assessment of PBC on UC (P < 0.05). While MR Egger’s intercept test did not reveal horizontal pleiotropy (P > 0.05), the global test of MR-PRESSO did (P < 0.05) (Supplementary Table S7). Examination of radial variants in the IVW and MR Egger plots identified outliers in the MR analysis (Supplementary Figure S2B). The distortion test of MR-PRESSO identified two significant outliers (rs522127, rs8067378) (Supplementary Table S7). Leave-one-out analysis indicated that genetic causal assessment was influenced by a single SNP (rs7775055) (Supplementary Figure S2C). MR-RAPS analysis supported a normal distribution (P > 0.05) (Supplementary Figure S2D; Supplementary Table S7). A second round of MR analysis was conducted after removing the two identified outliers. Random-effects IVW results confirmed a positive genetic causal relationship between PBC on UC (P = 0.004, OR = 1.069, 95% CI = 1.022–1.119). Among the four supplementary methods, MR Egger and the simple mode produced results contrary to the random-effects IVW findings, whereas the weighted median and weighted mode analyses supported the random-effects IVW results (Supplementary Figures S1, S3A). Consistency was maintained across three verification methods (Supplementary Figure S1). Heterogeneity persisted in the genetic causal assessment (P < 0.05). MR Egger’s intercept test again did not detect horizontal pleiotropy (P > 0.05), while the MR-PRESSO global test did (P < 0.05) (Supplementary Table S7). Radial variant analysis in IVW and MR Egger plots continued to reveal outliers (Supplementary Figure S3B). The MR-PRESSO distortion test detected one significant outlier (rs72678531) (Supplementary Table S7). Leave-one-out analysis showed that the genetic causal assessment was not affected by individual SNPs (Supplementary Figure S3C). MR-RAPS analysis continued to support a normal distribution (P > 0.05) (Supplementary Figure S3D; Supplementary Table S7).

Finally, a third round of MR analysis was conducted after the removal of one outlier. The random-effects IVW results indicated a positive genetic causality between PBC on UC (P < 0.001, OR 95% CI = 1.081 [1.037–1.127]). Among the supplementary methods, all except MR Egger corroborated the random-effects IVW findings (Figures 1, 2A). Additionally, the three verification methods were consistent with the random-effects IVW results (Figure 1). There was no evidence of heterogeneity or horizontal pleiotropy in the genetic causal assessment (P > 0.05) (Table 1). Although the radial variants of IVW and MR Egger identified outliers (Figure 2B), the MR-PRESSO distortion test detected no outliers (Table 1). The analysis was not influenced by single SNPs (Figure 2C), and showed a normal distribution (P > 0.05) (Figure 2D; Table 1).

After adjusting for LD, we identified 22 SNPs exhibiting strong correlation with PBC. All 22 SNPs are present in the GWAS summary data for CD. Notably, there were no palindromic SNPs among them. None of these 22 SNPs demonstrated an association with CD. It is crucial to mention that the confounding SNP rs11065979 was excluded from the analysis. Consequently, 21 SNPs were selected as IVs for the genetic causality assessment between PBC on CD (Supplementary Table S4).

Random-effects IVW analysis indicated a positive genetic causality between PBC on CD (P = 0.030, OR 95% CI = 1.114 [1.011–1.227]). Supplementary analyses using weighted median and weighted mode methods corroborated the random-effects IVW findings, whereas MR Egger and simple mode yielded inconsistent results (Supplementary Figures S1, S4A). Nonetheless, all three validation methods were in agree with the random-effects IVW results (Supplementary Figure S1). Cochran’s Q statistic for MR-IVW and Rucker’s Q statistic for MR Egger revealed significant heterogeneity (P < 0.05). Although MR Egger’s intercept test did not detect horizontal pleiotropy (P > 0.05), the global test of MR-PRESSO confirmed its presence (P < 0.05) (Supplementary Table S7). Outlier detection analyses identified outliers in the radial variants of IVW and MR Egger (Supplementary Figure S4B), and the MR-PRESSO distortion test also detected a significant outlier (rs10931468) (Supplementary Table S7). Leave-one-out analysis indicated that the genetic causal assessment was influenced by individual SNPs (rs79513546, rs35188261, rs7775055, rs2267407) (Supplementary Figure S4C). MR-RAPS analysis demonstrated that the genetic causal assessment was normally distributed (P > 0.05) (Supplementary Figure S4D; Supplementary Table S7).

Due to the inability of the sensitivity analysis to support the initial MR analysis results, a second round of MR analysis was conducted after excluding the outlier. Consistent with the initial findings, the random-effects IVW analysis continued to indicate a positive genetic causal relationship between PBC on CD (P = 0.002, OR 95% CI = 1.136 [1.047–1.232]). Despite the conflicting results from MR Egger and simple mode, the weighted median and weighted mode analyses supported the random-effects IVW findings (Figures 1, 3A). Furthermore, the three verification methods consistently indicated a positive genetic causal relationship between PBC and CD (Figure 1). No evidence of heterogeneity or horizontal pleiotropy was observed (P > 0.05) (Table 1). Although radial variants in IVW and MR Egger analyses identified outliers (Figure 3B), the MR-PRESSO distortion test did not detect any outliers (Table 1). The genetic causal assessment was not influenced by individual SNPs (Figure 3C) and adhered to a normal distribution (P > 0.05) (Figure 3D; Table 1).

At a correlation threshold of P < 5 × 10^-8, an insufficient number of SNPs were available for subsequent MR analysis. Consequently, we adjusted the correlation threshold to P < 5 × 10^-6. This adjustment yielded 30 SNPs that are strongly associated with UC, 13 of which are present in the GWAS summary data for PBC. No palindromic SNPs were identified. One SNP associated with PBC, rs12946510, was excluded from the analysis. No confounding SNPs were detected. Ultimately, we identified 12 SNPs to serve as IVs for the genetic causality assessment of UC on PBC (Supplementary Table S5).

The random-effects IVW analysis revealed no evidence of genetic causality between UC on PBC (P = 0.575, OR 95% CI = 0.942 [0.763–1.162]). This result was corroborated by four supplementary methods (Supplementary Figures S5, S6A). Notably, the penalized weighted median analysis suggested a negative genetic causality between UC on PBC (Supplementary Figure S5). However, the sensitivity analysis indicated heterogeneity in the genetic causal assessment (P < 0.05). While MR Egger’s intercept test suggested the absence of horizontal pleiotropy (P > 0.05), the global test of MR-PRESSO provided evidence for horizontal pleiotropy (P < 0.05) (Supplementary Table S7). Radial variant analysis in both IVW and MR Egger revealed outliers in the genetic causal assessment (Supplementary Figure S6B), and the distortion test of MR-PRESSO identified an outlier (rs3024493) (Supplementary Table S7). The leave-one-out analysis did not identify any single SNPs influencing the outcome of the genetic causal assessment (Supplementary Figure S6C). MR-RAPS analysis indicated a normal distribution (P > 0.05) (Supplementary Figure S6D; Supplementary Table S7). Subsequent MR analysis, after removing an outlier, still demonstrated no genetic causal relationship between UC on PBC (P = 0.174, OR 95% CI = 0.882 [0.735–1.057]). Except MR Egger, this finding was supported by three other complementary methods (Supplementary Figures S5, S7A). Among these verification methods, the penalized weighted median was consistent with MR Egger, indicated a negative genetic causal relationship between UC on PBC. The maximum likelihood and fixed effects IV analyses were consistent with the random-effects IVW results (Supplementary Figure S5). No heterogeneity was observed in the sensitivity analysis (P > 0.05). The MR Egger intercept test indicated the presence of horizontal pleiotropy (P < 0.05), whereas the global test of MR-PRESSO did not detect horizontal pleiotropy (P > 0.05) (Supplementary Table S7). Radial variant analyses in both IVW and MR Egger demonstrated the presence of outliers (Supplementary Figure S7B), but the distortion test of MR-PRESSO did not identify outliers (Supplementary Table S7). Importantly, the leave-one- out analysis showed that two individual SNPs (rs7936070, rs3197999) influenced the genetic causal assessment (Supplementary Figure S7C). MR-RAPS continued to demonstrate that the genetic causal assessment conformed to a normal distribution (P > 0.05) (Supplementary Figure S7D; Supplementary Table S7).

Finally, a third round of MR analysis was conducted after excluding two SNPs that impacted the genetic causal assessment results. The random-effects IVW analysis revealed a negative genetic causality between UC on PBC (P = 0.003, OR 95% CI = 0.788 [0.672–0.924]). Among the four supplementary methods employed, MR-Egger and weighted median analyses supported this finding, whereas simple mode and weighted mode analyses suggested no genetic causal relationship between UC on PBC (Figures 4, 5A). Nevertheless, all three validation methods provided evidence of a negative genetic causality between UC on PBC (Figure 4). Sensitivity analyses did not detect heterogeneity or horizontal pleiotropy (P > 0.05) (Table 1). Furthermore, there was no evidence of outliers (Figure 5B; Table 1). The genetic causality assessment remained robust when individual SNPs were considered (Figure 5C), and the results were normally distributed (P > 0.05) (Figure 5D; Table 1).

Under the initial correlation threshold of P < 5 × 10^-8, it was not possible to obtain sufficient SNPs for subsequent MR analysis. Therefore, the correlation threshold was adjusted to P < 5 × 10^-6, resulting in the identification of 22 SNPs strongly associated with CD, of which 4 SNPs were also present in the GWAS summary data for PBC. Notably, there were no palindromic SNPs, no SNPs related to PBC, and no confounding SNPs identified. Consequently, these 4 SNPs were utilized as IVs for the genetic causal assessment of CD on PBC (Supplementary Table S6).

The random-effects IVW analysis indicated no genetic causal relationship between CD on PBC (P = 0.431, OR 95% CI = 0.956 [0.856–1.068]). Four supplementary methods corroborated the findings of the random-effects IVW analysis (Figures 4, 6A). Similarly, three verification methods produced consistent results (Figure 4). Both Cochran’s Q statistic for MR-IVW and Rucker’s Q statistic for MR Egger indicated no heterogeneity in the genetic causal assessment of CD on PBC (P > 0.05). Moreover, the intercept test of MR Egger and the global test of MR-PRESSO did not reveal the presence of horizontal pleiotropy (P > 0.05) (Table 1). Due to the limited number of IVs, radial variants of IVW and MR Egger are not illustrated. The MR-PRESSO distortion test did not detect any outliers (Table 1). Leave-one-out analysis demonstrated that the genetic causality assessment was robust to the retention of individual SNPs (Figure 6B). Additionally, MR-RAPS analysis indicated a normal distribution (P > 0.05) (Figure 6C; Table 1).

Following LD pruning, we first identified 22 SNPs that were strongly associated with PBC. All of these variants were successfully retrieved from the IBD GWAS summary data, and none were palindromic. Subsequent checks confirmed that none of the 22 SNPs exhibited an independent association with IBD. We then proceeded to exclude one confounding SNP (rs11065979) and to remove four outlier variants (rs2267407, rs522127, rs72678531, rs8067378) in advance. After these quality-control steps, a final set of 17 SNPs was retained as IVs for estimating the causal effect of PBC on IBD (Supplementary Table S8). For the reverse direction, LD pruning yielded 9 SNPs that were strongly correlated with IBD. These variants were all available in the PBC GWAS summary statistics, with no palindromic or confounding SNPs detected, and none showed any association with PBC. Accordingly, all 9 SNPs were directly adopted as IVs for assessing the causal impact of IBD on PBC (Supplementary Table S9).

The random-effects IVW method revealed a significant positive genetic causal relationship of PBC on IBD (P < 0.001, OR 95% CI = 1.076 [1.042–1.110]). This finding was consistently supported by three supplementary methods and further corroborated by three additional verification approaches (Figures 7, 8A). Conversely, the reverse-direction analysis, using the same random-effects IVW model, showed no genetic causal effect of IBD on PBC (P = 0.357, OR 95% CI = 0.898 [0.714–1.129]). This null result was likewise confirmed by four supplementary methods and three verification procedures (Figures 7, 8C). Sensitivity analyses indicated the absence of heterogeneity, horizontal pleiotropy, and outlier variants, with the residuals conforming to a normal distribution (Table 1). Furthermore, leave-one-out analysis demonstrated that no single SNP disproportionately influenced the causal estimates (Figures 8B,D).