Detection of Genetic Overlap Between Rheumatoid Arthritis and Systemic Lupus Erythematosus Using GWAS Summary Statistics.

BACKGROUND: Clinical and epidemiological studies have suggested systemic lupus erythematosus (SLE) and rheumatoid arthritis (RA) are comorbidities and common genetic etiologies can partly explain such coexistence. However, shared genetic determinations underlying the two diseases remain largely unknown.
METHODS: Our analysis relied on summary statistics available from genome-wide association studies of SLE (N = 23,210) and RA (N = 58,284). We first evaluated the genetic correlation between RA and SLE through the linkage disequilibrium score regression (LDSC). Then, we performed a multiple-tissue eQTL (expression quantitative trait loci) weighted integrative analysis for each of the two diseases and aggregated association evidence across these tissues via the recently proposed harmonic mean P-value (HMP) combination strategy, which can produce a single well-calibrated P-value for correlated test statistics. Afterwards, we conducted the pleiotropy-informed association using conjunction conditional FDR (ccFDR) to identify potential pleiotropic genes associated with both RA and SLE.
RESULTS: We found there existed a significant positive genetic correlation (r g = 0.404, P = 6.01E-10) via LDSC between RA and SLE. Based on the multiple-tissue eQTL weighted integrative analysis and the HMP combination across various tissues, we discovered 14 potential pleiotropic genes by ccFDR, among which four were likely newly novel genes (i.e., INPP5B, OR5K2, RP11-2C24.5, and CTD-3105H18.4). The SNP effect sizes of these pleiotropic genes were typically positively dependent, with an average correlation of 0.579. Functionally, these genes were implicated in multiple auto-immune relevant pathways such as inositol phosphate metabolic process, membrane and glucagon signaling pathway.
CONCLUSION: This study reveals common genetic components between RA and SLE and provides candidate associated loci for understanding of molecular mechanism underlying the comorbidity of the two diseases.

Introduction

Rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE) are two frequent chronic rheumatic autoimmune diseases (Suresh, 2004; Cojocaru et al., 2011; Hu et al., 2019; Barbeira et al., 2020), with an incidence rate of 50 (Alamanos et al., 2006) or 8.6 (Stojan and Petri, 2018) per 100,000/year for RA or SLE in the United States, respectively, leading to growing disease burden worldwide among different ethnic, racial and age groups. Both RA and SLE have complicated etiologies, with genetic and hormonal factors (e.g., sex hormones, Tedeschi et al., 2013), environmental factors (e.g., cigarette smoking, Costenbader et al., 2004) and their interaction contributing to the onset and pathological processes of the two diseases (Salaman, 2003). Prior clinical and epidemiological studies have demonstrated that RA and SLE have overlapping clinical symptoms and increased familial aggregation (Alarcón-Segovia et al., 2005; Michou et al., 2008; Icen et al., 2009; Zhernakova et al., 2013; Acosta-Herrera et al., 2019; Ciccacci et al., 2019; James et al., 2019; Xiao et al., 2020), implying common susceptible mechanism and shared predisposition underlying two diseases.

One hypothesis to explain such comorbidity between RA and SLE is the common genetic etiology among autoimmune diseases (Cotsapas et al., 2011; Orozco et al., 2011; Suzuki et al., 2016; Márquez et al., 2017; Acosta-Herrera et al., 2019; Lehallier et al., 2019). Prior studies also provided evidence for shared genetic architecture of RA and SLE. For example, RA patients often have a higher incidence of HLA-DR4 genotype (chr6) compared to healthy controls (Stastny, 1978); meanwhile, correlation was well established between HLA class III (especially 6p21.3) and the susceptibility to develop SLE (Goldberg et al., 1976). Recently, genome-wide association studies (GWASs) have greatly advanced our knowledge of polygenic etiology of RA and SLE, and discovered a number of shared single nucleotide polymorphisms (SNPs) associated with them (Supplementary Table 1; Deng and Tsao, 2017; Okada et al., 2017). Examples of such shared loci include rs7574865 in STAT4 (Remmers et al., 2007), rs2230926 and rs10499194 in TNFAIP3 (Shimane et al., 2010), rs2476601 in PTPN22 (Orozco et al., 2005), rs5754217 in UBE2l3 (Orozco et al., 2011), rs9603612 nearby COG6 (Márquez et al., 2017), as well as multiple loci in NAB1, KPNA4-ARL14, DGQK, LIMK1, and PRR12 (Acosta-Herrera et al., 2019). Understanding these shared genetic determinants has significant implication for identifying important biomarkers and possesses the potential to develop novel therapeutic strategies for joint prediction, prevention, and intervention of RA and SLE.

However, the causal genes and pathways of RA and SLE remain largely unknown because, like many other diseases/traits (Manolio et al., 2009; Orozco et al., 2010), RA- or SLE-associated SNPs identified by GWASs explain only a very small fraction of phenotypic variance of RA or SLE (Julià Cano, 2011; Chen et al., 2017), suggesting that a large number of genetic variants with small to modest effect sizes (but still important) have yet not been discovered and that more pleiotropic genes would be discovered if increasing sample sizes (Wang et al., 2005; Tam et al., 2019). But increasing sample size is generally not feasible since recruiting and genotyping additional participants are expensive and time consuming. In addition, analyzing RA and SLE jointly with individual-level dataset is also difficult because of privacy concerns on data sharing (Pasaniuc and Price, 2016). Instead, summary statistics from large scale GWASs of RA and SLE are freely accessible (Okada et al., 2014; Bentham et al., 2015).

Therefore, a promising way is to apply genetic computational methods that efficiently analyze information contained in the existing pool of available GWAS summary statistics of RA and SLE for unveiling shared genetic contributors with pleiotropic effects more comprehensively. However, few of such studies on common genetic backgrounds of RA and SLE have been undertaken so far. To fill this literature gap, relying on summary statistics obtained from GWASs for SLE (N = 23,210) and RA (N = 58,284) (Okada et al., 2014; Bentham et al., 2015), in the present work we first evaluated the genetic correlation between the two diseases with cross-disease linkage disequilibrium score regression (LDSC) to quantify the extent of common genetic basis to which they share (Bulik-Sullivan et al., 2015). Next, we conducted an multiple-tissue eQTL (expression quantitative trait loci) weighted integrative analysis to aggregate the evidence of SNP-level associations into an integrative association at the gene level and then applied the recently proposed harmonic mean P-value (HMP) combination strategy (Wilson, 2019a) to combine a set of correlated P-values across various tissues into a single well-calibrated P-value. We further preformed the pleiotropy-informed association method using conditional false discovery rate (cFDR) (Andreassen et al., 2013; Smeland et al., 2020) to detect potentially pleiotropic genes. In total, we identified 14 genes that were associated with both RA and SLE, with four of them being likely newly novel pleiotropic genes. The flowchart of our data analysis is demonstrated in Figure 1.

FIGURE 1

Flowchart of data preparation and analysis for RA and SLE. First, two sets of summary statistics were included for the two diseases; a series of quality control procedures were implemented and the genetic correlation between RA and SLE was evaluated using LDSC. Next, a multiple-tissue eQTL weighted integrative analysis was performed to aggregate association signals at the SNP level into the gene level, following by the HMP combination across tissues. Finally, cFDR was carried out to identify associated genes with pleiotropic effects. In the analysis, the LD was estimated with genotypes of the 1000 Genomes Project.

Materials and Methods

GWAS Dataset for RA and SLE

We downloaded summary statistics (e.g., effect allele, effect size and P-value) of RA and SLE from public portal. Specifically, the GWAS of RA included 58,284 (14,361 cases and 43,923 controls) European individuals and 6,446,682 SNPs (Okada et al., 2014), while the GWAS of SLE included 23,210 (7,219 cases and 15,991 controls) individuals of European descent (Bentham et al., 2015) and 7,915,251 SNPs. Based on these summary statistics, we attempted to explore genetic overlap between RA and SLE at the gene level using novel statistical genetics approaches.

Estimation of Cross-Disease Genetic Correlation

We first applied cross-disease LDSC to assess the genetic correlation r between RA and SLE (Bulik-Sullivan et al., 2015). The software of LDSC (version v1.0.1) was available at https://github.com/bulik/ldsc and implemented with default parameter settings. Following prior studies (Bulik-Sullivan et al., 2015), we performed stringent quality control before the LDSC analysis: (i) removed non-biallelic SNPs and those with strand-ambiguous alleles; (ii) excluded duplicated SNPs and those having no rs labels; (iii) excluded SNPs that were located within two special genetic regions including major histocompatibility complex (chr6: 28.5-33.3 Mb) (Bulik-Sullivan et al., 2015) and chr8: 7.2-12.5 Mb (Price et al., 2008) due to their complicated LD structure; (iv) kept SNPs that were included in the 1000 Genomes Project phase III (1000G); (v) removed SNPs whose allele did not match that in the 1000G.

The LD scores were computed using genotypes of 4,098,768 common SNPs (minor allele frequency > 0.01 and the P-value of Hardy Weinberg equilibrium test > 1E-5) with a 10 Mb window on 503 European individuals in the 1000G (Consortium, 2015); and then regressed them on the product of Z score statistics of RA and SLE. The regression slope of LDSC provides an unbiased estimate for r even when the samples are overlapped between the two GWASs of diseases (Bulik-Sullivan et al., 2015).

Association Analysis by Integrating eQTL and GWAS Summary Statistics

Unlike prior studies which explored genetic overlap at the independent SNP level by using a pruning procedure (Lv et al., 2017; Peng et al., 2017; Hu et al., 2018a,b), we attempted to study common genetic component between RA and SLE at the gene level because gene is a more meaningful biological unit related to complex diseases compared with SNP. To do so, we performed the multiple-tissue eQTL weighted integrative analysis for a set of cis-SNPs located within a gene and produced a single P-value for the evidence of the significance of that gene (Gusev et al., 2016; Xu et al., 2017; Guo and Wu, 2018; Wu and Pan, 2018; Xue et al., 2020; Zhang et al., 2020). Specifically, for each tissue in turn and a set of predefined cis-SNPs of a gene of focus, we have:

where Z = (Z1, …, Z) is an m-vector of the Z score for cis-SNPs obtained from summary statistics, with m being the number of cis-SNPs and varying gene by gene across the whole genome; w is an m-vector of cis-SNP weights yielded from the k tissue of GTEx; R is the unknown LD among cis-SNPs and can be approximately estimated from reference panels such as the 1000G (Consortium, 2015). We performed our eQTL-weighted integrative analysis using the metaXcan software (Barbeira et al., 2018). For each gene its cis-SNPs were previously annotated by the authors of metaXcan and the eQTL weights were also in prior trained using the elastic net model with genotypes and gene expressions in tissues from the GTEx Project (GTEx Consortium, 2015). We downloaded tissue-specific eQTL weights from http://predictdb.org/ and run the integrative analysis in terms of the guideline of metaXcan (Barbeira et al., 2018).

Because both the RA and SLE GWASs were analyzed with pooled samples of male and female individuals; therefore, to avoid the influence of gender heterogeneity in gene expressions on association signals (Kassam et al., 2019; Lopes-Ramos et al., 2020; Oliva et al., 2020), we removed six gender-specific tissues in GTEx (i.e., breast mammary tissue, ovary, prostate, testis, uterus, and vagina), leaving 42 various gender-combined GTEx tissues for each disease (Supplementary Table 2). Moreover, due to population stratification or cryptic relatedness or overcorrection of test statistics (Cardon and Palmer, 2003; Freedman et al., 2004; van den Berg et al., 2019), the empirical null distribution in GWAS is sometimes inflated. To correct for such bias, following prior work (Price et al., 2010), we divided by the genomic inflation factor (λ) if it was greater than 1.05. We estimated λ as the ratio between the observed median of and the expected value of 0.456 (Devlin and Roeder, 1999; Dadd et al., 2009; Zeng et al., 2015b). Then, the P-value was easily calculated as the corrected test statistics asymptotically follow a chi-squared distribution with one degree of freedom (Gusev et al., 2016; Xu et al., 2017; Guo and Wu, 2018; Wu and Pan, 2018; Xue et al., 2020; Zhang et al., 2020). Afterwards, we obtained a set of various P-values (P, k = 1, 2, …, 42) for every gene across these tissues, with each representing the association significance of the gene associated with RA or SLE after integrating eQTLs.

To aggregate individual association evidence across tissues, we further applied the HMP combination method to generate a single well-calibrated P-value (Wilson, 2019a):

where ω represents the non-negative weight for each P with ∑_k  =  1^Kω_k   =  1 and assume that ω is independent of P; f denotes the Landau distribution probability density function. Note that, individual Ps often exhibited non-negligible positive dependence because they were implemented for the same gene following the similar logic (see below), existing methods such as the minimum P-value method (Conneely and Boehnke, 2007; Sun and Lin, 2020) and Fisher’s combination (Fisher, 1934), are either computationally intensive or rather difficult to implement because of the requirement of sampling-based algorithms for the valid null distribution or the unavailability of correlation structure (Ballard et al., 2010; Liu and Xie, 2019a,b; Zhang et al., 2019; Sun and Lin, 2020), especially when only summary-level datasets were available (Pasaniuc and Price, 2016). The advantage of HMP used here is that it has been theoretically demonstrated that the complicated positive dependency among P-values has little influence on the final pooled P-value (Wilson, 2019a). That is, T in the Eq. (2) still follows a Landau distribution asymptotically regardless of the correlation structure among these P-values. Consequently, we can yield the P-value for the test statistic T based on the right tail area of the Landau distribution as shown in (2). It has been also proven that under regularity conditions for the generalized central limit theorem the combined P-value by HMP is robust against the number of tests K and the selected weights (Wilson, 2019a). We implemented HMP with equal weights through the harmonic mean p package (version 3.0) in R (Wilson, 2019b). Using the HMP procedure, we generated two sets of P-values in the eQTL weighted gene-based association analyses of RA and SLE.

Pleiotropy-Informed Method Identifying Shared Genes Between RA and SLE

Finally, to leverage the pleiotropic information shared between RA and SLE to efficiently identify gene association signals, we utilized the cFDR method (Andreassen et al., 2013, 2014; Smeland et al., 2020) which extended the unconditional FDR (Benjamini and Hochberg, 1995) from an empirical Bayes perspective. The cFDR measures the probability of the association of the principal disease (e.g., RA) conditioned on the strength of association with the conditional disease (e.g., SLE):

where p and p are the observed HMP adjusted P-values of a particular gene of the principal disease (denoted by i) and the conditional disease (denoted by j), respectively; denotes the null hypothesis that there does not exist association between the gene and the principal disease. As the principal and conditional positions of the two diseases in cFDR are exchangeable, cFDR(p| p) is defined in a similar manner. Moreover, the conjunction conditional false discovery rate (ccFDR) is applied to identify genes with pleiotropic effect:

which is defined as the probability that a given gene has a false positive association with both the principal and conditional diseases, and provides an indicator for pleiotropy (Andreassen et al., 2013; Smeland et al., 2020).

Results

Estimated Genetic Correlation Between RA and SLE

After quality control, a total of 4,708,829 and 5,462,109 genetic variants were reserved for RA and SLE, respectively. The genome-wide SNP-based heritability is estimated to be 11.8% (se = 1.3%) for RA and 30.1% (se = 3.1%) for SLE with LDSC. Next, using cross-disease LDSC with 4,098,768 common SNPs, we observe that there exists a significantly positive genetic correlation between the two types of diseases (r = 0.404, P = 6.01E-10), providing statistical evidence supporting overlapped genetic foundation between RA and SLE. Overall, through the above genetic correlation analysis we reveal RA and SLE are genetically similar and share moderate to high overlap in genetic etiology. Therefore, it is worthy of further investigation into common genetic mechanisms through novel pleiotropy-informed statistical tools.

Associated Genes Identified by cFDR

In our gene-based association analysis, we assigned a set of genetic variants to predefined genes and obtained a total of 23,833 and 23,813 genes for RA or SLE, respectively. By integrating eQTLs and summary statistics, we generated for each gene of both RA and SLE across all the tissues and adjusted it if λ > 1.05 (Supplementary Table 3). Then, the P-values were yielded. Afterwards, we further performed the HMP procedure to aggregate the P-values of each gene across all the tissues into a single P-value for RA or SLE. As mentioned before, these P-values from various tissues are in highly positive correlation with each other (Supplementary Figure 1), implying the failure if applying Fisher’s method which instead requires mutually independent P-values from different experiment tests (Fisher, 1934; Rice, 2010). The Manhattan plots for RA and SLE are shown in Figures 2A,B, with some associated genes highlighted.

FIGURE 2

(A) Manhattan plot for RA; (B) Manhattan plot for SLE. The genes in both plots with P-value < 2.1E-06 (the significance level adjusted by Bonferroni’s method) are highlighted; (C) Q-Q plot of RA conditional on the nominal P-value of SLE; (D) Q-Q plot of SLE conditional on the nominal P-value of RA. RA, rheumatoid arthritis; SLE, systemic lupus erythematosus.

According to the results of HMP, we conducted the cFDR analysis. In our analysis the Q-Q plot of RA conditional on the nominal P-value of SLE illustrates the existence of enrichment at different significance thresholds of SLE (Figure 2C). The presence of leftward shift suggests that the proportion of true associations for a given P-value of SLE would increase when the analysis is limited to include more significant genes. On the other hand, in terms of the Q-Q plot of SLE conditional on the nominal P-value of RA (Figure 2D), we observe a more pronounced separation in different curves, implying that there exists a stronger enrichment for SLE given RA than that for RA given SLE.

We further formally analyze the two diseases jointly using cFDR and show the results of association signals in Table 1 and Supplementary Tables 4, 5. Briefly, we identify 76 RA-associated genes (Supplementary Table 4) and 33 SLE-associated genes (cFDR < 0.05) (Supplementary Table 5). Among these genes, 59 RA-associated (e.g., CCBL2, SLC10A4, and PLEKHA1) and 19 SLE-associated genes (e.g., INPP5B, SKP1, and TMEM80) are not implied in the original GWASs of RA and SLE (Okada et al., 2014; Bentham et al., 2015), and are likely newly candidate associated genes for each disease. These findings also confirm that our multiple-tissue eQTL weighted integrative gene-based association analysis has higher power compared to the conventional single SNP analysis, as shown in many prior studies (Gusev et al., 2016; Xu et al., 2017; Guo and Wu, 2018; Wu and Pan, 2018; Xue et al., 2020; Zhang et al., 2020).

TABLE 1

Potential pleiotropic genes associated with RA and SLE identified by ccFDR.

Gene	chr	Position	cFDRRA	cFDRSLE	ccFDR
INPP5B	1	38,326,369-38,412,729	3.91E-16	5.64E-3	2.54E-2
OR5K2	3	98,216,448-98,217,496	1.42E-2	9.12E-4	4.63E-2
RP11-2C24.5	16	30,832,389-30,833,431	1.86E-2	2.34E-16	3.72E-2
CTD-3105H18.4	19	12,490,560-12,494,501	1.37E-2	2.79E-4	4.45E-2
MAGI3	1	113,933,371-114,228,545	1.81E-30	1.60E-2	1.60E-2
INPP1	2	191,208,196-191,236,391	7.07E-3	5.79E-5	4.24E-2
PDHB	3	58,413,357-58,419,584	2.51E-5	5.66E-4	6.60E-3
ITPR3	6	33,588,522-33,664,351	9.83E-27	4.60E-4	2.76E-3
CLIP2	7	73,703,803-73,820,273	4.88E-3	7.32E-5	3.09E-2
IRF5	7	128,577,666-128,590,089	5.45E-4	4.27E-11	2.72E-3
ICAM5	19	10,400,655-10,407,453	1.90E-3	1.57E-4	1.09E-2
TYK2	19	10,461,209-10,491,352	1.01E-6	2.57E-5	5.39E-4
CCDC116	22	21,987,007-21,991,616	1.06E-2	3.38E-16	3.18E-2
RP11-387H17.4	17	39,927,742-39,939,601	4.14E-5	1.06E-4	2.17E-3

Pleiotropic Genes Identified by ccFDR

Across all these RA- and SLE-associated genes, 14 genes (i.e., pleiotropic genes) are commonly related to RA and SLE (ccFDR < 0.05) (Table 1), four of which (i.e., INPP5B, OR5K2, RP11-2C24.5, and CTD-3105H18.4) are possibly new genes. Interestingly, the SNP effect sizes of 12 out of 14 genes (except PDHB and RP11-387H17.4) are highly positively correlated between RA and SLE (Supplementary Figure 2), with an average Pearson’s correlation r of 0.579. For example, the SNP effect sizes of IRF5 have a maximal r of 0.862 (95% confidence intervals [CIs] 0.826-0.891), followed by INPP1 (r = 0.825, 95%CIs 0.792-0.853) and TYK2 (r = 0.755, 95%CIs 0.701-0.799). This observation indicates that the genetic effects of these pleiotropic genes on the two diseases in general show a consistent direction. In addition, we also calculated the genetic risk score (GRS) for each pleiotropic gene (Supplementary Figure 3). The GRS is generated as the product of SNP effect sizes of RA (or SLE) and genotypes available from 503 European individuals of the 1000G (Consortium, 2015). We applied the GRS as an overall measurement of the genetic effect for a given pleiotropic gene on each disease. We observe that most of these pleiotropic genes (except INPP5B and RP11-2C24.5) have substantial different average GRS on RA and SLE. For example, five (i.e., OR5K2, INPP1, PDHB, ITPR3, and ICAM5) show a higher overall genetic effect on SLE, while the rest (i.e., CTD-3105H18.4, MAG13, CLIP2, INF5, TYK2, CCDC116, and RP11-387H17.4) show a higher overall genetic effect on RA. These GRS results provide an insight into the magnitude of the genetic effects of these pleiotropic genes on the two diseases.

The result of enrichment analysis shows that these pleiotropic genes have marked enrichment patterns mainly in type I interferon (IFN) signaling pathway, membrane, and inositol phosphate metabolic process (Table 2). The IFN signaling pathway plays a major role in activation of both innate and adaptive immune systems that are related to RA (Wright et al., 2015) and SLE (Bezalel et al., 2014). IRF5 and TYK2 are shown to be involved in this pathway, in line with prior studies (Acosta-Herrera et al., 2019) and supporting the validity of our results. Some pleiotropic genes (e.g., INPP5B, MAGI3, ITPR3, TYK2, and ICAM5) are enriched in the biological process of membrane fraction. MAGI3, as a novel membrane-associated guanylate kinase, is also implicated in the Wnt/β-catenin pathway, which induces promotion of regulatory T cell responses (Norén et al., 2017) as well as immune tolerance and plays a critical role in mucosal tolerance and suppression of chronic autoimmune pathologies (Suryawanshi et al., 2016). INPP5B, enriched in the pathway of inositol phosphate metabolic process and membrane fraction, was recently reported to be associated with membranes through an isoprenyl modification near the C-terminus and regulated calcium signaling by inactivating inositol phosphates (Nakatsu et al., 2015). Sustained calcium signaling responses are prevalent in the immunological synapse of T cells of SLE patients (Nicolaou et al., 2010). In addition, the mobilization of calcium signaling may modulate the functions of inflammatory and immunity genes in RA patients (de Seabra Rodrigues Dias et al., 2017).

TABLE 2

Enriched pathways for potential pleiotropic genes of RA and SLE.

Pathway	Genes	P
Phosphatidylinositol signaling system	INPP5B, INPP1, ITPR3	4.00E-3
Inositol phosphate metabolic process	INPP5B, INPP1	2.20E-2
Membrane	INPP5B, MAGI3, ITPR3, TYK2, ICAM5	2.40E-2
Type I interferon signaling pathway	IRF5, TYK2	3.00E-2
Inositol phosphate metabolism	INPP5B, INPP1	7.00E-2
Glucagon signaling pathway	ITPR3, PDHB	9.70E-2

As another example, PDHB, enriched in the glucagon signaling pathway, was discovered to be related to the promotion of aerobic glucose metabolism together with oxidative stress (Kanda et al., 2015). Up-regulation of glucose metabolism was demonstrated to be associated with upon activation of immune cells such as Fibroblast-like Synoviocytes (FLSs) (Garcia-Carbonell et al., 2016). The involvement of FLSs in regulating the pathogenesis of RA was highlighted in recent work (Meng and Qiu, 2020), which supports the validity of our findings.

Discussion

It has been widely observed that RA and SLE have common pathological and clinical features (Manoussakis et al., 2004; Icen et al., 2009; Orozco et al., 2011; Toro-Domínguez et al., 2014; Márquez et al., 2017; Acosta-Herrera et al., 2019), which are partly attributable to common genetic foundation between the two diseases. However, the genetic overlap underlying RA and SLE remains elusive and a large proportion of genes related to them are yet not discovered (Ramos et al., 2011). Large-scale GWASs undertaken on RA and SLE offer an unprecedented opportunity to tackle this question. The objective of our study was to gain insight into genetic mechanism linking between RA and SLE using advanced bioinformatics approaches. By leveraging publicly available GWAS summary statistics, we identified a significant positive genetic correlation between RA and SLE, indicating that genetic variants associated with the risk of SLE would be also related to the risk of RA, or vice versa (Icen et al., 2009; Toro-Domínguez et al., 2014).

By integrating eQTLs and combing association evidence across tissues, our study ultimately discovered 14 potential pleiotropy genes, and four of them (i.e., INPP5B, OR5K2, RP11-2C24.5 and CTD-3105H18.4) were not directly reported in previous literature, providing new insight into shared genetic basis between RA and SLE. Furthermore, we found that the SNP effect sizes of these genes were positively correlated and that these genes were implicated in multiple auto-immune relevant pathways such as inositol phosphate metabolic process, membrane and glucagon signaling pathway. As compared to previous cross-phenotype studies of autoimmune diseases (Orozco et al., 2011; Márquez et al., 2017), our study differs from them in multiple aspects and has several strengths. First, previous studies in general attempted to detect pleiotropy loci at the SNP level, while our work aimed to identify shared loci with gene as the functional unit. The power for detecting single SNP association is limited because genetic variants often have weak effect sizes (Hindorff et al., 2009; Visscher et al., 2017), making the detection of common associated SNPs difficult even with large samples. In contrast, due to the aggregation of multiple weak association signals and the reduced burden of multiple testing, gene-based analysis often has higher power than its counterpart of single SNP analysis. Therefore, our analysis is biologically more meaningful and statistically more powerful as widely demonstrated by gene-based association studies (Kwee et al., 2008; Wu et al., 2010, 2011; Schifano et al., 2012; Huang and Lin, 2013; Ionita-Laza et al., 2013; Wang et al., 2013; Lee et al., 2014; Zeng et al., 2014a,b, 2015a).

Second, as shown in previous studies, disease-associated SNPs were more likely to be eQTLs (Nicolae et al., 2010), implying that the functional roles of associated SNPs were regulated through gene expression; thus, the power improvement of gene-based association studies would be further achieved by integrating eQTLs into the test (Su et al., 2018; Wu and Pan, 2018; Xue et al., 2020). To do so, we systematically evaluate predicted gene expressions in RA and SLE through integrating eQTLs of relevant tissues from the GTEx project and GWAS summary statistics. Third, we aggregated association evidence across various tissues by applying the HMP procedure (Wilson, 2019a) which is robust against positive dependency among P-values and can produce a single well-calibrated P-value for evaluating the association. Fourth, prior studies were performed as across-disease meta-analysis for RA and SLE aimed to detected genetic loci that were associated with at least one disease rather than simultaneously related to both the diseases (Orozco et al., 2011; Márquez et al., 2017; Acosta-Herrera et al., 2019). Compared to these studies, besides the HMP combination strategy, our work also applied the widely used pleiotropy-informed method of cFDR to formally discover shared genes and provided a solid statistical foundation for our analysis (Zhernakova et al., 2009; Andreassen et al., 2013, 2014; Li et al., 2015; Ellinghaus et al., 2016; Smeland et al., 2020).

Finally, there are some limitations of our study needed to state. First, we cannot replicate these identified pleiotropic genes with external datasets or via in vivo and in vitro experiments. Second, the limited sample size of reference panel for calculating LD among SNPs and the well-known tissue-specific genetic effects of eQTLs may undermine the power of our integrative analysis. Third, although prior evidence described before indicates these newly pleiotropic genes may underlie certain aspects of the pathogenesis of RA and SLE, their biological mechanisms are still largely unclear; therefore, further studies are required to characterize functional roles of these genes on RA and SLE.

Conclusion

This study reveals common genetic components between RA and SLE and provides candidate associated loci for understanding of molecular mechanism underlying the comorbidity of the two diseases.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

PZ conceived the idea for the study. PZ and HL interpreted the results of the data analyses. PZ, JZ, and HL wrote the manuscript with help from other authors. All authors obtained the data and performed the data analyses.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The handling editor declared a past co-authorship with one of the authors PZ