Using the Immunochip custom SNP array, which was designed for dense genotyping of 186 loci identified through genome-wide association studies (GWAS), we analyzed 11,475 individuals with rheumatoid arthritis (cases) of European ancestry and 15,870 controls for 129,464 markers. We combined these data in a meta-analysis with GWAS data from additional independent cases (n = 2,363) and controls (n = 17,872). We identified 14 new susceptibility loci, 9 of which were associated with rheumatoid arthritis overall and five of which were specifically associated with disease that was positive for anticitrullinated peptide antibodies, bringing the number of confirmed rheumatoid arthritis risk loci in individuals of European ancestry to 46. We refined the peak of association to a single gene for 19 loci, identified secondary independent effects at 6 loci and identified association to low-frequency variants at 4 loci. Bioinformatic analyses generated strong hypotheses for the causal SNP at seven loci. This study illustrates the advantages of dense SNP mapping analysis to inform subsequent functional investigations.

Rheumatoid arthritis is a common, complex disease affecting up to 1% of the adult population. It is an archetypal autoimmune disease, typified by the presence of serum autoantibodies, including antibodies directed against the Fc portion of immunoglobulins (rheumatoid factor) and against citrullinated peptides (anti-citrillunated peptide antibodies (ACPA)). Genetic studies of rheumatoid arthritis, including recent application of genome wide association studies (GWAS), have identified 32 risk loci among individuals of European ancestry, including HLA-DRB1, PTPN22, and other loci with shared autoimmune associations1, 2.

The Immunochip Consortium was formed to design a custom Illumina Infinium array that leveraged the remarkable genetic overlap of susceptibility loci identified across a range of autoimmune diseases. The custom array allows investigators to perform gene-finding and fine-mapping experiments in a co-ordinated manner. Full details have been described previously3. Briefly, the array consisted of all known single nucleotide polymorphisms (SNPs) from the 1000 Genomes Project as well as private resequencing efforts for 186 loci, known to be involved in 12 autoimmune diseases. For these loci there is the unique opportunity to fine map autoimmune disease associations. Additional SNPs were included as part of a deep replication effort. This not only provided the opportunity to identify novel rheumatoid arthritis associations with other autoimmune disease loci or with variants with suggestive statistical evidence for association from a previous meta-analysis of but also to refine the GWAS signal and reduce the number of potential causal variants in the 31 non-HLA confirmed loci.

We tested 129,464 polymorphic markers passing quality control, with a minor allele frequency >1%, in 11,475 cases (7,222 ACPA positive, 3,297 ACPA negative and 957 unassigned) and 15,870 controls (Table 1 and Supplementary Tables 1 and 2). We performed analysis on the total rheumatoid arthritis dataset, and also in subsets stratified by ACPA status (Supplementary Table 3). We also had access to GWAS data for an additional 2,363 ACPA positive cases and 17,872 controls, independent of the current study (Table 1). We observed strong evidence of association for the previously identified susceptibility loci (Table 2 and Supplementary Tables 3 and 4).

We identified fourteen novel rheumatoid arthritis loci for populations of European ancestry (TYK2, IRAK1, TLE3, RASGRP1, PADI4, IL6R, IRF8, ARID5B, IKZF3, RUNX1, POU3F1, RCAN1, CD5, GATA3) at genome-wide levels of significance (p<5×10−8)(Figure 1): 7 with Immunochip data alone (Table 2) and a further 7 when Immunochip data was combined with the GWAS meta analysis data (Table 2). These loci add 4% to the estimate of heritability explained by confirmed loci, bringing the total to 51%, of which HLA explains 36%. When we removed all known loci from the Immunochip data, we still observed evidence of an excessive number of nominally associated alleles, consistent with the possibility that there are many additional undiscovered alleles 4(Supplementary Figure 1). Interestingly, if a study-wide significance threshold of 9.0×10−7 is applied (calculated based on the number of effective independent tests when accounting for linkage disequilibrium (LD)), significant association is also observed at two additional loci; ELMO1 (rs75351767 pall=2.94×10−7) and BACH2 (rs72928038 Pall = 8.23×10−7) (Supplementary Table 3). A further 8 loci are implicated at suggestive levels of significance (p<1 ×10−5) in either the full or ACPA positive sub-group analysis including PTPN2 (rs62097857 Pall=4.4×10−6); TNIP1 (rs6579837 Ppos=1.7×10−6) and TNFSF4 (rs61828284 Ppos=5.4×10−6) (Supplementary Table 3).

Previously, we have fine-mapped MHC associations observed in GWAS data of partially overlapping samples by applying imputation of HLA classical alleles and amino acids5. The Immunochip platform includes denser SNP coverage within the MHC region which facilitates more accurate imputation. In a preliminary analysis applying the same imputation and fine-mapping approach to ACPA positive cases and controls typed on Immunochip, we observed the same associations that we reported previously. The most significant polymorphic nucleotide was again rs17878703, mapping to position 11 of the HLA-DRB1 peptide sequence (p<10−677). Testing individual amino acid positions within HLA-DRB1 revealed the strongest association at position 11 (p<10−745); conditioning on the position 11 effect we observed association at position 71 (p=6×10−60); finally conditioning on effects at both positions 11 and 71 we observed significant association at position 74 (p=7×10−19). Adjusting for all HLA-DRB1 alleles to identify independent effects outside this gene we observed significant associations at HLA-B corresponding to the presence of aspartate at position 9 in the peptide sequence (p=1×10−17). Adjusting for all HLA-DRB1 alleles and Asp-9 in HLA-B, we observed associations at HLA-DPB1 corresponding to the presence of phenylalanine at position 9 in the peptide sequence (p=1×10−17).

While it has been demonstrated that ACPA positive and ACPA negative disease has a different allelic association at the MHC and at PTPN226, previous studies have not been powered to address this issue definitively in additional non MHC loci. Here we analysed 3,297 ACPA negative cases and identify association at genome wide significance to ANKRD55 (rs71624119 p=5.2 ×10−12, OR=0.78) in addition to HLA (rs4143332 p=2.9×10−15, OR=1.37) (Supplementary Table 3). Strikingly, ANKRD55 has a similar effect as in ACPA positive disease. Comparing association in ACPA positive and negative subgroups we see that for the 45 non-HLA loci, around half show a significantly larger effect size in ACPA positive disease (comparison of OR p<0.05), 5 of these loci having a markedly stronger association with this form of disease (PTPN22, CCR6, CD40, RASGRP1 and TAGAP). Eleven loci show no statistical difference in association to either form of rheumatoid arthritis (Supplementary Table 5).This preliminary analysis indicates that differences in the serological subtype of disease may well be reflected in a difference in genetic pre-disposition potentially providing a basis for stratified medicine.

The majority of the 14 new loci associated with rheumatoid arthritis susceptibility, along with previously confirmed loci, were found to contain proteins strongly linked to immune function using GRAIL analysis (Supplementary Table 6 and Supplementary Figure 2), for example, CD5, IRF8 and TYK2. We also report novel association with IRAK1, previously associated with systemic lupus erythematosus (SLE)7. This is the first X chromosome locus association to rheumatoid arthritis, and is of relevance given the female predominance of both diseases (9:1 and 3:1 ratio of females: males in SLE and rheumatoid arthritis, respectively). Interestingly, this locus has been shown to occasionally escape X inactivation in female cells8. Three of the novel loci confirmed here for the first time in samples of European ancestry have previously been associated in either samples of East Asian ancestry (PADI4, ARID5B) or when using a multiethnic approach (IKZF3)9-11. The SNPs associated in this study are moderately correlated with those identified in samples of East Asian origin, PADI4 SNPs rs2240336 and rs766449 r2=0.25, D′=1; ARID5B SNPs rs12764378 and rs10821944 r2=0.52, D′=0.86. PADI genes are involved in the citrullination of peptides and as such are strong candidates for involvement in disease, given the presence of ACPA auto-antibodies. Although the association at PADI4 (rs2240336) is greater in ACPA positive disease (OR=0.88 P=6.49×10−9) compared to ACPA negative cases (OR=0.93, P=0.01) our formal test comparing OR did not show a statistically significant difference (P=0.14).

We applied conditional logistic regression to test for secondary effects within each locus. In 6 non-HLA loci (13%) (TNFAIP3, CD28, REL, STAT4, TYK2, RASGRP1) we observed additional independent association signals (Supplementary Figure 3). In total we observed 51 independent risk alleles in 45 non-HLA rheumatoid arthritis loci. To test the possibility that the two risk alleles tag an untyped SNP, we carried out haplotype analysis of the six loci but found no evidence for haplotype specific effects at any locus (Supplementary Table 7). At only four loci, REL, CD28, TYK2 and TNFAIP3, did we observe associations with low frequency variants (MAF<0.05) (Supplementary Table 8).

Out of the 46 rheumatoid arthritis loci, 39 were densely genotyped by Immunochip. For 12 loci we observed that the most strongly associated SNP was not tightly linked to the previously reported leading SNP at that locus, shifting the association signal (Supplementary Table 9). For the 39 confirmed non-HLA rheumatoid arthritis loci on Immunochip, dense mapping refines the association to a single gene for 19 loci (Supplementary Table 10).

Our analysis also identified 7 non-synonymous SNPs within exonic regions (Table 3), as well as a number showing strong regulatory potential (Supplementary Table 11), that are highly correlated (r2>0.9) with the lead SNP and which are strong candidates for the aetiological variant. The most associated SNP at the IL6R locus (rs8192284) is non-synonymous, shows high correlation with circulating IL6R levels and as well as being associated with a decrease risk of coronary heart disease12, 13, is in strong LD (r2=0.97, D′=1) with the SNP recently reported to be associated with asthma (rs4129267)14. Interestingly, the risk allele at the asthma associated SNP (OR=1.09, p=2.4×10−8) is protective for rheumatoid arthritis (OR=0.9, p = 1.3×10−8). The IL6R ligand, IL6, is the target of the biologic drug, tocilizumab, which has been shown to be an effective treatment for rheumatoid arthritis. Abatacept is another biologic drug, with therapeutic efficacy in clinical trials and which targets another rheumatoid arthritis susceptibility gene, CTLA4. These examples highlight the potential for targeting genes within risk loci.

From 38 rheumatoid arthritis associated SNPs or proxies accessed for eQTL analysis, 18 showed an eQTL effect on at least one probe, giving a total of 51 SNP-probe combinations with significant eQTL effect (Supplementary Table 13). From these 18 SNPs, 11 showed an independent or primary eQTL effect on one or more probes (20 SNP-probe combinations), whereas 7 SNPs were not significant after conditioning of the strongest eQTL signal in the locus.

Using a previously described approach, we assessed whether the 46 independent rheumatoid arthritis associated regions, defined by previously known and novel SNP associations discovered here, harboured genes that were specifically expressed in distinct immune cell-types16. We observed in a large expression data set of 223 sorted mouse immune cells17, that these regions contained genes that were most significantly more specifically expressed in CD4+ effector memory T-cell subsets (p<10−7) (Supplementary Figure 4).

Of the diseases sharing susceptibility loci with rheumatoid arthritis, systematic fine mapping has only been published, to date, for celiac disease3. Previously the two diseases were found to share 6 confirmed non-HLA loci (MMEL, REL, CD28/CTLA4, TNFAIP3, TAGAP and IL2/21) 2; Immunochip data now identifies an additional 4 confirmed loci common to both diseases (DDX6, STAT4, PRKCQ and IRAK1) and a further 4 potential rheumatoid arthritis loci (BACH2, p=8.2×10−7, ELMO1, p=2.9×10−7, PTPN2, p=4.4×10−6, PVT1, p=2×10−5) in common with confirmed celiac disease loci (Supplementary Table 14). Of the ten rheumatoid arthritis/celiac disease loci, 4 share the same lead SNP (CD28, IL2_21, TNFAIP3 and IRAK1) and a fifth (MMEL1) shares highly correlated SNPs (r2>0.88) and, for all of these variants, the risk allele is the same for both diseases. For two loci (PRKCQ and DDX6), the lead SNPs are only moderately correlated (r2>0.62) with the minor allele being protective in both diseases. The effects in STAT4 appear quite different with 3 independent effects in celiac disease and two different independent associations in rheumatoid arthritis. The strongest association signal for risk of celiac disease at TAGAP is with the minor allele of a SNP (rs182429) in moderate LD (r2=0.44) with the rheumatoid arthritis risk SNP (rs629326). Indeed, when considering overlap of rheumatoid arthritis susceptibility loci with other autoimmune diseases, only the PADI4 and CCL21 non-MHC loci currently show unique association, suggesting that they may be important in determining that the autoimmune reaction is directed at synovial joints.

In summary, through fine mapping on a custom made array designed to capture variation across a number of loci associated with autoimmune diseases, we have identified 14 novel European ancestry rheumatoid arthritis loci; refined the peak of association to a single gene at 19 loci, identified 7 SNPs which might potentially be functional, found independent effects at 6 loci and detected association with SNPs with low MAF (<0.05) at 4 loci. In one third of cases, imputation of GWAS signals without fine-mapping, would have implicated a different genetic region as being disease causal thus illustrating the importance of dense fine mapping analysis prior to embarking on expensive functional studies.

Methods

Genotyping

All samples were genotyped for the Immunochip custom array in accordance with Illumina protocols at six centres: UK (Sanger Centre, Hinxton, Cambridge, UK and the University of Virginia, USA), US and Spain (Feinstein Institute, New York, USA), Sweden EIRA (The Genome Institute, Singapore), Sweden Umea (Department of Medical Sciences, SNP&SEQ Technology Platform, Uppsala University Hospital, Uppsala, Sweden) and The Netherlands (Department of Genetics, University Medical Centre Groningen).

Genotype calling and quality control

Genotype calling was performed on all samples at The University of Manchester as a single project using the Genotyping Module (v1.8.4) of the GenomeStudio Data Analysis Software package. Initial genotype clustering was performed using the default Illumina cluster file (Immunochip_Gentrain_June2010.egt) and manifest file Immuno_BeadChip_11419691_B.bpm (NCBI build 36) using the GenTrain2 clustering algorithm. Poor performing samples (call rate < 0.90), labelled duplicates (selection informed by 10th percentile GenCall score (p10 GC)) and samples identified post-genotyping as inappropriate for inclusion were also excluded at this point (Supplementary table 1). Automated reclustering was performed on all remaining samples to calibrate clusters on the study sample set.

Poor quality assays were excluded prior to downstream quality control processes by extensive manual review of clustering performance. A subset of good quality SNPs was identified based on the ranking of quality metrics: cluster separation (<0.4), signal intensity (<1.0), call rate (<0.98) and allele frequency. In addition, SNPs that mapped to the Y chromosome or mitochondria, were non-polymorphic, were duplicates, or zeroed in the default Illumina cluster file were also excluded. This resulted in a dataset of 165,549 good quality SNPs (Supplementary Table 2).

To facilitate the meta-analysis and reduce differential missingness each of the six population datasets were processed as discrete entities. SNPs were excluded from each of the datasets with a call rate < 0.99 (cases or controls), a MAF < 0.01 or if they deviated from HWE (p < 5.7×10−7). Samples were excluded with a call rate < 0.99 or if they were identified as outliers based on autosomal heterozygosity (Supplementary Table 3 and 4). Samples were also excluded if they were considered to be outliers based on ethnicity inferred by principal component analysis (PCA). PCA was performed using EIGENSOFT v4.2 with HapMap phase 2 samples as reference populations on a subset of SNPs with a MAF > 0.05 and filtered to minimise inter-marker LD (excluding the MHC region, 23 regions of high LD and previously confirmed rheumatoid arthritis susceptibility regions) (Supplementary Figure 1). Cryptic relatedness was assessed within each dataset by calculating identity-by-descent (IBD) using PLINK v1.07 using the PCA SNP set. A single sample from any related pair (PI_HAT > 0.1875) was removed from the analysis (informed by call rate). In addition IBD was inferred across all six datasets to exclude cross-dataset related individuals (Supplementary table 5). The genomic control inflation factor (λGC) was calculated within each Immunochip dataset using SNPs included as deep replication for a study investigating the genetic basis for reading and writing ability (submitted by J.C. Barrett). This set of SNPs was filtered as described for the PCA SNP set, leaving a total of 1,469 SNPs distributed evenly across the genome. The λGC for the datasets was estimated at; 1.07 (UK), 1.03 (US), 0.97 (SE-E), 0.94 (SE-U), 1.12 (NL), and 1.10 (ES). Using the same SNPs to estimate λGC1000, where the factor is scaled to the equivalent of 1000 cases and 1000 controls, in the Immunochip meta-analysis resulted in a rescaled λ of 1.02 (1.23 without rescaling).

All novel findings remained significantly associated when including gender and λGC as a covariate in the analysis (Supplementary Table 15).

Immunochip meta-analysis

Association statistics were calculated in each dataset using logistic regression under an additive model (SNPs coded 0, 1 or 2 with respect to minor allele dosage) and incorporating the top ten principal components as covariates. Odds ratios and standard errors were combined across the six datasets using inverse-variance meta-analysis assuming a fixed effect.

Independent effects

Initial evidence for secondary effects was assessed at each of the previously known and newly identified loci using a forward stepwise logistic regression. The index SNP at each region was included as a covariate and the association statistics re-calculated for the remaining test SNPs. This process was repeated until no SNPs reached the minimum level of significance. The criteria for declaring an independent effect was defined as: p-value < 5×10−4, not highly correlated with index SNP, the conditioned p-value must not differ substantially from the unconditioned value. We next tested if the two-SNP fitted the risk at the locus significantly better than the one-SNP model using a likelihood ratio test.

The effect estimates for each two-SNP haplotype was calculated by including indicator variables for carriage of haplotypes. The indicator variables were constructed by phasing the genotype data for each region satisfying the above criteria were phased using the SHAPEIT algorithm18.

GWAS meta-analysis

GWAS case-controls collections were previously described1.Six collections were included in the present study: BRASS, CANADA, EIRA, NARAC1, NARAC2, WTCCC. After quality control and data filtering, the datasets were imputed using IMPUTE and haplotype-phased HapMap Phase 2 European CEU founders as a reference panel19.

We used IBS estimates to remove related samples across the Immunochip and GWAS collections, using GWAS genotype data instead of imputed data. In each of the twelve collections, we selected a set of SNPs with missing-genotype rate<0.5%, minor allele frequency>5% and Hardy-Weinberg PHWE>5×10−7. Then, we extracted SNPs that passed these filters and were shared between the 12 collections. After further LD pruning and resolving flipping issues, the data from the 12 collections were merged to calculate the IBS statistics. When related samples were identified (siblings or duplicates), the sample from the GWAS collection was removed to preferentially keep Immunochip data in the subsequent association analyses. Filtering and IBS calculation were performed using PLINK20. Two GWAS datasets, EIRA and NARAC1, were excluded because of strong overlap (>90% rheumatoid arthritis cases) with the Immunochip SE-E and US collections, respectively. This resulted in a total sample size of 13,838 rheumatoid arthritis cases and 33,742 controls, distributed in 10 collections (Table 1).

The software SNPTEST v2.2 was used to conduct logistic regression analysis of rheumatoid arthritis case-control status in each GWAS collection, conditioning upon the 5 first eigenvectors from PCA analysis, and after excluding SNPs with low statistical information (info score<0.7) or MAF<1%. We also excluded SNPs that were not represented in the filtered Immunochip data. The λGC for the individual datasets was estimated at; 1.04 (BRASS), 1.02 (CANADA), 1.04 (NARAC2) and 1.05 (WTCCC). There was a slight inflation in λGC in these cohorts when using the 1,469 SNPs included on Immunochip to investigate the genetic basis of reading and writing ability; 1.11 (BRASS), 1.15 (CANADA), 1.07 (NARAC2) and 1.05 (WTCCC).

We conducted an inverse-variance weighted meta-analysis to combine the results across the 10 collections. We also computed Cochran’s Q statistics and I2 statistics to assess heterogeneity across collections. Meta-analysis and heterogeneity statistics computation was adapted from the MANTEL program21.

Serological subtype statistical analysis

Multinomial logistic regression was applied to compute odds ratios (OR), 95% confidence interval and p-values for association between the minor allele at every locus and either ACPA-positive (ORACPA-positive) or ACPA-negative rheumatoid arthritis (ORACPA-negative) assuming additivity on the log-odds scale (i.e. every locus was coded as 0,1 or 2 corresponding to the copy number of the minor allele). The minor allele was defined according to the allele frequency in the total population, including cases and controls. To test for differences between ORACPA-positive and ORACPA-negative, the linear combination β+ - β−, where β+ is log (ORACPA-positive) and β− is log (ORACPA-negative) was calculated, along with its standard error. This enables a p-value for the difference in association to be calculated.

We performed an analysis of epistasis using the most significantly associated SNP from each of the 46 loci (Table 2). Logistic regression was performed in PLINK to model epistasis in each of the the six datasets with the top 10 PCs included as covariates. For each pair of SNPs, the likelihood ratio test was employed to compute the p-value of the interaction term for each dataset. Epistasis results were combined using METAL and Bonferroni corrected.

eQTL analysis

eQTL analysis was done on the peripheral blood of 1,469 unrelated individuals (1,240 samples run on the Illumina HT12v3 platform, 229 samples run on the Illumina H8v2 platform) from the United Kingdom and the Netherlands. Details of the eQTL analysis have been previously described22 . In short, we assessed the effect of all rheumatoid arthritis associated SNPs (Table 2) on expression of genes, located within 250kb left and right from the SNP (cis eQTLs).

All individuals from the eQTL study were genotyped on Illumina Hap300K platform and then imputed to HapMap 2 using Impute 2.0 software. Since not all SNPs from Illumina Immunochip platform were genotyped or imputed on the 1,469 eQTL samples, we used the following strategy (Supplementary Figure 5): First, we investigated whether the SNP is present in the eQTL data and had passed the QC for eQTL mapping (MAF >= 5%, HWE P-value >= 0.001, call rate >= 95%). From 50 rheumatoid arthritis-SNPs, 26 were present in HapMap imputed datasets and were directly assessed for eQTL effects (Supplementary Table 13). For the other 24 SNPs, not present in our HapMap imputed data, we checked whether the rheumatoid arthritis-SNP was available in 1000 genomes database. If so, we queried all SNPs within 10MB of the rheumatoid arthritis-SNP that were also present in the eQTL data and would pass eQTL QC measures, and picked the SNP with the highest LD present in HapMap after QC. The threshold of r2>0.8 for the LD was used. For 12 SNPs, no proxy was available with our criteria, and these SNPs were not included in the eQTL analysis. For the remaining 12 SNPs the best proxy SNP is included to the eQTL table (Supplementary Figure 5).

We also performed a cis-eQTL analysis for the top associated gene expression probe, as well as two conditional analyses: (1) conditioning on the effect of the rheumatoid arthritis-SNP (gSNP), and (2) conditioning on the effect of the top eQTL SNP (eSNP) (Supplementary Table 13).

The rheumatoid arthritis associated SNP was labelled as having a primary effect on gene expression if it was either the top eQTL in the locus, or was a good proxy of top eQTL SNP (r2>8). It was labelled as an independent eQTL if it showed an effect after conditioning on the primary eQTL. From 20 rheumatoid arthritis SNPs, that showed an eQTL effect, 13 had either an independent or primary eQTL effect on one or more probes (22 SNP-probe combinations). A further 7 SNPs were not significant after the conditioning of the strongest eQTL signal in the locus, suggesting that they are not primary eQTLs.

We thank Jeffrey Barrett and Chris Wallace for the SNP selection. We would like to thank the WTSI Genotyping Facility and in particular Emma Gray, Sue Bumpstead, Doug Simpkin and Hannah Blackburn. Genotyping of the United Kingdom Rheumatoid Arthritis Genetics samples was supported by the Arthritis Research UK grant reference number 17552 and by the Manchester Biomedical Research Centre. This work was made possible by funds from the Arthritis Foundation (PI SR) and the National Institutes Health (K08AR055688 to S.R. and 1R01AR062886-01 to P.I.W.d.B).Paul Gilbert prepared the UK samples. Genotyping of the Swedish Umea samples was performed by the SNP&SEQ Technology Platform in Uppsala, which is supported by Uppsala University, Uppsala University Hospital, Science for Life Laboratory - Uppsala and the Swedish Research Council (Contracts 80576801 and 70374401). This work was partially supported by the RETICS Program, RD08/0075 (RIER), from Instituto de Salud Carlos III, Spain. We acknowledge use of DNA from The UK Blood Services collection of Common Controls (UKBS-CC collection), which is funded by the Wellcome Trust grant 076113/C/04/Z and by US National Institute for Health Research program grant to the National Health Service Blood and Transplant (RP-PG-0310-1002). We acknowledge the use of DNA from the British 1958 Birth Cohort collection, which is funded by the UK Medical Research Council grant G0000934 and the Wellcome Trust grant 068545/Z/02. The NARAC and analysis of other U.S. patient and control collections at the Feinstein Institute were supported by the National Institutes of Health RO1-AR-4-4422, NO1-AR-2-2263; NO1-AR1-2256, RO1 AI068759, RC2AR059092-01 in addition to support from the Eileen Ludwig Greenland Center for Rheumatoid Arthritis and the family of Robert S. Boas.

Fehrmann RS,et al. Trans-eQTLs reveal that independent genetic variants associated with a complex phenotype converge on intermediate genes, with a major role for the HLAPLoS. GenetYear: 20117e100219721829388