Data on IBD was collected from the IEU GWAS database (https://gwas.mrcieu.ac.uk/), and GWAS summary statistics for IBD (ID: ebi-a-GCST004131), CD (ID: ebi-a-GCST004132) and UC (ID: ebi-a-GCST004133) were obtained from the European Bioinformatics Institute (EBI) database (32). Among them, IBD data includes a total of 59957 samples, including 25042 patients and 34915 healthy controls; The data of CD includes a total of 40266 samples, including 12194 patients and 28072 healthy controls; UC’s data includes a total of 45975 samples, including 12366 patients and 33609 healthy controls. The latest summary-level statistics for RA were obtained from the FinnGen GWAS results (https://r10.finngen.fi/) (33), which comprised a total of 12,555 patients with RA and 240,862 controls. All the participants were of European ancestry. We also selected the RA data (ID: ukb-b-9125) from the IEU database for validation.The flow chart of the present study is presented in Figure 1 .

LDSC (version 1.0.1) calculates the genetic correlation between two shapes to show the average degree of the sharing of genetic effects between the pairs of traits (30). Genetic correlation (rg), the most important objective of the present study, ranges from −1 to 1, with the “−” sign denoting a negative correlation and the “+” sign denoting a positive correlation. Per the LDSC manual, the first step is to convert all GWAS summary statistics into the LDSC format with the help of the default parameters of munge_sumstats.py provided in the package. The second step is to calculate genetic correlations using the -rg, -ref-ld-chr, and -w-ld-chr parameters in ldsc.py. The precalculated LD scores passed to the -ref-ld-chr and -w-ld-chr parameters can be downloaded from https://alkesgroup.broadinstitute.org/LDSCORE/. Moreover, information on European pedigrees from the 1000 Genomes Project (34), which was consistent with the European pedigrees of the GWAS samples, was used as a reference panel for chain imbalance.

In the local genetic correlation analysis, we investigated 1703 gene segments that were delineated and that should be independent of prespecified LDs. HESS, a novel computational tool, calculates local SNP heritability and measures the degree of similarity between the pairs of traits driven by genetic variation (27). However, the HESS results should be corrected using the Bonferroni method, for which a correlation was considered statistically significant at P < 0.05/1703 = 2.94E-05.

The specific procedure followed for the MR analysis is presented in Figure 2 . First, according to the correlation assumption of MR, we set the P-value to be less than 5 × 10−8 to obtain single nucleotide polymorphisms (SNPs) that are closely related to exposure. Second, LD was excluded by setting the r2 value to 0. 001 and kilobase pairs (kb) to 10 000 (35). Additionally, we manually input all exposure-related SNPs into Phenoscanner (http://www.phenoscanner.medschl.cam.ac.uk/) to check for any confounding variables related to outcomes. If confounding variables were present, the corresponding SNPs were deleted to satisfy the assumption that IVs is independent of confounding factors and outcomes.

Pleiotropy, heterogeneity, and sensitivity are important tools for the quality control of the MR analysis results. Based on the results of two polyvalence assessments, MR-Egger regression (36) and MR-PRESSO test (37), we identified abnormal SNPs, which were rejected. The results of inverse variance weighted (IVW) and MR Egger methods, including Cochran’s Q statistic and the corresponding P-value, where the P-value is the weighing of the presence and absence of heterogeneity (P < 0.05 denoting no heterogeneity), showed heterogeneity. The leave-one-out method eliminates each SNP in turn to calculate the combined effect of the remaining SNPs in order to determine whether the SNP exerts an unfavorable main effect on the overall causal relationship (38). Additionally, we analyzed the eligibility of the included SNPs by calculating the F-value (F = β2/SE2), where β is the effect value of the allele, and SE is the standard error (39). An F-value greater than 10 suggested the eligibility and absence of weak instrumental variables and vice versa (40).

TwosampleMR was used for the MR analysis, and IVW was selected as the stochastic model; this is the most critical method for evaluating analysis results. IVW evaluates the average causal effect between two traits, and the evaluation includes calculating the Wald ratio estimates for each SNP (41). The entire MR analysis is bidirectional.

The empirical Bayesian statistical framework of condFDR and conjFDR analyses was used to identify shared genetic risk loci for IBD and RA (42). First, conditional quantile–quantile (Q–Q) plots were generated, which showed overlapping SNP associations among different traits. The conditional Q–Q plot analysis sets SNPs against secondary traits as a reference and evaluates all SNPs against the primary traits. Second, condFDR was used as a reference to reorder the test statistics of the primary trait (e.g., IBD) on the basis of the strength of the association with the secondary trait (e.g., RA) to enhance the feasibility of SNP identification (29). Next, the primary and secondary traits were inverted to perform a reverse condFDR analysis. The loci of the primary and secondary traits were considered significant at condFDR ≤ 0.01. Finally, the maximum value of condFDR was set as the value of conjFDR to ensure that the FDR associated with the two traits was more accurate. The entire analysis was performed by directly eliminating SNPs around the major histocompatibility complex expansion region, which improved the accuracy of the FDR estimates (43). The detailed procedure for the conjFDR analysis is available at https://github.com/precimed/pleiofdr.

FUMA (44), an online analysis website (available at https://fuma.ctglab.nl/), allows for the positional mapping and functional annotation of new, shared, and specific locus. We considered candidate SNPs with conjFDR < 0.05 and LD r2 < 0.6 with each other as independent significant SNPs. Among them, those SNPs with LD r2 < 0.1 were designated as lead SNPs. All the candidate SNPs in LD r2 ≥ 0.6 with a lead SNP demarcated the boundaries of a genomic locus. We merged loci separated by less than 250 kb. Therefore, any candidate SNP located within the boundaries of a genomic locus was defined to belong to a single independent genomic locus. We also employed the Combined Annotation Dependent Depletion (CADD), a framework that objectively integrates diverse annotations into a unified quantitative score, which is pivotal in evaluating functional, deleterious, and pathogenic genetic variations. A CADD score exceeding 12.37 signifies harmfulness (45).

We investigated pathway enrichment for the mapped genes. Gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to identify common trends related to the function of differentially expressed genes (46), which was achieved by performing GO annotations of genes from the R package org.Hs.eg.db (version 3.1.0). The genes were mapped to the background set, and the enrichment analysis was performed using the R package clusterProfiler (version 3.14.3). The minimum gene set included 5 entities and the maximum gene set included 5000 entities, and the P-value and FDR value were set to less than 0.05 and less than 0.1, respectively. For this, the enrichment analysis tool of the Sangerbox platform (http://vip.sangerbox.com/) was used (47).

The LDSC–SEG analysis was performed by considering the tissue gene expression data as a reference and combining the GWAS data of both IBD and RA to identify tissues associated with IBD and RA (30). In this analysis, the t-statistic of each gene expressed in different tissues was calculated, the t-statistic scores of all genes were ranked in ascending order, and the top 10% of the genes were considered as the gene set corresponding to the lesion tissue. Further, the gene sets were divided into strong and weak specific expression gene sets per their intensity. In each gene set, a 100-kb window was set between the front and back of the transcribed region of each gene in order to establish the annotation information of the tissue genome. Finally, GWAS summary statistics were performed to assess the role of focused genome annotation in trait heritability. GTEx (v8) offers insights into 53 distinct tissue types, encompassing data on SNP mutations linked to quantitative traits of gene expression across various tissues (48). S-LDSC provides preprocessed reference comparison panels for 53 tissues of the GTEx project in this version, which allows researchers to more conveniently evaluate the correlation between tissues and traits. We downloaded and used pre-calculated GTEx project gene expression data as genomic annotations for each tissue to reveal the potential tissue origins of the comorbidity inheritance between IBD and RA. The detailed procedure of the LDSC–SEG analysis is available at https://github.com/bulik/ldsc/wiki/Cell-type-specific-analyse.

Utilizing Python 3.11.5, we conducted a GWAS multi-trait analysis (MTAG) between IBD and RA to ascertain the risk SNPs significantly associated with both traits (31). MTAG facilitates the joint analysis of GWAS statistical summaries for different traits without potential sample overlap between studies. It is predicated on the establishment of a variance-covariance matrix with shared effect sizes across traits (31). The GWAS derived from MTAG analysis between the two traits was submitted to the FUMA (44) online analysis platform for identification of the associated genetic risk loci.

The LDSC analysis results showed 2380 SNPs that were significantly correlated to IBD and 773 SNPs that were significantly correlated to RA.2076 and 1042 on the CD and UC side, respectively. Excluding the intercept, the probability of inheritance for IBD was 31.3%, whereas that for RA was 2.73%, and a positive correlation was observed between the two (rg = 0.1279, P = 0.0191).The correlation of RA with CD (rg = 0.1242, P= 0.0347) was also positive, but there was no significant association in terms of UC (rg = 0.1006, P= 0.0756). The results of the validation were consistent with the above performance ( Table 1 ).

The local genetic correlation mapping analysis indicated a local genetic overlap between IBD and RA, and the specific site was on two chromosomes, 6 and 12. CD aspects are expressed on chromosomes 1 and 12 and UC on chromosome 6 ( Figures 3A-C ). Localized correlations in ukb-b-9125 are enriched on chromosomes 1 and 6 ( Supplementary Figures 1A-C ).

In the MR analysis of IBD on RA orientation, 81 SNPs were studied. The IVW results showed a positive causal effect of IBD on RA risk (P = 0.0039). The leave-one-out analyses showed no potentially affecting SNP-driven causal associations, indicating the reliability of the present results ( Figure 4 ). Additionally, MR-PRESSO showed no horizontal multidirectionality (P = 0.3873). The reverse MR analysis showed no significant causal effect of RA on IBD.There was no causal relationship between IBD subtypes and RA ( Supplementary Table 1 ). There were too few reverse dependent instrumental variables, and only positive MR was performed for IBD and ukb-b-9125, with results suggesting no causality.( Supplementary Table 1 ). The F-values for all instrumental variables were greater than 10, which suggested the absence of a bias in the weak instrumental variables ( Supplementary Table 1 ).

The conditional Q–Q plots for IBD and RA showed that when the association significance of one disease increased, the other disease significantly shifted to the left from the expected zero line, which suggested the presence of genetic enrichment and a shared genetic background between IBD (including CD and UC) and RA ( Figures 5A-F ). The same conclusion was reached with respect to ukb-b-9125 ( Supplementary Figures 2A-F ). The ConjFDR analysis showed that all genes between the two traits overlapped in a way that the overlapping improved the qualitative confidence and guaranteed that the identified loci belonged to both diseases. At conjFDR < 0.01, there were 20 genetic risk loci shared by IBD GWAS and RA GWAS ( Figure 6A ; Supplementary Table 2 ), of which 17 showed the same direction of effect ( Supplementary Table 2 ). For CD and UC, the genetic risk loci were 21 and 13, respectively, of which 16 and 13 showed the same direction of effect ( Figures 6B, C ; Supplementary Tables 3 , 4 ). Ukb-b-9125 had fewer genetically risk loci between ukb-b-9125 and IBD, CD, and UC, corresponding to 1, 3, and 0, respectively ( Supplementary Figures 3A-C ).

At concFDR < 0.01, a total of 92 mapping genes ( Supplementary Table 5 ) and 1575 candidate SNPs were obtained for IBD and RA, which mainly corresponded to ncRNA_intronic (n = 544, 34.8%), intronic (n = 530, 34.0%), and intergenic (n = 327. 20.9%) functions ( Supplementary Table 6 ). Among the 20 genetic risk loci, the functional attributes intronic、intergenic、ncRNA_intronic、downstream, and ncRNA_exoni corresponded to 9, 5, 3, 2, and 1 SNPs, respectively ( Supplementary Table 2 ). Among these 20 enetic risk loci, the combined annotation-dependent depletion value of the SNP rs353341 exceeded the threshold of 12.37, which indicated that it was deleterious (45). There are 2061 and 1311 candidate SNPs for CD and UC, respectively ( Supplementary Tables 7 , 8 ). For mapping genes, there are 103 for CD and 67 for UC ( Supplementary Tables 9 , 10 ).

The GO annotation and KEGG pathway enrichment analyses were performed on the 92 mapped genes in IBD, and the results showed that in terms of biological processes, the top three aspects of gene enrichment were signaling regulation, protein metabolism regulation, and cellular response to cytokine stimulus ( Figure 7A ). In terms of cellular components, the cytosol, membrane raft, and membrane microdomain were the main directions of gene enrichment ( Figure 7B ). In terms of molecular functions, genes were mainly enriched for protein kinase binding and kinase binding ( Figure 7C ). Conversely, the notable pathways analyzed by KEGG were Epstein–Barr virus (EBV) infection, Janus kinase signal transducer and activator of transcription (JAK–STAT) signaling pathway, and Th1 and Th2 cell differentiation ( Figure 7D ). The results of the enrichment analysis for CD and UC can be seen in Supplementary Figures 4 , 5 .

Based on the tissue expression data from GTEx, we performed additional LDSC–SEG analyses to identify disease-associated tissues. At P < 0.05, IBD was associated with the whole blood, lung, spleen, small intestine terminal ileum, EBV-transformed lymphocytes, adipose visceral (omentum), nerve tibial, transverse colon, sun-exposed skin (lower leg), and uterus ( Figure 8A ; Supplementary Table 11 ). Tissues that showed a significant correlation with RA were the spleen, whole blood, EBV-transformed lymphocytes, and small intestine terminal ileum ( Figure 8B ; Supplementary Table 12 ). Crohn’s disease is closely associated with whole blood, spleen, lungs, terminal ileum of the small intestine, EBV-transforming lymphocytes, fatty viscera (omentum), fallopian tubes, and uterus ( Figures 8C ; Supplementary Table 13 ). Tissues significantly associated with UC include the lung, terminal ileum of the small intestine, whole blood, spleen, transverse colon, tibial nerve and bladder ( Figures 8D ; Supplementary Table 14 ). LDSC-SEG results for ukb-b-9125 indicated significant associations with spleen tissue as well as the terminal ileum of the small intestine; also showing significant association with bladder and mammary gland tissue ( Figures 8E ; Supplementary Table 15 ). The spleen and small intestine terminal ileum were common to both diseases, which suggested that IBD and RA might have a common tissue origin. The relevant specific analyses and their results are presented in Supplementary Tables 11 - 15 .

After performing MTAG analysis and Fuma annotation for both IBD and RA, a total of 94 genetic risk loci were identified ( Supplementary Table 16 ). The conjfdr and MTAG analysis revealed the intersection of five genes (FOXP1, RP11–95M15.1, RP11–89M16.1, PTPN2, and UBE2L3) ( Supplementary Figure 6A ). For CD, we found 80 genetic risk loci ( Supplementary Table 17 ), with the corresponding intersection genes being AC020743.4, RP11–89M16.1, PTPN2, and UBE2L3 ( Supplementary Figure 6B ). Upon analyzing UC, only 49 genetic risk loci were identified ( Supplementary Table 18 ), with the intersection genes being RP11–95M15.1 and IRF5 from the two analyses ( Supplementary Figure 6C ).