Genetic architecture study of rheumatoid arthritis and juvenile idiopathic arthritis

Background Rheumatoid arthritis and juvenile idiopathic arthritis are two types of autoimmune diseases with inflammation at the joints, occurring to adults and children respectively. There are phenotypic overlaps between these two types of diseases, despite the age difference in patient groups. Methods To systematically compare the genetic architecture of them, we conducted analyses at gene and pathway levels and constructed protein-protein-interaction network based on summary statistics of genome-wide association studies of these two diseases. We examined their difference and similarity at each level. Results We observed extensive overlap in significant SNPs and genes at the human leukocyte antigen region. In addition, several SNPs in other regions of the human genome were also significantly associated with both diseases. We found significantly associated genes enriched in 32 pathways shared by both diseases. Excluding genes in the human leukocyte antigen region, significant enrichment is present for pathways like interleukin-27 pathway and NO2-dependent interleukin-12 pathway in natural killer cells. Discussion The identification of commonly associated genes and pathways may help in finding population at risk for both diseases, as well as shed light on repositioning and designing drugs for both diseases.


INTRODUCTION
Rheumatoid arthritis (RA) is a symmetric polyarticular arthritis that primarily affects the small diarthrodial joints of the hands and feet, while juvenile idiopathic arthritis (JIA) is caused by unknown etiology and persists at least 6 weeks in children under the age of 16, which does not contain other known conditions (Firestein, 2003;Prakken, Albani & Martini, 2011). The prevalence rate of RA varies from 0.41 to 0.54% from 2004 to 2014 based on US administrative health insurance claims databases (Hunter et al., 2017), which is observably greater than the prevalence rate of JIA ranging from 0.0038 to 0.40% according to a systematic literature review including 29 articles (Thierry et al., 2014).
Phenotypically, RA and JIA are similar in some aspects. They show some common symptoms and physical signs such as joint pain and swelling, limited joint mobility and deformity, morning stiffness, elevated rheumatoid factor, fever, etc. Some of the subtypes of JIA, such as polyarticular JIA, are particularly similar to RA. However, with distinct clinical and pathological features of these two diseases being noted, they have been defined as separate diseases by International League of Associations for Rheumatology (Petty et al., 2004). In particular, JIA is heterogeneous with variable clinical presentation and outcome. It has been classified into seven subtypes, including oligoarticular JIA (persistent/ extended), polyarticular JIA with negative rheumatoid factor (RF), polyarticular JIA with positive RF, psoriatic JIA, enthesitis related arthritis, systemic JIA and undifferentiated JIA (Nigrovic, Raychaudhuri & Thompson, 2018). RA is more homogeneous but with a poorer outcome.
It has long been recognized that both RA and JIA are related to autoimmune and inflammatory disorders (Ravelli & Martini, 2007;Scott, Wolfe & Huizinga, 2010). Population-based heritability estimates for RA and JIA are both 60% approximately (Macgregor et al., 2000;Prahalad, 2006). Single-nucleotide polymorphism (SNP)-based heritability for RA has been reported to be around 50% (Speed & Balding, 2014;Speed et al., 2012), slightly lower than that of JIA estimated to be 73% (Li et al., 2015b). Certain alleles in the HLA region are strong genetic predisposition factors for RA and JIA. It has been reported that for both RA and JIA, the odds ratio of HLA region is about 2.8, while that of most non-HLA loci is only 1.1 to 1.4. (Nigrovic, Raychaudhuri & Thompson, 2018;Hersh & Prahalad, 2015;Viatte, Plant & Raychaudhuri, 2013) The genetic predisposition of JIA is attributable to HLA class II molecules (HLA-DRB1, HLA-DPB1), HLA class I molecules and non-HLA genes. The clinical presentation of RF-positive JIA resembles that of RA, and they share the HLA-DRB1 epitope (De Silvestri et al., 2017;Hinks et al., 2018). The HLA-DRB1*04 confers a protective role in JIA before the age of 6, while it renders an increased risk of RA (Nigrovic, Raychaudhuri & Thompson, 2018). The immunopathogenesis of RA has become clear in recent years, but the pathogenesis of JIA remains unknown (Firestein & McInnes, 2017;Mellins, Macaubas & Grom, 2011).
With the rapid development of genomic technology, a large number of genetic variants associated with RA or JIA have been identified. To date, genome-wide association studies (GWASs) have identified a large number of variants associated with RA and JIA respectively. A total of 789 RA-associated variants from 52 studies and 129 JIA-associated variants from 11 studies have been reported in GWAS Catalog (association testing P-value <1 × 10 −5 ) (Buniello et al., 2019), including 101 loci associated with RA and around 30 loci associated with JIA at genome-wide significant level. We aimed to compare the genetic architecture of RA and JIA at multiple levels systematically.
In this study, we conducted gene, pathway and network analyses of RA and JIA using robust and computational efficient methods based on their summary GWAS statistics. We compared genetic difference and similarity between RA and JIA, identified their shared genetic signature. Significant overlap in genes and pathways were observed between these two diseases.

Data collection
RA genetic loci information came from GWAS summary statistics of a trans-ethnic study (Okada et al., 2014) including a total of 29,880 RA cases and 73,758 controls of European and Asian ancestries. Summary statistics was downloaded from GWAS catalog (https://www.ebi.ac.uk/gwas/) (Buniello et al., 2019). All RA patients met the RA diagnostic criteria established by the American College of Rheumatology in 1987(Arnett et al., 1988, or were confirmed by a professional rheumatologist (Okada et al., 2014).
JIA genetic loci information came from two resources. First, summary statistics of our previous GWAS on JIA (Finkel et al., 2016) was included in the current study. Our JIA GWAS is composed of discovery and replication cohorts with 1166 JIA cases and 9500 unrelated controls of European ancestry totally. Summary statistics of meta-analysis on the discovery and replication cohorts were used in our current study. Second, JIA variants revealed in published GWASs (Behrens et al., 2008;Cobb et al., 2014;Finkel et al., 2016;Hinks et al., 2009;Hinks et al., 2013;Li et al., 2015a;Ombrello et al., 2017;Thompson et al., 2012) were extracted from GWAS catalog (Buniello et al., 2019).

Gene-based association analysis
A gene-based association analysis for RA and JIA was performed using fastBAT method (Bakshi et al., 2016) implemented in GCTA v1.91.7 (Yang et al., 2011 respectively, based on GWAS summary statistics of RA or JIA and linkage disequilibrium (LD) information from EUR population in the 1000 Genomes Project (The Genomes Project Consortium et al., 2015). Each gene region was defined as its transcript region and 50kb upstream/downstream, and the threshold for LD pruning was set to r 2 -value >0.9, following the default setting of fastBAT. The gene list of human genome used by fastBAT method contains 24765 annotated genes (Bakshi et al., 2016), thus the genome-wide significant threshold for gene based tests was set at 0.05/24765 = 2 × 10 −6 . JIA SNPs in GWAS catalog was also mapped to genes according to its report (Buniello et al., 2019).

SNP-level comparison
A total of 26,285 SNPs (Table S1) in RA study and 105 SNPs (Tables S2, S3) in JIA study reached genome-wide significance threshold P-value <5 × 10 −8 , and these two diseases shared 47 significant SNPs. Among these SNPs, 37 were located in the human leukocyte antigen (HLA) region on chromosome 6. The rest 10 SNPs were located in or close to 9 genes (Table 1). Interestingly, 8 SNPs located in the HLA region showed opposite direction of effects, which meant risk allele of JIA could be protective allele for RA and vice versa.

Gene-based comparison
To increase statistical power and to consider the combined effects of SNPs in genes, we conducted gene-disease association analyses, based on SNP-level summary statistics and taking into account of LD between SNPs. Several methods have been developed for computing gene-level associations based on SNP-level summary statistics, such as the commonly used PLINK (Purcell et al., 2007) set-baesd test and software VEGAS (Versatile Gene-based Association Study) (Liu et al., 2010), which are permutation and simulation-based approaches respectively. Both methods rely on resampling which is computationally intensive. Here, we adopted the fastBAT method which was a robust set-based association test computing the P-value of a gene with a number of SNPs from an approximated distribution (Bakshi et al., 2016). 431 genes located at 50 loci reached genome-wide significance in the RA dataset, including 17 known loci (Acosta-Herrera et al., 2019;Buniello et al., 2019;Eyre et al., 2012;Plenge et al., 2005;Raychaudhuri et al., 2009;Zhu et al., 2016) and 33 novel loci which should be examined in future replication studies (Table S4).
However only genes in the HLA region showed genome-wide significant association with JIA, which was likely due to the limited power of our previous GWAS (Table S5). A total of 75 significant genes or regions in the HLA were shared by JIA and RA (Table  S6). Then we checked whether significant genes in RA contained additional genome-wide significant SNPs in JIA reported in GWAS catalog. Not surprisingly, one RA significant gene in the HLA region and 8 genes outside the HLA region containing genome-wide significant SNPs for JIA (Table 2) were observed. Because the fastBAT method conducted LD-pruning before combining SNP statistics, the top SNP showed in Table 2 may not be the one with the best P-value in original GWAS. Table 1 Genome-wide significant SNPs shared by RA and JIA (P-value < 5 × 10 −8 ). The raw data of genome-wide significant SNPs of RA are presented in Table S1; and the raw data of genome-wide significant SNPs of JIA are shown in Tables S2 and S3.   Notes. Chr, chromosome; Start-End, start and end boundaries of the gene region on human genome build UCSC hg19 (NCBI GRCh37); RA, rheumatoid arthritis; JIA, juvenile idiopathic arthritis; Pval, gene-level P-value based on fastBAT method; TopSNP, the top associated GWAS SNP; TopSNP_Pval, smallest single-SNP GWAS P-value in the gene region.

Pathway-level comparison
GWAS pathway analysis consider either competitive null hypothesis or self-contained null hypothesis. Many methods for GWAS pathway analysis have been developed, but they are still subjected to the issues of low power and being influenced by some free parameters. The recently developed GSA-SNP2 package (Yoon et al., 2018) uses the random set model to compute pathway enrichment with decent type I error control by integrating the gene scores adjusted by the number of SNPs mapped to each gene and removing high inter-gene correlated adjacent genes in each pathway. It does not require any key free parameters concurrently. We applied this method to our analyses. RA or JIA associated genes were enriched in numerous canonical pathways at a threshold of Q-value <0.05. A total of 32 enriched pathways were shared by RA and JIA, which mostly were immune-related pathways, such as allograft rejection, type 1 diabetes mellitus, graft versus host disease,

Notes.
Pathway, abbreviation for each enriched pathway; Database, database from which the pathways were extracted; Size, total number of genes in each pathway; RA, rheumatoid arthritis; JIA, juvenile idiopathic arthritis; Count, the number of RA/JIA-significant genes falling into each pathway; Pval, P-value of each pathway; Qval, Q-value of each pathway based on the trend curve adjusted gene scores.
antigen processing and presentation, autoimmune thyroid disease, asthma, etc. (Table  S7). Most of these significant pathways were driven by genes in the HLA region. In order to explore the role of loci outside the HLA region for these two diseases, we performed pathway enrichment analysis again after removing loci in the HLA region based on their genomic coordinates. The HLA region was defined as chr6:28,477,797-33,448,354 (GRCh37/hg19). Pathways such as interleukin(IL)-27 pathway and NO2-dependent IL-12 pathway in natural killer (NK) cells were significantly enriched even after the HLA region loci were removed (Table 3). Global networks were visualized at a threshold of gene-score <0.005 (Figs. 1&2). We observed the common hub role of several genes such as TYK2. The networks before removing the HLA region were shown in Figs. S1 and S2.

DISCUSSION
Despite the phenotypic similarity between JIA and RA, systematic comparison of genetic similarity and distinction between these two types of diseases are lacking. Large scale GWASs of RA and JIA respectively render us ability to conduct such comparison and to identify potential common mechanism in disease pathogenesis, which may help repositioning and designing treatment strategies.
To systematically compare the genetic architecture of the two diseases, we performed gene-level, pathway-level analyses and conducted comparison at each level. Not only did we observe a large amount of overlaps in the HLA region as expected, but we also observed several SNPs and genes which significantly associated with both diseases outside the HLA region. Among them, the risk alleles of several SNPs were different between the two diseases, which meant that a certain allele may play a risk role in one disease but a protective role in the other. These SNPs might be related to the differences in pathogenesis and phenotype between JIA and RA. As we did not perform genome-wide imputation analysis due to unavailability of individual-level data, the number of genome-wide significant SNPs shared by these two diseases was actually underestimated.
Due to the limited sample size of our JIA data, we could not perform analysis for each subtype of JIA with enough statistical power. However, the heterogeneity of JIA and the genetic basis of its subtypes are worth noting. Some HLA alleles show different directions of effects on different subtypes of JIA and RA. For instance, HLA-DRB1*8, HLA-DRB1*11 and HLA-DRB1*13 are risk alleles of seronegative JIA, but do not exhibit association with seropositive polyarticular JIA and seronegative RA, and these HLA alleles render protective effect for seropositive RA. In particular, DRB1*11 is also a risk allele of systemic JIA, while the other two alleles are not associated with this JIA subtype (Nigrovic, Raychaudhuri & Thompson, 2018). As for alleles outside the HLA region, certain SNPs in genes PTPN22 and STAT1/STAT4 do not show association with systemic JIA, but confer risk for most other subtypes of JIA and RA (Nigrovic, Raychaudhuri & Thompson, 2018 Figure 2 The global network of JIA after the HLA region being removed (Q-value < 0.05, gene-score < 0.005). The PPI network was constructed among proteins encoded by the significant JIA-associated genes excluding those in the HLA region. The nodes in the figure represent the proteins and the connections between nodes indicate protein-protein interactions. The size of each node suggests the degrees of the connection between the node and the others. Full-size DOI: 10.7717/peerj.8234/ fig-2 on Immunochip (Hinks et al., 2018;Onuora, 2018). Further analysis of the genetic nature of different subtypes of JIA and RA would be helpful to optimize the classification of the two diseases, and may lead to more effective treatment and better prognosis.
We observed significant enrichment of NO2-dependent IL12 pathway and IL27 pathway for both RA and JIA. Macrophages release IL-12 which plays an important role in activation of NK cells and induces cytotoxicity with nitric oxide (Liu et al., 2005). NK cells are regarded as a bridge between innate and adaptive immunity, serving as a key regulator in the pathogenesis and development of autoimmune diseases (Gianchecchi, Delfino & Fierabracci, 2018). It has been reported that high percentages of NK cells and their activity were found in synovial fluid of active RA patients at advanced stage (Yamin et al., 2019), and dysfunction of NK cells was also observed in patients with systemic-onset JIA and its complication (Grom et al., 2003). NO2-dependent IL12 pathway plays a unique role in the activation of NK cells by macrophage. The enrichment of this pathway in our analyses implies the potential role of abnormal IL-12-mediated activation of NK cell in the pathogenesis of RA and JIA. IL-12 has long been considered as a therapeutic target of arthritis and other autoimmune and inflammatory disorders (Hasko & Szabo, 1999;Siebert et al., 2015). As a member of the IL-12 family, IL-27 induces T cell differentiation and causes immunosuppressive effects by inhibiting the development of Th17 cells (Yoshida & Miyazaki, 2008). Previous studies have suggested that IL-27 is another key modulator of autoimmunity and elevation of IL-27 signaling may be inhibitory to some autoimmune diseases, such as multiple sclerosis or uveitis (Amadi-Obi et al., 2007). Our results suggest that such therapeutic approach may be also applied to the management of RA and JIA.

CONCLUSION
Our study identified genetic similarities and differences between RA and JIA at multiple levels. We observed a number of genes being associated with both diseases especially in the HLA region, and distinct genetic loci were found as well. Such systematic comparison and further functional characterization of these genetic loci and signaling pathways may lead to the identification of common drug targets for both diseases or drug repositioning, and may also contribute to the precision treatment of each disease.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by National Natural Science Foundation of China (No.81771769), Tianjin Natural Science Foundation (No.18JCYBJC42700), Startup Funding from Tianjin Medical University. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: National Natural Science Foundation of China: No.81771769. Tianjin Natural Science Foundation: No.18JCYBJC42700. Tianjin Medical University.