Validation of GWAS-Identified Variants for Anti-TNF Drug Response in Rheumatoid Arthritis: A Meta-Analysis of Two Large Cohorts

We aimed to validate the association of 28 GWAS-identified genetic variants for response to TNF inhibitors (TNFi) in a discovery cohort of 1361 rheumatoid arthritis (RA) patients monitored in routine care and ascertained through the REPAIR consortium and DANBIO registry. We genotyped selected markers and evaluated their association with response to TNFi after 6 months of treatment according to the change in disease activity score 28 (ΔDAS28). Next, we confirmed the most interesting results through meta-analysis of our data with those from the DREAM cohort that included 706 RA patients treated with TNFi. The meta-analysis of the discovery cohort and DREAM registry including 2067 RA patients revealed an overall association of the LINC02549 rs7767069 SNP with a lower improvement in DAS28 that remained significant after correction for multiple testing (per-allele ORMeta=0.83, P Meta=0.000077; P Het=0.61). In addition, we found that each copy of the LRRC55 rs717117G allele was significantly associated with lower improvement in DAS28 in rheumatoid factor (RF)-positive patients (per-allele ORMeta=0.67, P=0.00058; P Het=0.06) whereas an opposite but not significant effect was detected in RF-negative subjects (per-allele ORMeta=1.38, P=0.10; P Het=0.45; P Interaction=0.00028). Interestingly, although the identified associations did not survive multiple testing correction, the meta-analysis also showed overall and RF-specific associations for the MAFB rs6071980 and CNTN5 rs1813443 SNPs with decreased changes in DAS28 (per-allele ORMeta_rs6071980 = 0.85, P=0.0059; P Het=0.63 and ORMeta_rs1813443_RF+=0.81, P=0.0059; P Het=0.69 and ORMeta_rs1813443_RF-=1.00, P=0.99; P Het=0.12; P Interaction=0.032). Mechanistically, we found that subjects carrying the LINC02549 rs7767069T allele had significantly increased numbers of CD45RO+CD45RA+ T cells (P=0.000025) whereas carriers of the LINC02549 rs7767069T/T genotype showed significantly increased levels of soluble scavengers CD5 and CD6 in serum (P=0.00037 and P=0.00041). In addition, carriers of the LRRC55 rs717117G allele showed decreased production of IL6 after stimulation of PBMCs with B burgdorferi and E coli bacteria (P=0.00046 and P=0.00044), which suggested a reduced IL6-mediated anti-inflammatory effect of this marker to worsen the response to TNFi. In conclusion, this study confirmed the influence of the LINC02549 and LRRC55 loci to determine the response to TNFi in RA patients and suggested a weak effect of the MAFB and CNTN5 loci that need to be further investigated.


INTRODUCTION
Rheumatoid Arthritis (RA) is a complex and chronic disease marked by symptoms of inflammation and pain in the joints that eventually lead to joint destruction, loss of function and disability. These symptoms of inflammation are mostly driven by certain central cytokines that modulate both cellular and humoral immune responses in the synovial fluid and synovium of patients (1). Although RA remains as a chronic and incurable autoimmune disease that occurs in as much as 0.5-1% of the general population (2), the introduction of biological agents to target deregulated cytokines has substantially improved the signs and symptoms of the disease (3). Among these cytokines, tumor necrosis factor alpha (TNFa) has attracted most attention as it has been found to be deregulated in patients with autoimmune diseases including RA (4). It has been reported, for instance, that TNFa activates macrophages, synoviocytes, chondrocytes, and osteoclasts in a dose-dependent manner (5) and that high levels of circulating TNFa correlate with disease activity and disease progression (6).
The increasing number of biological agents approved by the FDA and the increased prevalence of the disease all around the world (7) have placed a substantial economic burden for health care systems. Although the introduction of biosimilars in clinical practice reduced the cost of these treatments in many countries (8), there is still an unmet need to optimize biologic therapies, avoiding unnecessary adverse effects risks and reducing costs (9).
The interplay between genetics and drug response has been the subject of intense investigations during last decades. Response to biologics has been shown to vary between individuals and that a large proportion of patients show no clinical improvement (10). Given the high cost of these drugs and the potential impairment of non-responding patients, the identification of genetic biomarkers associated with drug response to specific biological agents would help to know which patients might benefit from a particular treatment. However, to date, only a few genome-wide association studies (GWAS) (11)(12)(13)(14)(15)(16)(17) or well powered candidate gene association studies have been conducted (18)(19)(20)(21)(22)(23)(24)(25)(26). We are far from being able to optimize drug dosing or prioritize drug combinations based on genetic findings. In fact, attempts to validate the association of most of the genetic markers identified in these association studies have failed (27), which confirms the limited application of genetic findings in a clinical setting. Considering that the validation of previous GWAS findings is an essential step to tailor treatments for RA and to approach personalized medicine, we aimed to validate the association of GWAS-identified variants for response to TNF inhibitors (TNFi) in a two-stage nested casecontrol association study including a cohort of 1361 anti-TNF naïve RA patients ascertained through the REPAIR consortium and DANBIO registry and an independent replication cohort of 706 RA patients treated with TNFi from the DREAM registry. We also investigated whether the effect of selected markers on the response to TNFi could be modified by rheumatoid factor (RF) status, and whether genetic variants could influence immune responses and affect the serological concentration of 108 plasmatic inflammatory proteins, 7 serum steroid hormones or counts of 91 blood-derived immune cell populations.

Study Populations and Response to TNFi
The discovery population consisted of 1361 RA patients ascertained through the REPAIR consortium and the DANBIO registry (Table 1) (28,29). RA patients fulfilled the 1987 revised American College of Rheumatology (ACR) (30) and/or the ACR/ EULAR 2010 classification criteria (31). In order to further replicate the most interesting results, we validated the association with response to anti-TNF drugs of those SNPs showing a P<0.05 in the discovery cohort in 706 Dutch RA patients treated with TNFi from the DREAM (Dutch RhEumatoid Arthritis Monitoring) registry (Supplementary Table 1). The study followed the Declaration of Helsinki. Study participants were of European origin and gave their written informed consent to participate in the study, which was approved by the ethical review committee of participant institutions: Virgen de las Nieves University Hospital (2012/89); Santa Maria Hospital-CHLN (CE 877/121.2012); University Clinical Hospital of Santiago de Compostela (2013/156); Wroclaw Medical University (KB-625/2016); Radboud university medical center (2011/299) and by the Regional Ethics Committee of Central Denmark Region (S-20120113). A detailed description of the discovery population has been reported elsewhere (19,20,22,24). All RA patients were naïve for TNFi and response to TNFi for each patient in all study populations was calculated using the change in disease activity score (DAS28CRP) between baseline and 6 months after treatment. Overall and RF-stratified linear regression analyses adjusted for age, sex and country of origin were used to determine the association between GWAS-identified SNPs and changes in DAS28. RA patients with missing values either for DAS28 (in at the time points of interest) or RF were not included in the analysis.

DNA Extraction, SNP Selection, Genotyping, and Quality Control
Genomic DNA from RA patients was extracted from blood samples using the QIAamp DNA Blood Mini kit (Valencia, CA, EEUU) according to manufacturer's instructions. The single-nucleotide polymorphisms (SNPs) were selected through an extensive literature search of relevant GWAS and metaanalyses published by February 2019 using publicly available online databases. Additional criteria were potential functionality and linkage disequilibrium between the reported SNPs. Biological function was predicted according to the data publicly available in the integrated Regulome database (www. regulomedb.org), and eQTL data browsers (www.gtexportal.org/ home/ and https://genenetwork.nl/bloodeqtlbrowser/). A total of 28 SNPs in 25 genes were selected for genotyping in the discovery cohort ( Table 2). Genotyping of selected SNPs was performed using KASP ® probes according to manufacturer's instructions (LGC Genomics, Hoddesdon, UK). For quality control,~5% of DNA samples were randomly included as duplicates and concordance between duplicate samples was ≥99.0%. Replication of the most interesting association was conducted in the DREAM registry (n=706) following a similar quality control genotyping strategy.

Hardy-Weinberg Equilibrium, Genetic Association Analysis, and Meta-Analysis
Deviation from Hardy-Weinberg Equilibrium (HWE) was tested in the control group (responders and moderate responders according to the EULAR response criteria) by chi-square (c2), and linear regression analysis adjusted for age, sex and country of origin was used to assess the associations of the GWAS-identified polymorphisms with absolute changes in DAS28 assuming logadditive, dominant and recessive models of inheritance. Those SNPs with the lowest P-value in the discovery population according to each genetic model were advanced for replication in the DREAM cohort and meta-analysis of the discovery and replication populations using a fixed effect model was performed to validate the association observed. I2 statistic was used to assess heterogeneity between studies. Correction for multiple testing was performed using the Bonferroni method but also considering the two inheritance models tested. Given that log-additive and dominant models showed a high degree of collinearity, the significant threshold for the meta-analysis was set to 0.00089 considering log-additive/dominant and recessive inheritance models (0.05/28SNPs/2models). Overall statistical power was calculated using Quanto (v.12.4) assuming a log-additive model and a baseline risk of 30% for response to TNFi (32,33).

Cell Isolation, Differentiation and Cytokine Quantitative Trait Loci, and Hormone Analysis in Relation to the GWAS-Identified Variants for Response To TNFi
With the aim of determining whether those SNPs associated with response to TNFi had a role in modulating immune responses, we performed in vitro stimulation experiments and measured cytokine production (IFNg, IL1Ra, IL1b, IL6, IL8, IL10, TNFa, IL17, and IL22) after stimulation of peripheral blood mononuclear cells (PBMCs), whole blood or monocyte-derived macrophages (MDMs) from 408 healthy subjects of the 500FG cohort from the Human Functional Genomics Project (HFGP) with LPS (1 or 100 ng/ml), PHA (10mg/ml), Pam3Cys (10mg/ml), CpG (ODN M362; 10mg/ml) and B. burgdorferi and E. coli, as experimental model for cytokine production capacity. Given the sex disparities in the prevalence and course of RA and the impact of steroid hormones in modulating immune responses, we also evaluated the correlation of SNPs with serum levels of 7 steroid hormones (androstenedione, cortisol, 11-deoxy-cortisol, 17hydroxy progesterone, progesterone, testosterone and 25 hydroxy vitamin D3) in a subset of the 500FG cohort without hormonal replacement therapy or oral contraceptives (n=280). After log transformation, cytokine or serum steroid hormone levels were correlated with the SNPs of interest using a linear regression model with age and sex as co-factors in R (http:// www.r-project.org/). This analysis led to cytokine quantitative trait loci (cQTL) and hormone quantitative trait loci (hQTL). Significance thresholds were set to be 0.000463 and 0.00357 (0.05/6stimulants/9cytokines or 0.05/7hormones and 2 inheritance models) for cQTL and hQTL, respectively.

Correlation Between GWAS-Identified Polymorphisms and Cell Counts of 91 Blood-Derived Immune Cell Populations and Serum/Plasmatic Proteomic Profile
We also investigated whether selected polymorphisms had an impact on blood cell counts by analyzing a set of 91 manually annotated immune cell populations and genotype data from the 500FG cohort that consisted of 408 healthy subjects (Supplementary Table 2). Cell populations were measured by 10-color flow cytometry (Navios flow cytometer, Beckman IFNK||C9orf72 Coulter) after blood sampling (2-3 hours) and cell count analysis was performed using the Kaluza software (Beckman Coulter, v.1.3). In order to reduce inter-experimental noise and increase statistical power, cell count analysis was performed by calculating parental and grandparental percentages, which were defined as the percentage of a certain cell type within the cellpopulations one or two levels higher in the hierarchical definitions of cell sub-populations (34). Detailed laboratory protocols for cell isolation, reagents, gating and flow cytometry analysis have been reported elsewhere (35) and the accession number for the raw flow cytometry data and analyzed data files are available upon request to the authors (http://hfgp.bbmri.nl).
A proteomic analysis was also performed in serum and plasma samples from the 500FG cohort. Circulating proteins were measured using the commercially Olink ® Inflammation panel (Olink, Sweden) that resulted in the measurement of 103 different biomarkers (Supplementary Table 3). Proteins levels were expressed on a log2-scale as normalized protein expression values, and normalized using bridging samples to correct for batch variation. Considering the number of proteins (n=103) and cell populations (n=91) tested, P-values of 0.00049 and 0.00055 were set as significant thresholds for the proteomic and cell-level variation analysis, respectively.

RESULTS
A total of 1361 anti-TNF patients were included in the discovery population. The mean age of the RA patients was 52±14 and they showed a female/male ratio of 3.4 (1050/310). Sixty-seven percent of the RA patients were positive for RF and 64% had anti-citrullinated protein antibodies (ACPA). The median disease duration was of 12.92 years and the disease activity score 28 (DAS28CRP) calculated at patient recruitment was of 5.91 ( Table 1).

Association of GWAS-Identified SNPs With
Response to Anti-TNF Drugs

Functional Characterization of the Most Interesting Findings
Considering these results, we attempted to shed some light into the functional consequences of the overall or RF-specific effects of the LINC02549 rs7767069 , LRRC55 rs717117 , MAFB rs6071980 and CNTN5 rs1813443 SNPs to modulate the response to TNFi. Interestingly, our functional experiments showed that, when considering the total number of leukocytes as reference, the LINC02549 rs7767069 polymorphism significantly correlated with increased numbers of CD45RO+CD45RA+ T cells in blood (P=0.00047; Figure 1A). Subjects carrying the LINC02549 rs7767069T allele (associated with poor response to TNFi in RA patients) had significantly increased numbers of CD45RO+CD45RA+ T cells (P=0.000025), which suggested that this genetic marker might influence the response to TNFi by mediating the number of this specific T cell subset in blood and, thereby contribute to inflammation. In addition, we observed that those subjects carrying two copies of the LINC02549 rs7767069T allele showed significantly increased serum levels of soluble scavenger receptors CD5 and CD6 when compared with those carrying the A/T or A/A genotypes (P=0.00037 and P=0.00041; Figures 1B, C). These results also suggested a functional role of the LINC02549 rs7767069 SNP in RA likely through the CD5/CD6-mediated modulation of T cells and certain subsets of B cells that control multiple processes including cellular adhesion and migration across endothelial and epithelial cells, antigen presentation by B cells and the subsequent proliferation of T cells.
Furthermore, although we could not stratify our functional analyses by RF status because of the healthy nature of the blood donors, we found that carriers of the LRRC55 rs717117G allele showed significantly decreased levels of IL6 production after stimulation of PBMCs with either B. burgdorferi (P=0.00046; Figure 2A) or E. coli (P=0.00044; Figure 2B), which suggested an implication of the LRRC55 locus in the modulation of the response to TNFi by regulating IL6 production and likely IL6mediated T cell differentiation into effector Th2 cells. Functional data from Haploreg also showed that the LRRC55 rs717117 variant correlates with mRNA P2RX3 expression levels, a well-known gene involved in controlling T cell proliferation. Finally, our functional experiments revealed that carriers of the MAFB rs6071980 C allele showed decreased levels of Chemokine (C-C motif) ligand 23 (CCL23; P=0.0060; Figure 3A) and increased levels of serum Fibroblast growth factor 19 (FGF-19; P=0.0034; Figure 3B). Whereas FGF-19 protein modulates inflammation by mediating IL6 production, CCL23 has been implicated in monocyte recruitment during inflammation and it has been previously shown to positively correlate with drop in DAS28 after treatment with TNFi. Although the effect of the MAFB SNP to modulate either serum FGF-19 or CCL23 levels did not remain significant after correction for multiple testing, these results might indicate a weak, but still functional, effect of the MAFB locus in modulating response to anti-TNF drugs.

DISCUSSION
Although treatment of RA patients using monoclonal anti-TNF drugs has been a particularly successful approach to control inflammation and to prevent joint destruction and the appearance of bone erosions, non-responsiveness is prevalent and no effective biomarkers for drug response prediction have been consistently identified (36). This comprehensive validation study aimed at confirming the association of GWAS-identified variants with response to TNFi and to shed some light into the biological mechanisms underlying the most interesting associations. For that purpose, we conducted a two-stage case control study including 2067 RA patients treated with anti-TNF drugs ascertained through the REPAIR consortium but also DANBIO and DREAM registries. The most significant result was the overall association of the LINC02549 rs7767069 SNP with a poor response to anti-TNF drugs. The meta-analysis of the discovery and replication cohorts showed that each copy of the LINC02549 rs7767069T allele significantly decreased the improvement in DAS28 by 17% after the treatment with a TNFi. Importantly, the association of the LINC02549 rs7767069 variant with poor response to TNFi was significant in the two populations analyzed and remained significant after correction for multiple testing, which confirmed a role of the LINC02549 locus in the modulation of response to anti-TNF drugs. LINC02549 (Long Intergenic Protein Coding RNA 2549) is an RNA gene that is affiliated with the lncRNA class, which represents a large proportion of the human transcriptome. LINC02549 maps to chromosome 6 and it is expressed in resting T cells and CD4 activated T cells. Although its function is still largely unknown, our data suggest that it might exert a role in determining the number of circulating CD45RO+CD45RA+ T cells, which are a subset of cells frequently found in the synovial fluid of both chronic arthritis (37) and RA patients (38). According to the results of Koch et al. (1990), CD45RA+ CD45RO+ T lymphocytes are mostly detected in perivascular regions, which suggest that these lymphocytes might access the RA synovial tissue via the synovial vasculature (38) and that, once there, they could play a role in promoting synovial tissue inflammation mainly by the induction of memory immune responses. Therefore, it seems to be plausible to suggest that the negative impact of the LINC02549 rs7767069 SNP on the response to TNFi might be mediated by its role in modulating numbers of CD45RA+CD45RO+ T lymphocytes that could migrate to the synovial tissue and promote inflammatory responses and, thereby hamper the control of inflammation during treatment with anti-TNF drugs. In support of this hypothesis, we found that carriers of two copies of the LINC02549 rs7767069T allele also showed significantly increased levels of soluble scavenger receptors CD5 and CD6 (sCD5 and sCD6) in serum that are proteins highly expressed in regulatory T cells and a specific subset of B cells (CD5+ or B1a) (39). Although the origin of these soluble scavenger receptors in RA is poorly understood, it has been suggested that they are shed in the serum by proteases from the surface of activated lymphocytes that subsequently infiltrate synovium structures (40). In fact, increased serum levels of sCD5 and sCD6 has been found in subjects diagnosed with RA (41-44) but also other autoimmune diseases such as primary Sjögren's syndrome (42,45,46), systemic inflammatory response syndrome (47), multiple sclerosis (44) or dermatitis (48). Although the functional role of both soluble scavengers in autoimmune diseases is still under investigation, it is well established that sCD5 and sCD6 are regulators of T cell functions and induce autoreactivity. It is known that they are required for the initiation, differentiation and maintenance of T cell immune responses (49,50) but also T cell migration and extravasation to the synovial tissue (51). Furthermore, it has been reported that both sCD5 and sCD6 are involved in the modulation of TCR and BCR signaling and determinate T-and B-cell survival (52) and Th17 differentiation (53,54). Furthermore, clinical trials using humanized anti-CD6 mAbs have provided valuable information regarding the potential targeting of CD6 for the treatment of RA but also psoriasis and potentially other T cell-driven autoimmune diseases (55)(56)(57). Recent investigations have also suggested that genetic alterations within the CD6 gene associated with clinical outcome of several autoimmune diseases (44,58) and correlated with the response to TNFi (59), which pointed to a role of these soluble scavenger receptors in modulating response to anti-TNF drugs. Considering these findings, we hypothesize that, besides its effect on modulating number of the CD45RA+CD45RO+ T lymphocytes, the LINC02549 rs7767069 SNP might negatively influence the response to anti-TNF drugs by stimulating directly or indirectly the production of sCD5 and sCD6 and thereby inducing long-term T cell-mediated immune responses. Another interesting result that remained significant after correction for multiple testing was the RF-specific association of the LRRC55 rs717117 SNP with lower changes in DAS28 after the treatment with TNFi. The meta-analysis of the discovery and replication cohorts showed that RF-positive patients carrying the LRRC55 rs717117G allele have a significantly decreased drop in DAS28 after treatment with a TNFi, whereas an opposite but not statistically significant effect was observed in RF-negative RA patients. Noticeably, functional experiments showed that, after stimulation of PBMCs from healthy subjects with B. burgdorferi and E. coli bacteria, carriers of the LRRC55 rs717117G allele showed significantly decreased production of IL6 when compared to those carrying the most common genotype. Although functional experiments could not be stratified by RF because of the healthy nature of blood donors, these results suggested a role of the LRRC55 locus in modulating IL6-mediated immune responses. On the other hand, functional data from Haploreg also suggested an implication of the LRRC55 rs717117 variant in controlling P2RX3-mediated T cell proliferation.
LRRC55 gene maps on chromosome 11 and it encodes for the leucine-rich repeat-containing protein 55, a protein that belongs to the LRRC superfamily that include hundreds of proteins mainly expressed in brain. Several LRRC proteins have been linked to the regulation of ion channels (60) but it has been also demonstrated that LRRC proteins are also implicated in modulating immune responses against bacterial pathogens (61) and modulate cell trafficking of membrane receptors such as toll-like receptors (62). Although the interplay between LRRC55 and IL6 has not been demonstrated, our experimental data suggest that the LRRC55 rs717117 SNP modulates IL6 production in response to bacteria and, therefore, might be involved in other IL6-dependent immune processes that could worsen the response to TNFi.
It is widely known that IL6 can induce both antiinflammatory and pro-inflammatory immune responses, which depend entirely on the signalling pathway triggered. Whereas anti-inflammatory responses are mostly mediated by the classic signalling cascade (through binding to the transmembrane IL6 receptor), pro-inflammatory responses and chronic inflammation are mediated by trans-signalling (through binding to the soluble IL6 receptor) or by the interaction of IL6R with gp130 (63,64). Considering our functional data, it is conceivable to suggest that the LRRC55 rs717117 SNP might affect the response to TNFi by decreasing IL6 production and thus inhibiting the classical IL6-dependent anti-inflammatory pathway and dysregulating pro-inflammatory responses. In support of this hypothesis, several mouse models have shown that the activation of the IL6 classic signaling pathway is essential for the activation of STAT3-mediated signaling pathways which reduce inflammation and induce the regeneration of the affected tissues (65). In addition, it has been reported that IL6 is one of the earliest factors that trigger the differentiation of naive T cells into effector Th2 cells in vitro and that, when absent, aggravates the development of the inflammatory processes (64).
Finally, although the genetic association of the MAFB rs6071980 SNP with lower response to TNFi did not remain significant after correction for multiple testing, we found that the MAFB rs6071980 SNP correlated with higher levels of serum FGF-19 and decreased levels of CCL23. Given that FGF-19 is a master protein involved in the inhibition of intestinal inflammation (66,67) and CCL23 has been positively correlated with the DAS28 score in RA patients (68), we think that it would worth to investigate more in detail the impact of this SNP on drug response in future studies. In addition, it might be interesting to further analyze the weak association of the CNTN5 rs1813443 SNP with poor response to TNFi. However, given that we could not find any significant impact of this marker on immune responses, blood cell counts or serum inflammatory proteins or steroid hormones, we are prone to think that this SNP might not have a relevant role in modulating response to TNFi.
This study has both strengths and weaknesses. Among the strengths we can highlight the use of large and well-characterized RA patient populations that allowed the development of a wellpowered overall association analysis but also to investigate the effect modification by RF status. Overall, we had 80% of power to detect an OR of 1.18 (a=0.00089) for a SNP with a frequency of 0.25. On the other hand, it is worth mentioning the comprehensive analysis of the functional effect of the most interesting genetic variants on modulating immune responses, which was performed using a large sample size for this kind of studies. We analysed cQTL and hQTL data but also counts of 91 blood-derived cell populations and serum levels of 103 immunological proteins. An important limitation of this study was the impossibility to adjust linear regression analyses for potential confounding factors including concomitant treatments that might influence the response to TNFi. In addition, given the healthy nature of the subjects included in the HFGP cohort, we could not control our functional experiments by RF status.
In conclusion, this study validates the overall or RF-specific association of LINC02549 and LRRC55 loci with the response to TNFi and provides new insights into the functional role of these polymorphisms in modulating immune responses and response to anti-TNF drugs.

DATA AVAILABILITY STATEMENT
All data used in this project have been meticulously catalogued and archived in the BBMRI-NL data infrastructure (https://hfgp. bbmri.nl/) using the MOLGENIS open-source platform for scientific data (69). This allows flexible data querying and download, including sufficiently rich metadata and interfaces for machine processing (R statistics, REST API) and using FAIR principles to optimize Findability, Accessibility, Interoperability and Reusability (70,71). Genetic data from the discovery and DANBIO populations can be accessed at ftp.genyo.es and data from the DREAM registry are available at https://www.synapse. org/#!Synapse:syn3280809/wiki/194735 and https://www.synapse. org/#!Synapse:syn3280809/wiki/194736.

FUNDING
This study was supported by grants PI17/02276 and PI20/01845 from Fondo de Investigaciones Sanitarias (Madrid, Spain) and by intramural funds of GENYO and FIBAO foundation (Granada, Spain). This study was also supported by the Novo Nordisk Fonden (NNF15OC0016932, VA) and Knud og Edith Eriksens Mindefond (VA) and Gigtforeningen (A2037, A3570, VA). JS and KB-K were supported by the grant No. 2016/21/B/NZ5/ 01901 from the National Science Centre (Poland). MGN was supported by a Spinoza grant from the Netherlands Organization for Scientific Research. YL was supported by an ERC Starting Grant (948207) and the Radboud University Medical Centre Hypatia Grant (2018) for Scientific Research. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
We thank all participants who have agreed to participate in this study. Authors also thank Marıá Dolores Casares, Ángeles Molina, Carmen Oloriz for the collection of Spanish samples and Hans Jurgen Hoffmann, Marianne Thomsen, Vibeke Østergaard Thomsen, Malene Rohr Andersen, Lise Lotte B. Laursen, Helle Jørgensen, Ram Benny Christian Dessau, Niels Steen Krogh, Ulla Vogel, Paal Skytt Andersen, Ivan Brandslund, Steffen Bank, Frederik Trier Møller, Nikolai Toft and Niels Møller Andersen for the participation in collection and purification of Danish samples. We also thank the Danish Departments of Rheumatology for their implication in the collection of clinical data from RA patients included in the DANBIO cohort and the Danish Rheumatologic Biobank. Likewise, we would like to thank Teun van Herwaarden for steroid hormone measurements in serum samples from subjects ascertained through the HFGP initiative.