The foundational data for this investigation was retrieved from publicly available summary-level datasets from GWAS. Univariate Mendelian Randomization (UVMR), Multivariable Mendelian Randomization (MVMR), genetic correlation, and colocalization analyses were used to elucidate the causal interplay between autoimmune disorders and outcome phenotypes.

The selection of Instrumental Variables (IVs) for exposure was grounded in a tripartite criterion: i) the nominated genetic determinant, earmarked as the instrumental variable, must display a robust affiliation with the exposure; ii) the genetic determinant must not be intertwined with any potential confounders; and iii) the influence of genetic determinants on the outcome is channeled exclusively through its interaction with the exposure, thus eliminating the prospect of secondary routes (13). The architectural blueprint of the MR is illustrated in Figure 1 and Table 1 , along with Supplementary File 1 , which provides a comprehensive exposition of the summary statistics data repositories.

It is imperative to note that all encompassed GWAS investigation procured endorsements from the relevant academic oversight committees. Given that our study was based on a secondary analysis of publicly disclosed datasets, further ethical vetting was not required.

In order to preserve the integrity of our Mendelian Randomization approximations, the chosen Single Nucleotide Polymorphisms (SNPs) were obligated to align with the ensuing benchmarks:

For autoimmune diseases, all from large abstract-level GWAS studies, ulcerative colitis (UC) and Crohn’s disease (CD) from de Lange KM et al. (16), celiac disease (CeD) from Trynka et al. (17), multiple sclerosis (MS) from International Multiple Sclerosis Genetics Consortium (IMSGC) (18), RA from Ha E et al. (19), systemic lupus erythematosus (SLE) from Bentham J et al. (20), type 1 diabetes (T1D) from Forgetta V et al. (21), psoriasis (PsO) from Tsoi LC et al. (22), primary sclerosing cholangitis (PSC) from International PSC Study Group (IPSCSG) (23), primary biliary cirrhosis (PBC) from Cordell HJ et al. (24), ankylosing spondylitis (AS) from Cortes A et al. (25), vitiligo (ViT) from Jin Y et al. (26), and for the outcome phenotype BE from FinnGen (R9) Consortium (27).

We further obtained genetic associations for Body Mass Index (BMI) from the Genetic Investigation of Anthropometric Traits (GIANT) consortium (28), smoking from GWAS and Sequencing Consortium of Alcohol and Nicotine use (GSCAN) (29), triglycerides (TG), Low Density Lipoprotein Cholesterol (LDL-C) and High Density Lipoprotein Cholesterol (HDL-C) from Global Lipids Genetics Consortium (GLGC) (30), 25-hydroxyvitamin D (25OHD) levels from Manousaki D et al. (31).

For the UVMR study, the Wald ratio test was used for exposure with only one instrument, and the multiplicative random-effects inverse-variance-weight (IVW) method was implemented for the causative assessment of multiple IVs (comprising two or more). This approach was further enhanced by incorporating both the MR-Egger and weight median techniques. The weightage in IVW is directly related to each SNP’s Wald ratio estimate and inversely correlated with the variance estimate of each SNP’s Wald ratio (32). When all genetic markers are judged valid, IVW provides estimates that are both consistent and efficient. Conversely, the weight median method stands out when over half of the genetic markers are deemed questionable, and the MR-Egger approach is adopted when all genetic markers are refutable (33). Stringent adjustment for multiple comparisons was performed using the False Discovery Rate (FDR). Following this adjustment, a P-value < 0.05 was considered indicative of a significant causal relationship. However, instances where the raw P-value was below 0.05, but the FDR-adjusted P-value exceeded this threshold were regarded as tentative.

Given the potential confounding effects of factors, such as BMI, smoking habits, lipid profiles (LDL-C, HDL-C, and TG), and 25OHD levels on the progression from exposure to outcome, subsequent MVMR analyses were performed. This study aimed to accurately quantify the direct causative effects of exposure on the results. When juxtaposed with the UVMR paradigm, the primary supposition of MVMR focus on genetic variability associated with one or more exposures, whereas the succeeding assumptions harmonize with the UVMR framework (34). A refined investigation was undertaken to ascertain the magnitude of mediation by certain factors. The initial step was to obtain the MR effect projections for exposure in relation to the outcome phenotypes using the IVW approach. Thereafter, multivariate MR analysis was performed to ascertain the impact of nine mediating factors on the outcome while concurrently considering exposure attributes. The indirect influence of the exposure was determined by multiplying the resulting estimates for each outcome. Finally, the division of the mediation effect by the overarching effect provided insight into the relative contribution of the mediators to the overall outcome.

The LDSC regression, specifically tailored for GWAS summary data, serves as a robust approach for dissecting genetic correlations across complex diseases and traits. Notably, LDSC efficiently differentiates genuine polygenic signals from potential confounders such as cryptic relatedness and population stratification (35). A consequential genetic correlation, both statistically and quantitatively robust, signifies that an overarching phenotypic correlation is not merely attributable to environmental confounders (35). The LDSC tool, accessible at (https://github.com/bulik/ldsc), was used to scrutinize the genetic intersections between exposure and an array of outcome phenotypes.

Within the framework of UVMR analysis, several tests were conducted to validate its rigor and authenticity. The heterogeneity of the selected genetic variants was assessed using Cochran’s Q test, wherein a P-value of < 0.05 indicated pronounced discrepancies among the scrutinized SNPs (36). Employing the MR-Egger regression (37), this investigation discerned the potential for directional pleiotropy within the MR context. MR-Egger’s intercept, with a P-value < 0.05, signified the presence of consequential directional pleiotropy despite the inherent limitations of this methodology (38). The MR Pleiotropy Residual Sum and Outlier (MR-PRESSO) approach was used to identify probable outliers and delve into horizontal pleiotropy, which was inferred when the global p-value was less than 0.05 (39). By excluding such outliers, the data correction was refined. An ensuing leave-one-out analysis elucidated the impact of singular SNPs on collective outcomes (40).

R² was calculated using the formula 2×MAF×(1-MAF) ×beta², where MAF denotes the minor allele frequency for each designated SNP. The cumulative values provided a coefficient essential for power computation (41). The determination of statistical potency was anchored on the mRnd platform (42) and is accessible at https://shiny.cnsgenomics.com/mRnd/.

The SNPs of each autoimmune disease were screened according to the genetic instrument selection process described above. Power calculations for bidirectional univariable MR analyses between autoimmune diseases including CD, CeD, MS, RA, SLE, UC, T1D, PsO, PSC, PBC, AS, ViT and BE, were performed. The study reported F- statistics exceeding 60 for all instrumental variants, signifying a robust reduction in bias from weak instruments. The SNPs selected as IVs ranged from 15 to 83, accounting for an explained variance of 2.59% to 1535.64% ( Supplementary Table 1 ).

LDSC genetic correlation analyses were conducted to estimate the genetic correlation between different autoimmune diseases and BE. LDSC analysis revealed significant genetic correlations between BE and CD (r_g = 0.220, P = 0.037), RA (r_g = 0.210, P = 0.021), and UC (r_g = 0.247, P = 0.023) ( Supplementary Table 2 ). However, no genetic correlation was found with other autoimmune diseases (P > 0.05). The SNP-based liability-scale heritability (h²) ranged from 0.1% to 232.99%. Additionally, the genetic correlation between each autoimmune disease and BE was analyzed ( Figure 2 ; Supplementary Table 3 ).

A scatter plot illustrates the causal relationship between each autoimmune disease and BE ( Supplementary Figure 1 ). After adjusting for multiple comparisons, the primary IVW analysis provided strong evidence for two causal relationships ( Figure 3 ). Specifically, for each standard deviation (SD) increase in genetically predicted RA, there was a 10.3% increase in the incidence of BE (odds ration [OR] = 1.103, 95% CI 1.055-1.154, P = 1.75×10^-5, FDR = 5.25×10^-5). Furthermore, for every SD increase in CeD, the incidence of BE was reduced by 5.1% (OR = 0.949, 95% CI 0.902-0.999, P = 0.044, FDR = 0.044). We also observed suggestive evidence corresponding to a 3% increase in BE incidence with T1DM (OR = 1.033, 95% CI 1.001-1.066, P = 0.042, FDR = 0.063). Additionally, we had 96%, 100%, and 92% statistical power to detect the associations of CeD, RA, and T1D with BE, with OR values of 1.103, 0.949, and 1.033, respectively ( Supplementary Table 1 ). No other causal relationship evidence was found (P > 0.05, FDR > 0.05) ( Table 2 ). Furthermore, MVMR analysis showed that RA was an independent risk factor for BE, whereas mediator MR analysis did not identify any mediating factors ( Figure 4 ).

To avoid excessive bias effects, Cochran’s Q test was performed to analyze the sensitivity of the MR results, and no evidence of heterogeneity was observed (P>0.05). Moreover, no horizontal pleiotropy was identified using the MR-Egger intercept test (P>0.05) or the MR-PRESSO global test (P>0.05). These analyses confirmed the robustness of the findings ( Table 3 ). Leave-one-out analysis did not reveal any horizontal pleiotropy and further confirmed that the causal relationship was not influenced by any individual SNP ( Supplementary File 1 ).