A genome-wide cross-trait analysis identifying shared genetic basis and causal relationships between Hunner-type interstitial cystitis and autoimmune diseases in East Asian populations

Abstract

Purpose:

Epidemiological studies have demonstrated the clinical link between Hunner interstitial cystitis (HIC) and autoimmune diseases (ADs), suggesting potential shared genetic bases for their comorbidity. We aimed to investigate the shared genetic architecture and causal relationships between HIC and ADs.

Methods:

We conducted a genome-wide cross-trait study with ~170000 individuals of East Asian ancestry to investigate the shared architecture between HIC and ADs. Bidirectional Mendelian randomization (MR) was used to assess potential causal relationships and a multi-trait analysis of GWAS (MTAG) was conducted to identify their associated pleiotropic loci. Fine-mapping analysis narrowed candidate gene susceptibility loci and colocalization analysis was performed to identify shared variants at specific locus. Lastly, transcriptome-wide association (TWAS) and functional analysis were utilized to explore potential shared gene-tissue associations.

Results:

Through bidirectional MR analysis, we observed a positive causal effect of AIH(OR_IVW=1.09, P_IVW=1.00×10^-3) and RA (OR_IVW=1.47, P_IVW<1.00×10^-4) on HIC and a negative causal effect of UC on HIC (OR_IVW=0.89, P_IVW< 1.00×10^-4). Furthermore, we unveiled a robust positive causal effect of HIC on T1D(OR_ConMix=1.05, P_ConMix=1.77×10^-3). Cross-trait meta-analysis identified a total of 64 independent SNPs associated with HIC and ADs. Functional analysis revealed that the identified variants regulated gene expression in major tissues belonging to the autoimmune system.

Conclusions:

Our findings might offer insights into the shared underlying etiology of HIC and ADs.

1 Introduction

Hunner interstitial cystitis (HIC) is a rare and challenging chronic inflammatory bladder disease of uncertain etiology. It is characterized by persistent bladder pain and lower urinary tract symptoms (LUTS), such as urinary frequency and urgency, along with the presence of Hunner ulcers on cystoscopy (1). HIC is a subtype of interstitial cystitis (IC)/bladder pain syndrome (BPS), which is a broadly defined chronic pelvic pain syndrome involving the urinary system. IC encompasses various potential etiologies and clinical phenotypes, with a global prevalence of approximately 10.6 cases per 100,000 individuals (2). The diagnosed prevalence in females is about five times higher than in males. Among IC cases, HIC accounts for 3.5% to 50% of all cases. The pathogenic mechanisms underlying HIC remain unclear, although previous studies have suggested that it involves complex interactions between multiple mechanisms, including neural, endocrine, and immune factors (1, 3). Immunoglobulin and complement deposition, aggregation of restricted light-chain plasma cells, and upregulation of pro-inflammatory genes/molecules involved in innate and adaptive immune responses have been detected in the bladder tissues of HIC patients (4–8).

Furthermore, IC may be a systemic disease and is often comorbid with various autoimmune diseases (ADs). Numerous studies have demonstrated an increased prevalence of multiple ADs in IC patients, including rheumatoid arthritis (RA), systemic lupus erythematosus (SLE), Sjögren’s syndrome (SS), inflammatory bowel disease, and autoimmune thyroid diseases (9–12).

However, due to the observational nature of traditional epidemiological studies, methodological limitations still exist in current research on the comorbidity between IC and ADs. The underlying pathological and physiological mechanisms, as well as genetic associations, particularly specific shared genetic factors and potential genetic causal effects between IC and ADs, require further investigation. Genome-wide association studies (GWAS), functional genomics research, and integrative analyses offer new avenues for studying the genetic architecture of complex diseases. These approaches can identify candidate pathogenic genes or tissue/cell types and provide insights into the development of disease-related biomarkers and targeted therapeutics.

In this study, we conducted a large-scale genome-wide cross-trait association study with ~170000 individuals of East Asian ancestry to investigate the shared architecture between HIC and ADs, including atopic dermatitis (AD), autoimmune hepatitis (AIH), allergic rhinitis (AR), asthma (AS), contact dermatitis (CD), Graves’ disease (GD), Hashimoto’s thyroiditis (HT), hypothyroidism (HY), hyperthyroidism (HYPE), myasthenia gravis (MG), pollinosis (PO), psoriasis vulgaris (PV), rheumatoid arthritis (RA), sarcoidosis (SA), systemic lupus erythematosus (SLE), Sjögren’s syndrome (SS), type 1 diabetes mellitus (T1D), ulcerative colitis (UC), uveitis (UV). We aimed to not only assess the genetic correlation and potential causal relationship between HIC and ADs but also to identify pleiotropic loci associated with joint phenotypes. We hope our findings can help better delineate the biological implications of shared genetic architecture between HIC and ADs.

2 Materials and methods

2.1 Study population, design, and data summary

The workflow of our analysis was shown in Figure 1. In brief, there were three main parts in our study: causal inference analysis, cross-trait meta-analysis and post-GWAS analysis between HIC and the 19 autoimmune disorders.

Figure 1

The GWAS summary statistic of HIC includes a total of 153 cases and 46,087 controls (13). The 153 cases were recruited at Tokyo University Hospital in Japan between 2018 and 2020. DNA samples of controls were obtained from the Biobank Japan Project (BBJ). The autoimmune diseases included AD(N_case/N_control=2472/142192, AIH(N_case/N_control=85/166529), AR(N_case/N_control=7897/153666), AS(N_case/N_control=13015/162933), CD(N_case/N_control=247/161777), GD(N_case/N_control=2809/172656), HT(N_case/N_control=537/172656), HY(N_case/N_control=1114/172656), HYPE(N_case/N_control=994/172656), MG(N_case/N_control=81/178630), PO(N_case/N_control=18593/153666), RA(N_case/N_control=5348/173268), SA(N_case/N_control=220/177667), SLE(N_case/N_control=317/175937), SS(N_case/N_control=303/175599), T1D(N_case/N_control=1219/132032), UC(N_case/N_control=314/178375), UV(N_case/N_control=125/174600), the GWAS summary statistics of were obtained from NBDC Human Database, with the Dataset ID hum0197.v3.gwas.v1 (14). We used ANNOVAR (15) to annotate variants of GWAS summary statistics based on ‘hg19 avsnp150’. Additional details for each dataset can be found in Supplementary Table S1 and Supplementary Notes.

2.2 Causal inference analysis

To identify independent genetic instruments, we utilized the PLINK clumping function with the following parameters: clump-p1 = 5e-8, clump-p2 = 0.01, clump_kb=500Kb, and clump_r2 = 0.2. This allowed us to determine the top loci that were independent of each other. To ensure statistical power due to the limited number of instrumental variables (n<10), we employed more lenient criteria for instrument variable selection: clump-p1 = 5e-6, clump-p2 = 0.01, clump_kb=500Kb, and clump_r2 = 0.2. Furthermore, we applied Steiger filtering to the instrumental variables and excluded instruments with F-statistics<10.

We utilized several MR methods to examine the causal relationships between each ADs and HIC. Our primary MR analysis was the contamination mixture (ConMix) approach (16), which explicitly modeled multiple potential causal estimates and inferred multiple causal mechanisms associated with the same risk factor that affects the outcome to different degrees. Additionally, we also applied several sensitivity analyses to validate our results. The MR-PRESSO (17) was employed to remove outliers and ensure efficient use of valid IVs. MR-Egger regression (18) provided estimates after the correction of pleiotropy. The weighted-median (WM) estimator approach, as a median of the weighted estimates, provides a consistent effect even if half of the IVs are pleiotropic (19). The median-based method(MBE) proceeds by constructing a kernel-weighted density of the variant-specific estimates, and taking the maximum point of this density as the point estimate. A confidence interval is obtained by bootstrapping (20). Finally, we employed the inverse-variance weighted (IVW) method (21), which is a robust approach. We corrected multiple testing for MR P-values by the Bonferroni method, and a P-value of 0.00263 (0.05/19) was considered as the significant level. A P-value less than 0.05 is considered as the threshold indicating statistical significance. ADs that showed possible causal relationships with HIC were included in the following analyses.

We also performed several sensitivity analyses to assess the robustness of our results to potential violations of several MR assumptions. a) Heterogeneity was estimated by the Cochran Q test of IVW and MR-Egger; b) The horizontal pleiotropy was estimated using MR-Egger’s intercept; c) The influential outlier IVs due to pleiotropy was identified using MR-PRESSO’s outlier test. After removing the outlier instrumental variables (IVs) identified by MR-PRESSO, we conducted the MR analysis again.

The same approach was taken for the reverse MR which was used to eliminate spurious results due to reverse causation. Generally, all the analyses were conducted using R software 4.2.0. The MR-PRESSO method was performed using the “MRPRESSO” package. The IVW, MR–Egger, WM, ConMix and MBE methods were performed using the “MendelianRandomization” package. The forest plot of single SNP, funnel plot and scatter plot were performed using the “TwoSampleMR” package.

2.3 Fine-mapping credible set analysis

We performed statistical fine-mapping using FINEMAP (22). We computed LD in each locus using R package ‘LDlinkR’ (23)based on genome build GRCh37 and the East Asian population of 1000 Genome project population. We defined a fine-mapping region as the 3Mb (± 1.5Mb) window around each lead variant. This window size is based on recommendations for fine-mapping and colocalization analyses. We allowed up to 10 causal variants per window and extracted the posterior inclusion probabilities (PIP) of each variant using each method independently. The variants with PIP>0.90, along with having LD r²>0.2 with the lead variant, are considered the final candidate causal variants. We applied 3DSNP, a comprehensive database for human noncoding variants annotation, to annotate these causal variants (24).

2.4 Colocalization analysis

We extracted summary statistics for variants within 500 kb(± 250kb) of the index SNP at each of the shared loci between HIC and ADs and performed colocalization analysis between HIC and each ADs trait using R ‘coloc’ package (25) to calculate the probability that the two traits shared a common genetic causal variant. We caculated the posterior probability that the 2 traits were associated with different causal variants(H3) or that the 2 traits were associated and shared 1 common causal variant(H4). The posterior probability for H3(PPH3) or H4(PPH4) that was greater than 0.5 was considered colocalized (26).

2.5 Cross-trait meta-analysis

We then implemented a cross-trait meta-analysis of GWAS summary data using Multi-Trait Analysis of GWAS (MTAG) (27), a method for joint analysis of summary statistics from GWASs of different traits, to identify pleiotropic loci with strong signals associated with ADs and HIC. By analyzing multiple traits together, this approach increases the statistical power of detecting genetic associations for each trait. The MTAG estimator is a variant of the IVW meta-analysis that utilizes summary statistics from single-trait GWASs and generates trait-specific associations statistics. The resulting P-values can be considered as P-values from a single-trait GWAS. We applied PLINK clumping function parameters: -clump-p 5×10^-8, -clump-r2 0.2 -clump-kb 500, to determine top loci. The variant with the lowest p-value was defined as the sentinel variant. A P-value of 5×10^-8 was considered as a genome-wide significance for cross-trait meta-analysis and the significant SNP should also meet a requirement that p value of 5×10^-3 from both single traits.We then performed functional annotation by Functional Annotation of Variants-Online Resource (FAVOR) (28), an open-access variant functional annotation portal for cross-trait meta-analysis.

2.6 eQTL mapping, tissue enrichment analysis, and pathway analysis

To map the shared SNPs between HIC and ADs traits to specific genes which they show a significant eQTL association with, we conducted the eQTL mapping analysis using the Functional Mapping and Annotation(FUMA) website (29), incorporating the SNP2GENE function with the cis-eQTLs reported by DICE (30), which were identified in 13 immune cell types isolated from 106 leukapheresis samples provided by 91 healthy subjects, and the eQTLs reported by van der Wijst et al (31), which were identified from 25,000 peripheral blood mononuclear cells (PBMCs) from 45 donors.

To delve into the biological implications of the shared genes among the ten trait pairs, we performed GTEx tissue enrichment analysis within the clumping region for each trait identified by MTAG, using the FUMA website based on 54 tissue types sourced from GTEx(version 8) (29). Additionally, we utilized the FUMA website to assess the enrichment of independent loci for each trait pairing and to explore shared genes between ADs and HIC, examining their association with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms to elucidate relevant biological pathways. The method for multiple testing correction was BH with an adjusted p-value cutoff (0.05).

2.7 Transcriptome-wide association

To explore the potential shared gene-tissue associations between ADs and HIC, we performed a TWAS using FUSION (R package), based on 49 GTEx (version 8) multi-tissue expression weights. FUSION adopts a Bayesian sparse linear mixed model (BSLMM) (32) that combines Bayesian variable selection (BVSR) (33) and linear mixed model (LMN) (34) with the normal mixture prior assumption to train weights between observed gene expressions and cis-acting genetic variants with reference dataset. This method tests the association between predicted gene expression and phenotypes of interest. Besides, we applied Benjamini-Hochberg correction for each trait’s all gene-tissue pairs on TWAS P-values, accounting for multiple tests (false discovery rate < 0.05).

2.8 Extraction of TPM expression matrix of genes of interest from the GEO database

The microarray datasets GSE11783 (6), GSE55235 (35), GSE206364 (36), and GSE181674 (37) were extracted from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Specifically, the GSE11783 dataset for HIC is based on the GPL570 platform, the GSE55235 dataset for RA is based on the GPL96 platform, the GSE206364 dataset for AIH is based on the GPL20301 platform, and the GSE181674 dataset for T1D is based on the GPL21290 platform. We selected the TPM expression matrix of the genes of interest and performed a Wilcoxon rank-sum test to compare gene expression levels between the control and disease groups. The inclusion details of the samples in each dataset are provided in Supplementary Table S2.

3 Results

3.1 Causal inference between HIC and ADs

We explored the genetic relationship between HIC and ADs using Bidirectional MR analyses to gain additional insights into the genetic connections between these diseases. Bidirectional MR analyses were first addressed to investigate the causal relationship between HIC and ADs. The details were presented in Supplementary Tables S3-4; Figure 1. Based on forward MR analysis, we found that ADs might have a potential causal effect on HIC. Specifically, genetically predicted AIH (OR_IVW=1.09, P_IVW=1.00×10^-3) and RA (OR_IVW=1.47, P_IVW<1.00×10^-4) exhibited significant positive causal effect on HIC. This result was further validated in sensitivity analyses using other MR methods. The forward MR analysis using the ConMix approach revealed potential positive causal effects of SA (OR_ConMix=1.16, P_ConMix=1.63×10^-2) and T1D (OR_ConMix=1.37, P_ConMix=9.22×10^-3) on HIC despite not meeting the threshold of FDR correction. However, the reverse MR analysis revealed significant positive causal effects of HIC on T1D (OR_ConMix=1.05, P_ConMix=1.77×10^-3), suggesting that HIC is likely a risk factor of T1D (Figure 2).

Figure 2

Genetically predicted UC (OR_IVW=0.89, P_IVW< 1.00×10^-4) has a significant negative causal effect on the risk of HIC, which was agreed by other MR methods (P_WM<1.20×10^-4, P_ConMix=3.07×10^-3) and sensitivity analysis using MR-PRESSO (OR_PRESSO=0.89, P_PRESSO=2.23×10^-5). Additionally, a forward MR analysis using the ConMix and Egger approach revealed a potential negative causal effect of GD on HIC (OR_ConMix=0.81, P_ConMix=3.56×10^-2; OR_Egger=0.80, P_ConMix=2.22×10^-2).

Reverse MR analysis revealed a positive causal effect of HIC on RA. However, after correction using the outlier test in MR-PRESSO, the aforementioned causal effects are no longer significant. Additionally, the reverse analysis revealed potential negative causal effects of HIC on HT (OR_ConMix=0.95, P_ConMix=4.98×10^-2) and SLE (OR_ConMix=0.92, P_ConMix=3.58×10^-2) after removing the outliers identified by the PRESSO analysis (Figure 2). The detailed results of the sensitivity analysis for MR can be found in Supplementary Table S5 and Supplementary Figures S1-116.

3.2 Cross-trait meta-analysis between HIC and ADs

After investigating the causal relationships between HIC and ADs, we conducted cross-trait meta-analyses to identify individual SNPs underlying the joint phenotypes based on MTAG(all these SNPs fulfilled P_single<5×10^-3, P_MTAG<5×10^-8) to combine the association evidence for HIC with ADs (Table 1 and Figure 3).

Figure 3

We identified a total of 64 independent SNPs. For HIC and AIH, in total, we identified 1 independent SNPs. The most significant SNP (rs564176274, P_MTAG= 9.40× 10^− 9, P_HIC= 1.62 × 10^− 5, P_AIH= 3.23 × 10^− 5) was located at the intergenic region, which was near gene HLA-DQB1.For HIC and RA, we identified 62 shared independent SNPs. The most significant SNP (rs117530403, P_MTAG= 2.91×10^-131, P_HIC= 3.29×10^-3, P_RA= 2.53×10^-150) was located near gene HLA-DRB1.Notably, there were three significant SNPs located in the exonic regions of HLA-DQB1, HLA-DQB2 and HLA-C. For HIC and T1D, we identified 1 shared independent SNPs. The most significant SNP (rs9268831, P_MTAG= 2.52×10^-8, P_HIC= 1.81×10^-3, P_T1D= 4.26×10^-6) was located near gene HLA-DRA. For HIC and UC, We did not identify any SNP that meets the significance threshold.

3.3 Fine-mapping to identify potential causal variants and colocalization analysis

In fine-mapping potential causal variants underlying HIC and ADs shared signals detected in the GWAS meta-analysis, we nominate 5 putative causal variants at 1 locus, 286 causal variants at 62 loci, 5 causal variants at 1 locus shared between AIH, RA, T1D and HIC, respectively (Supplementary Tables S6-8). Further colocalization analysis was conducted to ascertain whether the genetic variants driving the association in 2 traits are the same or different.We found 41 loci colocalized at different causal variants within 500 kb(± 250kb) of the lead SNP. The only shared loci between HIC and AIH and T1D demonstrated colocalization at distinct causal variants (Supplementary Table S9).

3.4 Tissue-specific enrichment analysis, pathway analysis, and eQTL mapping

The GTEx enrichment analysis independent tissue expression was significantly enriched (FDR < 0.05) for expression of cross-trait associated genes for HIC-RA, which included whole blood, spleen, brain, and small intestine. However, no significant tissue enrichment was identified for HIC-AIH, HIC-UC, and HIC-T1D meta-analysis (Supplementary Figures S117-120).

In terms of GO, we observed that pathways were primarily associated with antigen processing and presentation, immune processes mediated by lymphocytes, leukocytes, B cells, and T cells, production of immunoglobulins, and regulation of cell cytotoxicity, among others. Additionally, in KEGG pathways, we identified enrichment of shared loci between HIC and autoimmune diseases in immune response processes such as antigen processing and presentation, allograft rejection, intestinal immune network for IgA production, natural killer cell-mediated cytotoxicity, as well as autoimmune diseases including graft-versus-host disease, type 1 diabetes mellitus, autoimmune thyroid disease, viral myocarditis, asthma, systemic lupus erythematosus, and leishmania infection(Supplementary Tables S10-13, Supplementary Figures S121-128).

To explore whether the significant SNPs shared between HIC and ADs had a downstream functional impact, we extracted the eQTL data corresponding to the significant SNPs found in each HIC and ADs trait pair from the immune cell eQTL data and the PBMCs eQTL data. Among the 62 shared SNPs between HIC and RA (Supplementary Table S14), we identified 10 SNPs associated with the expression of nearby genes. Notably, increased expression was observed for the symbol genes BTN2A2, BTN3A1, FLOT1, DDR1, HLA-DQA2, HLA-DQB2, HLA-DRB5, HLA-DQA1, and HLA-DPA1. For HIC and T1D (Supplementary Table S15), we identified 1 SNP (rs9262670) associated with the increased expression of HLA-DQA2 and HLA-DQB2 in monocytes, CD4 T cells, CD8 T cells and B cells.

3.5 Shared genes between HIC and ADs from TWAS

We then examined shared TWAS genes between HIC and ADs in specific tissues. After BH correction, we identified 5 TWAS-significant genes between HIC with RA, 4 genes for UC and 2 for T1D (Supplementary Tables S16-18). Notably, the cross-trait meta-analysis of HIC-RA has identified HLA-DOA as genome-wide significant (rs3135196, P_MTAG= 2.44×10^-11, P_HIC= 2.67.00×10^-6, P_RA= 2.41×10^-7). HLA-DPB1 was also significant in the cross-trait meta-analysis of HIC and RA (rs12528890, P_MTAG= 2.01×10^-8, P_HIC= 1.91×10^-7, P_RA= 2.20×10^-4). Notably, although PSMB8 is not significant in the MTAG analysis, both PSMB8 and PSMB9 belong to the gene encoding the proteasome 20S subunit. Additionally, PSMB9 is significant in the MTAG analysis of HIC-RA.

3.6 The transcriptomic expression levels of shared genes between HIC and ADs.

We further validated the differences in the expression of shared genes between the disease and control groups at the transcriptomic level (Figure 4). After excluding pseudogenes and different transcripts, we found that among the shared genes identified through cross-trait meta-analysis, HLA-DQB1, a shared gene between HIC and AIH, was upregulated only in AIH, but the difference was not statistically significant in HIC (Supplementary Figure S129). HLA-DRA was significantly upregulated in HIC, RA, and T1D. Seven shared genes between HIC and RA trait pairs, including BTN3A1, C2, HLA-C, HLA-DQA1, PSMB9, MICB, and TAPBP, were significantly upregulated in both HIC and RA (P<0.05). However, HLA-DOA was upregulated only in HIC, while HLA-DQB2, HLA-DQB1, and HCP5 were upregulated only in RA. Additionally, ABCF1 and TNXB were significantly downregulated only in RA.

Figure 4

4 Discussion

We conducted, to the best of our knowledge, the first comprehensive assessment of the shared genetics between HIC and ADs in the East Asian population by analyzing large-scale GWAS summary data using multiple statistical genetic approaches. Through bidirectional MR analysis, we investigated the causal relationships between HIC and several ADs. Notably, we observed a significant positive causal effect of AIH and RA on HIC, which is consistent with previous epidemiological findings. A nationwide population-based study conducted in Taiwan also indicated an association between IC/BPS and the development of RA. Furthermore, our study unveiled a robust positive causal effect of HIC on T1D. In contrast, in the MR analysis, we observed a significant negative causal effect of UC on HIC, which diverges from previous epidemiological reports (12).

Moreover, several ADs, including GD, SA, and SLE, demonstrated statistical significance in the bidirectional MR analysis with HIC. However, these results did not pass the multiple corrections. Interestingly, although previous epidemiological reports have suggested a common comorbidity between Sjögren’s syndrome and IC, we did not observe a significant correlation in our causal inference. This discrepancy might be due to the relatively small sample size of GWAS studies on Sjögren’s syndrome.

We also identified 64 independent loci shared between HIC and AIH, RA, and T1D at genome-wide significant level. We highlighted HLA region (several sentinel SNPs) for its significant role in between HIC and ADs. HLA region harbors more than 200 genes located close to each other on chromosome 6, one of the most extensively studies regions in human genome that contains abundant pleiotropy for many complex diseases, especially involved in the immune-related process (38). The genome-wide association study has identified that three amino acid positions in human leukocyte antigen HLA-DQB1 and one amino acid position in HLA-DPB1 were associated with the increased risk of HIC, which revealed that genetic contributions to HIC risk that may be associated with class II MHC molecule antigen presentation. Notably, there is no significant difference in the transcriptional levels of HLA-DQB1 and HLA-DPB1 between the HIC and control groups, suggesting that the associated risk SNPs may contribute to disease by affecting protein function rather than regulating gene transcription levels. Furthermore, Tseng et al. compared global gene expression profiles in bladder epithelial cells between patients with HIC and normal controls, and observed upregulations of major histocompatibility complex (MHC) class I(HLAF) and class II (HLA-F, HLA-DQB1, HLA-DRB1, HLA-DPA1, HLA-DOA, HLA-DMA and HLA-DRA) molecules in bladder epithelial from IC and ulcerative IC area.

In this study, we also confirmed that HLA-DQB1 is the most significant locus in the MTAG analysis of HIC-AIH. Additionally, we have uncovered distinct potential susceptibility loci between HIC and different ADs, with the majority of these loci being located within the HLA region mentioned above. Specifically, we emphasized several genes, including overlapping loci across trait pairs such as HLA-DRB1, HLA-DRA and two genes that were significant in both MTAG and TWAS analysis, namely HLA-DOA and PSMB9.

HLA-DRA is the most significant loci shared between HIC and T1D, and it also exhibited significance between HIC and RA. HLA-DRA is one of the HLA class II alpha chain paralogues, playing a crucial role in antigen presentation (39). Notably, DRA lacks polymorphisms in the peptide binding region and acts as the sole alpha chain for DRB1, DRB3, DRB4, and DRB5. The HLA-DR4 (DRB1*0405–DQB1*0401) and HLA-DR9 (DRB1*0901–DQB1*0303) haplotypes were primarily associated with T1D in East Asian populations (40). Genetic variation at the HLA-DRB1 gene is associated with RA, and HLA-DRB1*0404, *0405 have also been found to exhibit strong associations with RA in Asians (41).

HLA-DOA, a member of the HLA class II alpha chain paralogues, exhibited a significant association between HIC and RA, as well as in TWAS analyses of HIC-RA, HIC-T1D and HIC-UC. HLA-DOA forms a heterodimer with HLA-DOB. The heterodimer, HLA-DO, is localized in lysosomes of B cells and regulates HLA-DM-mediated peptide loading on MHC class II molecules (42). Okada et al. found that HLA-DOA, a non-classical HLA gene, was an independent risk factor on ACPA-positive RA and demonstrated a cis-eQTL effect of the causal variant in Japanese population (43). Furthermore, HLA-DOA has also been confirmed as a susceptibility locus within MHC with a moderate contribution to T1D that is independent of HLA-DRB1 locus (44).

PSMB9 and PSMB8 encode beta subunits of the proteasome and are located within the class II region of the MHC. The proteasome is primarily responsible for cleaving class I MHC peptides in an ATP/ubiquitin-dependent process within a non-lysosomal pathway (45). Li et al. identified PSMB9 as a diagnostic marker for RA using a machine learning approach (46). Furthermore, we also observed AL662789.1 and AL671883.3, which are human DNA sequences derived from clone XXbac-254F23 on chromosome 6 and clone CH501-248L24 on chromosome 6, respectively. These sequences were found to be shared by more than one pair of HIC and ADs. However, there is currently limited research on these genes.

This study possesses several notable strengths. Firstly, most of the cross-trait studies now focus on the European population, and this is the first analysis to identify the shared genetic architecture of HIC and ADs using a large-scale observational GWAS dataset consisting exclusively of East Asian samples, after identifying ten GWAS sources (details were shown in Supplementary Table S1). Furthermore, we utilized multi-omics statistical methods such as MTAG and TWAS to identify novel genes and pathways associated with both HIC and ADs. The novel genes we discovered may serve as potential drug targets for the treatment of the disease, although further validation is required. Notably, they could offer new diagnostic and therapeutic avenues for patients with HIC, particularly those with comorbid autoimmune disorders.

We would like to acknowledge several potential limitations in our study. Firstly, the limited sample size of participants with each mental disorder and the restriction to individuals of East Asian ancestry have limited the statistical power of our GWAS analysis, which may limit the generalizability of our findings to other ancestral populations. Secondly, the lack of GWAS data for other subtypes of IC, such as non-ulcerative IC, has hindered our ability to investigate and differentiate the genetic associations between different subtypes of interstitial cystitis and immune-related disorders. Thirdly, it is important to consider common non-genetic risk factors for the occurrence of hypersensitivity HIC and ADs, such as medication and environmental factors. Our current study focused solely on evaluating shared genetic factors between HIC and ADs, and future research should investigate shared environmental factors between these conditions.

5 Conclusion

Our study provides evidence of a genetic correlation and causality between HIC and ADs. We identified genetic loci associated with both HIC and autoimmune disorders, as well as potential causal relationships between disease trait pairs, thereby enriching our understanding of HIC and shedding light on the shared genetic etiology of HIC and autoimmune disorders. In addition, we discovered multiple potential common biological mechanisms that can enhance our knowledge of the link between HIC and ADs. These discoveries open up new avenues for future research on functional validation, disease prevention, and clinical treatment strategies.

References

• 1

  HommaYAkiyamaYTomoeHFurutaAUedaTMaedaDet al. Clinical guidelines for interstitial cystitis/bladder pain syndrome. Int J Urol. (2020) 27:578–89. doi: 10.1111/iju.14234

  • CrossRef
  • Google Scholar

• 2

  PatnaikSSLaganàASVitaleSGButticèSNoventaMGizzoSet al. Etiology, pathophysiology and biomarkers of interstitial cystitis/painful bladder syndrome. Arch Gynecol Obstet. (2017) 295:1341–59. doi: 10.1007/s00404-017-4364-2

  • CrossRef
  • Google Scholar

• 3

  van de MerweJPNordlingJBouchelouchePBoucheloucheKCervigniMDahaLKet al. Diagnostic criteria, classification, and nomenclature for painful bladder syndrome/interstitial cystitis: an ESSIC proposal. Eur Urol. (2008) 53:60–7. doi: 10.1016/j.eururo.2007.09.019

  • CrossRef
  • Google Scholar

• 4

  OchsRLSteinTWPeeblesCLGittesRFTanEM. Autoantibodies in interstitial cystitis. J Urol. (1994) 151:587–92. doi: 10.1016/S0022-5347(17)35023-1

  • CrossRef
  • Google Scholar

• 5

  MaedaDAkiyamaYMorikawaTKunitaAOtaYKatohHet al. Hunner-type (Classic) interstitial cystitis: A distinct inflammatory disorder characterized by pancystitis, with frequent expansion of clonal B-cells and epithelial denudation. PloS One. (2015) 10:e0143316. doi: 10.1371/journal.pone.0143316

  • CrossRef
  • Google Scholar

• 6

  GamperMViereckVGeissbühlerVEberhardJBinderJMollCet al. Gene expression profile of bladder tissue of patients with ulcerative interstitial cystitis. BMC Genomics. (2009) 10:199. doi: 10.1186/1471-2164-10-199

  • CrossRef
  • Google Scholar

• 7

  AkiyamaYMaedaDKatohHMorikawaTNiimiANomiyaAet al. Molecular taxonomy of interstitial cystitis/bladder pain syndrome based on whole transcriptome profiling by next-generation RNA sequencing of bladder mucosal biopsies. J Urol. (2019) 202:290–300. doi: 10.1097/JU.0000000000000234

  • CrossRef
  • Google Scholar

• 8

  GamperMViereckVEberhardJBinderJMollCWelterJet al. Local immune response in bladder pain syndrome/interstitial cystitis ESSIC type 3C. Int Urogynecol J. (2013) 24(12):2049–57. doi: 10.1007/s00192-013-2112-0

  • CrossRef
  • Google Scholar

• 9

  YuehHZYangMHHuangJYWeiJCC. Risk of autoimmune diseases in patients with interstitial cystitis/bladder pain syndrome: A nationwide population-based study in Taiwan. Front Med (Lausanne). (2021) 8:747098. doi: 10.3389/fmed.2021.747098

  • CrossRef
  • Google Scholar

• 10

  UedaYTomoeHTakahashiHTakahashiYYamashitaHKanekoHet al. Interstitial cystitis associated with primary Sjögren’s syndrome successfully treated with a combination of tacrolimus and corticosteroid: A case report and literature review. Mod Rheumatol. (2016) 26:445–9. doi: 10.3109/14397595.2014.895283

  • CrossRef
  • Google Scholar

• 11

  KellerJJLiuSPLinHC. A case-control study on the association between rheumatoid arthritis and bladder pain syndrome/interstitial cystitis. Neurourol and Urodynamics. (2013) 32(7):980–5.

  • Google Scholar

• 12

  de MerweJPYamadaTSakamotoY. Systemic aspects of interstitial cystitis, immunology and linkage with autoimmune disorders. Int J Urol. (2003), 10 Suppl:S35–38. doi: 10.1046/j.1442-2042.10.s1.10.x

  • CrossRef
  • Google Scholar

• 13

  AkiyamaYSoneharaKMaedaDKatohHNaitoTYamamotoKet al. Genome-wide association study identifies risk loci within the major histocompatibility complex region for Hunner-type interstitial cystitis. Cell Rep Med. (2023) 4:101114. doi: 10.1016/j.xcrm.2023.101114

  • CrossRef
  • Google Scholar

• 14

  SakaueSKanaiMTanigawaYKarjalainenJKurkiMKoshibaSet al. A cross-population atlas of genetic associations for 220 human phenotypes. Nat Genet. (2021) 53:1415–24. doi: 10.1038/s41588-021-00931-x

  • CrossRef
  • Google Scholar

• 15

  WangKLiMHakonarsonH. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. (2010) 38:e164. doi: 10.1093/nar/gkq603

  • CrossRef
  • Google Scholar

• 16

  BurgessSFoleyCNAllaraEStaleyJRHowsonJMM. A robust and efficient method for Mendelian randomization with hundreds of genetic variants. Nat Commun. (2020) 11:376. doi: 10.1038/s41467-019-14156-4

  • CrossRef
  • Google Scholar

• 17

  VerbanckMChenCYNealeBDoR. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet. (2018) 50:693–8. doi: 10.1038/s41588-018-0099-7

  • CrossRef
  • Google Scholar

• 18

  BowdenJDavey SmithGBurgessS. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol. (2015) 44:512–25. doi: 10.1093/ije/dyv080

  • CrossRef
  • Google Scholar

• 19

  BowdenJDavey SmithGHaycockPCBurgessS. Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol. (2016) 40:304–14. doi: 10.1002/gepi.2016.40.issue-4

  • CrossRef
  • Google Scholar

• 20

  HartwigFPDavey SmithGBowdenJ. Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption. Int J Epidemiol. (2017) 46:1985–98. doi: 10.1093/ije/dyx102

  • CrossRef
  • Google Scholar

• 21

  YuanSLarssonSC. Assessing causal associations of obesity and diabetes with kidney stones using Mendelian randomization analysis. Mol Genet Metab. (2021) 134:212–5. doi: 10.1016/j.ymgme.2021.08.010

  • CrossRef
  • Google Scholar

• 22

  BennerCSpencerCCHavulinnaASSalomaaVRipatiSPirinenM. FINEMAP: efficient variable selection using summary data from genome-wide association studies. Bioinf (Oxford England). (2016) 32(10):1493–501. doi: 10.1093/bioinformatics/btw018

  • CrossRef
  • Google Scholar

• 23

  MyersTAChanockSJMachielaMJ. LDlinkR: an R package for rapidly calculating linkage disequilibrium statistics in diverse populations. Front Genet. (2020) 11:157. doi: 10.3389/fgene.2020.00157

  • CrossRef
  • Google Scholar

• 24

  QuanCPingJLuHZouGLuY. 3DSNP 2.0: update and expansion of the noncoding genomic variant annotation database. Nucleic Acids Res. (2022) 50(D1):D950–5.

  • Google Scholar

• 25

  GiambartolomeiCVukcevicDSChadtEEFrankeLHingoraniADWallaceCet al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PloS Genet. (2014) 10:e1004383. doi: 10.1371/journal.pgen.1004383

  • CrossRef
  • Google Scholar

• 26

  ChenDWangXJiaJHuangT. Sleep and alzheimer’s disease: shared genetic risk factors, drug targets, molecular mechanisms, and causal effects. Front Genet. (2022) 13:794202. https://www.researchsquare.com/article/rs-853181/v1.

  • Google Scholar

• 27

  TurleyPWaltersRKMaghzianOOkbayALeeJJFontanaMAet al. Multi-trait analysis of genome-wide association summary statistics using MTAG. Nat Genet. (2018) 50:229–37. doi: 10.1038/s41588-017-0009-4

  • CrossRef
  • Google Scholar

• 28

  ZhouHArapoglouTLiXLiZZhengXMooreJet al. FAVOR: functional annotation of variants online resource and annotator for variation across the human genome. Nucleic Acids Res. (2023) 51:D1300–11. doi: 10.1093/nar/gkac966

  • CrossRef
  • Google Scholar

• 29

  WatanabeKTaskesenEvan BochovenAPosthumaD. Functional mapping and annotation of genetic associations with FUMA. Nat Commun. (2017) 8:1826. doi: 10.1038/s41467-017-01261-5

  • CrossRef
  • Google Scholar

• 30

  SchmiedelBJSinghDMadrigalAValdovino-GonzalezAGWhiteBMZapardiel-GonzaloJet al. Impact of genetic polymorphisms on human immune cell gene expression. Cell. (2018) 175:1701–1715.e16. doi: 10.1016/j.cell.2018.10.022

  • CrossRef
  • Google Scholar

• 31

  van der WijstMGPBruggeHde VriesDHDeelenPSwertzMAFrankeL. Single-cell RNA sequencing identifies cell type-specific cis-eQTLs and co-expression QTLs. Nat Genet. (2018) 50:493–7. doi: 10.1038/s41588-018-0089-9

  • CrossRef
  • Google Scholar

• 32

  ZhouXCarbonettoPStephensM. Polygenic modeling with bayesian sparse linear mixed models. PloS Genet. (2013) 9:e1003264. doi: 10.1371/journal.pgen.1003264

  • CrossRef
  • Google Scholar

• 33

  GuanYStephensM. Bayesian variable selection regression for genome-wide association studies and other large-scale problems. Ann Appl Statistics. (2011) 5:1780–815, 36. doi: 10.1214/11-AOAS455

  • CrossRef
  • Google Scholar

• 34

  YuJPressoirGBriggsWHVroh BiIYamasakiMDoebleyJFet al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat Genet. (2006) 38:203–8. doi: 10.1038/ng1702

  • CrossRef
  • Google Scholar

• 35

  WoetzelDHuberRKupferPPohlersDPfaffMDrieschDet al. Identification of rheumatoid arthritis and osteoarthritis patients by transcriptome-based rule set generation. Arthritis Res Ther. (2014) 16:R84. doi: 10.1186/ar4526

  • CrossRef
  • Google Scholar

• 36

  KriegermeierAHyonALeCuyerBHubchakSLiuXGreenRM. Inositol-requiring enzyme 1α/X-box protein 1 pathway expression is impaired in pediatric cholestatic liver disease explants. PloS One. (2022) 17:e0279016. doi: 10.1371/journal.pone.0279016

  • CrossRef
  • Google Scholar

• 37

  Campbell-ThompsonMButterworthEABoatwrightJLNairMANasifLHNasifKet al. Islet sympathetic innervation and islet neuropathology in patients with type 1 diabetes. Sci Rep. (2021) 11:6562. doi: 10.1038/s41598-021-85659-8

  • CrossRef
  • Google Scholar

• 38

  SivakumaranSAgakovFTheodoratouEPrendergastJGZgagaLManolioTet al. Abundant pleiotropy in human complex diseases and traits. Am J Hum Genet. (2011) 89:607–18. doi: 10.1016/j.ajhg.2011.10.004

  • CrossRef
  • Google Scholar

• 39

  MaternBMOlieslagersTIVoorterCEMGroenewegMTilanusMGJ. Insights into the polymorphism in HLA-DRA and its evolutionary relationship with HLA haplotypes. HLA. (2020) 95:117–27. doi: 10.1111/tan.13730

  • CrossRef
  • Google Scholar

• 40

  KawabataYIkegamiHKawaguchiYFujisawaTShintaniMOnoMet al. Asian-specific HLA haplotypes reveal heterogeneity of the contribution of HLA-DR and -DQ haplotypes to susceptibility to type 1 diabetes. Diabetes. (2002) 51:545–51. doi: 10.2337/diabetes.51.2.545

  • CrossRef
  • Google Scholar

• 41

  LiuWXJiangYHuQXYouXB. HLA-DRB1 shared epitope allele polymorphisms and rheumatoid arthritis: a systemic review and meta-analysis. Clin Invest Med. (2016) 39:E182–203. doi: 10.25011/cim.v39i6.27487

  • CrossRef
  • Google Scholar

• 42

  PetersdorfEWStevensonPMalkkiMStrongRKSpellmanSRHaagensonMDet al. Patient HLA germline variation and transplant survivorship. J Clin Oncol. (2018) 36:2524–31. doi: 10.1200/JCO.2017.77.6534

  • CrossRef
  • Google Scholar

• 43

  OkadaYSuzukiAIkariKTeraoCKochiYOhmuraKet al. Contribution of a non-classical HLA gene, HLA-DOA, to the risk of rheumatoid arthritis. Am J Hum Genet. (2016) 99:366–74. doi: 10.1016/j.ajhg.2016.06.019

  • CrossRef
  • Google Scholar

• 44

  SantinICastellanos-RubioAAransayAMGutierrezGGaztambideSRicaIet al. Exploring the diabetogenicity of the HLA-B18-DR3 CEH: independent association with T1D genetic risk close to HLA-DOA. Genes Immun. (2009) 10:596–600. doi: 10.1038/gene.2009.41

  • CrossRef
  • Google Scholar

• 45

  KimMSerwaRASamlukLSuppanzIKodrońAStępkowskiTMet al. Immunoproteasome-specific subunit PSMB9 induction is required to regulate cellular proteostasis upon mitochondrial dysfunction. Nat Commun. (2023) 14:4092. doi: 10.1038/s41467-023-39642-8

  • CrossRef
  • Google Scholar

• 46

  LiZChenYZulipikaerMXuCFuJDengTet al. Identification of PSMB9 and CXCL13 as immune-related diagnostic markers for rheumatoid arthritis by machine learning. Curr Pharm Des. (2022) 28:2842–54. doi: 10.2174/1381612828666220831085608

  • CrossRef
  • Google Scholar