Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci

RA is a common autoimmune disease that affects up to 1% of the general adult population worldwide¹. Approximately two-thirds of patients are seropositive for rheumatoid factor (RF) or anti-cyclic citrullinated peptide (anti-CCP) autoantibodies². Genetic studies in autoantibody positive RA among subjects of European ancestry have identified multiple major histocompatibility complex (MHC) region alleles and 25 confirmed RA risk alleles in 23 non-MHC loci^3–15. These alleles explain about 23% of the genetic burden of RA¹¹, indicating that additional alleles remain to be discovered.

To identify novel RA risk alleles, we conducted a GWAS meta-analysis of 5,539 autoantibody positive RA cases and 20,169 controls of European ancestry (Table 1). This study expands upon our previous GWAS meta-analysis of 3,393 cases and 12,462 controls¹¹ by (a) adding a new GWAS dataset of 483 RA cases recruited from the Boston area (Brigham Rheumatoid Arthritis Sequential Study, BRASS) and 1,449 shared controls, (b) adding 513 cases and 431 controls recruited from Sweden (Epidemiological Investigation of Rheumatoid Arthritis, EIRA), (c) incorporating a recently published GWAS of 2,418 cases and 4,504 controls recruited from North America (Canada and North American Rheumatoid Arthritis Consortium-III, NARAC-III)¹³, and (d) adding additional shared controls to the NARAC-III dataset. For each of six GWAS case-control collections, we removed SNPs and individuals that failed quality control, matched case-control samples using principal components analysis (PCA) to minimize bias due to stratification, and imputed¹⁶ genome-wide to infer genotypes at additional central European (CEU) HapMap SNPs. We used logistic regression to conduct a GWAS for 2.56 million SNPs in each collection, corrected for residual inflation using genomic control¹⁷, and combined results across collections by inverse variance-weighted meta-analysis (P_GWAS)¹⁸. We found little evidence of systematic bias, as indicated by the genomic control inflation factor λ_GC (1.04)¹⁷, quantile-quantile (Q–Q) plots, results for markers highly differentiated across Europe¹⁹, and comparisons with an alternative analysis using PCA eigenvectors as covariates (Supplementary Figures 1–4; Supplementary Table 1). See Methods and Supplementary Note for full details of the analysis.

Our GWAS meta-analysis found support for RA risk loci previously confirmed among individuals of European ancestry, consistent with power for our study design (Table 2, Supplementary Figure 5). Four of the 25 confirmed non-MHC risk alleles obtained P_GWAS<5×10⁻⁸ (PTPN22, CTLA4, TNFAIP3, and CD40 loci); the remaining 21 confirmed alleles obtained P_GWAS≤0.002. We found modest evidence for association at PADI4 (rs2240340, P_GWAS=0.01, odds ratio (OR)=1.06) and FCRL3 (rs3761959, P_GWAS=0.001, OR=1.08), but not CD244 (rs3753389, P_GWAS=0.26, OR=1.03), loci implicated in RA patients of Asian ancestry^3,20,21 (Supplementary Table 1 online).

We used three criteria to select 34 independent SNPs for replication after removing confirmed RA risk alleles. First, we selected all 11 SNPs with P_GWAS<10⁻⁶. Second, we selected 9 SNPs (of 116 tested) with P_GWAS<0.0001 and identified as functionally related to known RA risk loci by GRAIL²², a bioinformatic analysis that identifies connections among genes in published abstracts (Supplementary Figure 6 online). Third, we selected 14 SNPs (of 104 tested) with P_GWAS<0.01 and previously known to be associated with other autoimmune diseases^23–31, as we found evidence for enrichment of autoimmune SNPs in our RA GWAS (Supplementary Figure 7 online). See Supplementary Note online for additional details about our selection of SNPs for replication.

These 34 SNPs were genotyped in an independent set of 6,768 autoantibody positive RA cases and 8,806 matched controls of European ancestry (Table 1). As in our GWAS, all cases were seropositive for disease-specific autoantibodies (anti-CCP or RF). For each SNP, we tested for replication of the GWAS meta-analysis association by calculating a one-tailed P-value (P_replication) for the same allele conferring risk in both analyses. Additionally, we conducted overall association tests across all 41,282 samples (GWAS meta-analysis plus replication, P_overall), and considered P_overall<5×10⁻⁸ to be a reproducible level of significance for RA risk association.

Of the 34 SNPs tested, 12 obtained a Bonferroni-corrected P_replication<0.0015 (=0.05/34), and an additional 9 obtained P_replication<0.05. Ten of these achieved genome-wide significance in the combined analysis, P_overall<5×10⁻⁸, indicating that these are validated RA risk alleles (Table 3). The other 11 SNPs had P_replication<0.05, but P_overall>5×10⁻⁸, indicating highly suggestive but not genome-wide significant association (Table 4). Results for all 34 SNPs in replication and combined analyses are shown in Supplementary Table 2 online.

Three of the 10 newly validated RA risk loci have not been implicated in any genetic studies of RA or other autoimmune disease. The SNPs are located at chromosome 2p14 (rs934734, P_overall=5.3×10⁻¹⁰), 5q11 (rs6859219, P_overall=9.6×10⁻¹²) and 5q21 (rs26232, P_overall=4.1×10⁻⁸). Replication sample odds ratios (for minor alleles) were 1.13, 0.85 and 0.93, respectively (Table 3).

Although no single gene can be declared causal, we labeled RA risk loci with the names of the most compelling candidate gene(s) from each region of linkage disequilibrium (LD) based upon GRAIL analysis and/or knowledge of RA pathogenesis^2,22. At 2p14, the most significant SNP (rs934734) is located within intron 1 of the sprouty-related, EVH1 domain containing 2 (SPRED2) gene (Figure 1A), which has been shown to regulate CD45+ hematopoietic cells via the Ras/MAP kinase pathway³². At 5q11, the most significant SNP (rs6859219) is located in ANKRD55, an ankyrin repeat domain-containing gene of unknown function (Figure 1B). A more compelling immunological candidate, interleukin 6 signal transducer (IL6ST), lies ~150 kilobases (kb) proximal but outside the region of LD with associated SNPs. The IL6ST protein product, gp130, functions as a part of the receptor complex for the inflammatory cytokine IL6³³. At 5q21, there is no obvious biological candidate gene; SNP rs26232 lies within the intron of predicted gene C5orf30 (Figure 1C).

Our study provides the first convincing evidence that 4 loci implicated in other autoimmune diseases are also associated with risk of RA. Of these, three of the 4 SNPs were selected for replication based on obtaining P_GWAS<10⁻⁶, regardless of their previously reported associations with autoimmune disease. These 3 SNPs are located at chromosome 3p14 (rs13315591, near PXK, P_overall=4.6×10⁻⁸), 4p15 (rs874040, near RBPJ, P_overall=1.0×10⁻¹⁶) and 6q27 (rs3093023, in CCR6, P_overall=1.5×10⁻¹¹). The 3p14 SNP lies 187 kb proximal to a SNP associated with systemic lupus erythematosus (SLE)²⁵; the SLE SNP, intronic in the PXK gene, is only weakly associated with RA risk in our study (rs6445975, P_GWAS = 0.03, r² = 0.15 and D’ = 0.75 with rs13315591). The 4p15/RBPJ SNP is in complete LD (r² = 1) with a SNP associated with risk of type 1 diabetes (T1D, rs10517086)²³, and the same allele confers risk in both diseases. The 6q27/CCR6 SNP rs3093023 is in LD with a SNP associated with Crohn’s disease²⁴, which is only weakly associated with RA risk in our study (Crohn’s SNP rs2301436, P_GWAS = 0.045, r² = 0.48 and D’ = 0.80).

The IRF5 SNP rs10488631 (P_GWAS =2.8×10⁻⁶), chosen because of its association in SLE^34,35, was convincingly associated with autoantibody positive RA in our study (P_replication=1.2×10⁻⁶, P_overall=4.2×10⁻¹¹). Graham et al.³⁶ proposed that the complex SLE risk associations in IRF5 are explained by three independent groups of SNPs, each consisting of SNPs in tight LD with each other. In our dataset, these groups are represented by rs10488631 (group 1), rs729302 (group 2) and rs4728142 (group 3). In addition to the group 1 SNP, rs10488631, we found evidence for association with group 3 SNP rs4728142 (P_GWAS=7.2×10⁻⁶), which is in LD with a variant that alters IRF5 polyadenylation and expression³⁶. However, we found no evidence for association for group 2 (rs729302, P_GWAS=0.15). Conditional analyses indicated that the group 1 and group 3 effects are independent (Supplementary Table 3 online). We note that group 3 SNP rs3807306 has been suggested to be associated with autoantibody negative RA³⁷, and our results suggest that this SNP is associated with autoantibody positive RA (P_GWAS=4.5×10⁻⁶).

Our GWAS meta-analysis identified new alleles exhibiting independent effects from previously implicated RA risk alleles at two other loci, CCL21¹¹ (rs951005, P_overall=3.9×10⁻¹⁰) and IL2RA¹⁵ (rs706778, P_overall=3.3×10⁻¹¹). Our new CCL21 SNP has little LD with the previous RA risk-associated SNP rs2812378¹¹ (P_GWAS=1.6×10⁻⁷; r²=0.06 and D’=0.83). Conditional analyses indicated that the associations are partly independent (Supplementary Table 3 online); the new SNP rs951005 remains significant (P = 1.7×10⁻⁵) conditional on the previously known SNP rs2812378; conditional on rs951005, rs2812378 remains nominally significant (P = 0.01).

SNPs at IL2RA are known to be associated with multiple autoimmune diseases including RA¹⁵, T1D²³ and multiple sclerosis (MS)^38,39. The IL2RA SNP identified in our GWAS (rs706778, P_overall=3.3×10⁻¹¹) is in partial LD with a SNP previously implicated in RA, T1D and MS (rs2104286, r²=0.25 and D’=1)¹⁵ and another SNP implicated in T1D and MS (rs11594656; r²=0.24 and D’=1)²³. In fact, T1D risk is conferred by a haplotype tagged by rs2104286(T)-rs11594656(T)³⁸. In our GWAS dataset, rs706778 correlated strongly with the T1D risk haplotype (r²=0.83), and conditional and haplotype analyses demonstrated that risk of RA is better explained by this haplotype [rs2104286(T)-rs11594656(T)] than by rs706778 (Supplementary Tables 3 and 4 online).

A SNP at the AFF3 locus has been previously implicated in RA¹⁴, with equivocal evidence for association with T1D^23,40. Our study provides definitive evidence that this locus is associated with risk of autoantibody positive RA (rs11676922, P_overall=1.0×10⁻¹⁴, OR=1.14).

Eighteen of the 34 SNPs that we tested in replication (based on the criteria above) are in LD with validated autoimmune risk alleles (Figure 2). Of these, 5 demonstrated genome-wide (P_overall<5×10⁻⁸) and 7 showed suggestive (P_replication<0.05) evidence of association in our study; only 6 showed no evidence of association in replication. While not meant as a complete comparison of all known RA and other autoimmune disease risk alleles, these results underscore the emerging overlap in genetic bases of RA and other autoimmune diseases^41,42. Although additional replication in large sample collections is required to confirm the suggestive associations, many likely represent true RA risk alleles.

Several of the SNP associations further implicate T-cells in RA pathogenesis. RBPJ (recombination site binding protein J; also known as CSL) is a transcription factor within the Notch signaling pathway. RBPJ-deficient mice have no T-cell development, while early B-cell development is maintained⁴³. CCR6 is a cell surface protein that distinguishes Th17 cells from other CD4+ helper T-cells⁴⁴. Synoviocytes from arthritic joints of mice and RA patients produce CCL20, a CCR6 ligand. Intriguingly, anti-CCR6 monoclonal antibodies substantially inhibit mouse arthritis, suggesting that CCR6 could be a therapeutic target in RA⁴⁴. CD247 (P_replication=0.006, P_overall=1.6×10⁻⁶; Table 3) encodes the T-cell receptor zeta chain, a subunit of the T-cell receptor-CD3 complex that when mutated causes an inflammatory arthritis in mouse^45,46. The zeta chain plays an important role in coupling antigen recognition to several intracellular signal-transduction pathways. These results build upon knowledge of T-cell activation-differentiation in RA pathogenesis evidenced by genetic associations at HLA-DRB1 (presents antigens to T-cells), PTPN22 (alters T-cell thresholds), CTLA4 (co-stimulator for T-cell activation), IL2RA/CD25 (mediates IL2-dependent T-cell responses), and STAT4 (transcription factor in Th1 cell differentiation), among others.

In conclusion, we find convincing evidence for association with risk of autoantibody positive RA at 10 loci, seven of which represent novel risk loci for RA. We estimate that the now >30 validated non-MHC RA risk alleles explain 3.9% of the total disease variance, with 0.67% due to the new RA risk alleles reported here. It is clear that additional risk alleles remain to be identified, as genetic discoveries explain only 16% of disease variance (including an estimated 12% for the MHC⁴⁷), whereas more than half the risk of autoantibody positive RA is thought to be genetic^47,48. The number of SNPs with suggestive evidence of association in our study (Table 4) further indicates that many more common risk alleles with modest effect sizes remain to be discovered. In addition to common variants, the roles of rare variants, copy number variants and epigenetic modifications will need to be explored with newer genomic technologies.