Evaluation of 12 GWAS-drawn SNPs as biomarkers of rheumatoid arthritis response to TNF inhibitors. A potential SNP association with response to etanercept

Research in rheumatoid arthritis (RA) is increasingly focused on the discovery of biomarkers that could enable personalized treatments. The genetic biomarkers associated with the response to TNF inhibitors (TNFi) are among the most studied. They include 12 SNPs exhibiting promising results in the three largest genome-wide association studies (GWAS). However, they still require further validation. With this aim, we assessed their association with response to TNFi in a replication study, and a meta-analysis summarizing all non-redundant data. The replication involved 755 patients with RA that were treated for the first time with a biologic drug, which was either infliximab (n = 397), etanercept (n = 155) or adalimumab (n = 203). Their DNA samples were successfully genotyped with a single-base extension multiplex method. Lamentably, none of the 12 SNPs was associated with response to the TNFi in the replication study (p > 0.05). However, a drug-stratified exploratory analysis revealed a significant association of the NUBPL rs2378945 SNP with a poor response to etanercept (B = -0.50, 95% CI = -0.82, -0.17, p = 0.003). In addition, the meta-analysis reinforced the previous association of three SNPs: rs2378945, rs12142623, and rs4651370. In contrast, five of the remaining SNPs were less associated than before, and the other four SNPs were no longer associated with the response to treatment. In summary, our results highlight the complexity of the pharmacogenetics of TNFi in RA showing that it could involve a drug-specific component and clarifying the status of the 12 GWAS-drawn SNPs.


Introduction
Rheumatoid arthritis (RA) is a systemic autoimmune disease that until the late 1990s led to permanent disability, low life quality and increased mortality [1]. The development of targeted drugs, pioneered by TNF inhibitors (TNFi), transformed this poor clinical evolution. Now, it is possible to obtain long-term clinical remission or low disease activity in an important proportion of patients [1,2]. The remaining patients (about 30%) will not appropriately respond to a specific drug although they may respond to another. Therefore, biomarkers for prediction of the response will improve the benefits and avoid the unnecessary costs and side effects of the targeted drugs [3,4].
The goal of predicting the response to treatment in RA patients has been pursued in many research areas [3,4]. One of these areas has been genetics, where candidate-gene and genomewide studies (GWAS) have been performed [5,6]. They have been primarily concentrated on the response to three TNFi: infliximab, adalimumab, and etanercept, as the most widely used biologic Disease Modifying Anti-Rheumatic Drug (bDMARD). The initial studies were focused on candidate genes, with many addressing the TNFα gene [7,8]. These studies were small, probably expecting polymorphisms with an important influence in the drug effect [6,9]. Unfortunately, their findings were not reproducible showing the initial expectations were too optimistic [6,8,10-12]. More recently, several large studies have been reported including many hundreds or thousands of RA patients [12][13][14][15][16][17]. They have demonstrated promising SNPs that are associated with the response to TNFi at various levels of evidence. Some appeared in candidate-gene studies, as the PTPRC rs10919563 SNP, which approached the GWAS-level of significance combining three large studies [15][16][17]. Others have been highlighted in GWAS [11][12][13][14]18,19], like the four SNPs we attempted to validate in a previous work [20], and the 12 SNPs that we have selected now.
We have drawn these 12 SNPs from the three largest published GWAS [12][13][14]. Two of them included the same � 2700 patients that were analyzed according to different protocols [12,14], while the third GWAS counted with 1278 patients [13]. The 12 SNPs fulfilled the requirements of replicability established on the respective GWAS, although none of them reached the GWAS-level of significance (p < 5 x10 -8 ). Nevertheless, the CD84 rs6427528 was associated with p = 8 x10 -8 , but only with the response to etanercept, not with the response to infliximab or adalimumab [14]. This result signaled the possibility of drug-specific biomarkers within the response to the TNFi. Indeed, other studies have shown drug-specific genetic [19,[21][22][23] and protein biomarkers [24]. This specificity could be consequence of the known differences in structure, pharmacokinetics and interactions between the three TNFi [25,26]. Therefore, we have addressed the replication of the 12 SNPs considering the three TNFi together and separately. In addition, we have completed the SNPs assessment by meta-analysis to combine our results with the data from previous studies.

Patients
A total of 788 patients with RA according to the American College of Rheumatology classification criteria [27] were included. They were either of self-reported Spanish European ancestry (n = 731) recruited in 15 Spanish Rheumatology Units, or of Greek European ancestry recruited in two Greek hospitals (n = 57). All provided blood samples for DNA extraction and their informed written consent. The study was conducted according to the principles of the Declaration of Helsinki (2013)  . All the patients were treated with a TNFi as the first bDMARD between 2000 and 2011. The indication of treatment, the choice of drug, and the control of clinical evolution were performed with independence of this study during the standard care of the patients. Evaluations included Disease Activity Score 28 (DAS28) at the start of treatment and at 3, 6, and 12 months. DAS28 is a composite index of RA activity including the number of tender joints and swollen joints (28 joints maximum), erythrocyte sedimentation rate, and global patient health status assessment [28]. Patients with baseline DAS28 < 3.2 (i.e. showing low activity, n = 13), and samples failing most genotypes (n = 20) were excluded from further analysis. The remaining 755 patients were distributed as follows: 397 treated with infliximab, 155 treated with etanercept, and 203 treated with adalimumab. Their clinical characteristics are detailed in Table 1. We have sufficient information to analyze response to treatment of 452 patients at 3 months, 689 patients at 6 months and 531 at 12 months. The corresponding raw data are provided in Supporting information S1 Table. Genotyping assays Twelve SNPs (Table 2) were selected because of their reported association with response to TNFi in published GWAs [12][13][14]. These SNPs were genotyped with a multiplex single-base extension technology (SNaPshot Multiplex Kit from Applied Biosystems, Foster City, CA) starting with DNA amplification in a multiplex PCR reaction (KAPA2G fast HotStart, Kapa Biosystems, Woburn MA). Ten per cent of the samples were re-genotyped for quality control. Primers and probes used for these analyses are available in Supporting Information S2 Table.

Statistical analyses
We have used reproducibility of genotypes, genotype call rate, the Hardy-Weinberg equilibrium (HWE) and coincidence with SNP frequencies in the HapMap Toscani in Italy (TSI) collection [29] as quality control measures. The response to TNFi was assessed primarily as change in DAS28 (ΔDAS28 = DAS28 baseline −DAS28 follow-up ) at 6 months of follow-up. Additionally, we have also considered ΔDAS28 at the 3 and 12 month points and the responder (good + moderate) versus non-responder classification according to the European League Against Rheumatism (EULAR) criteria [30]. The EULAR criteria divide patients into three classes based on change in DAS28 from baseline (ΔDAS28) and DAS28 at the time of evaluation: good responders are those with ΔDAS28 � 1.2 and DAS28 � 3.2; non-responders are all patients with ΔDAS28 � 0.6 and those with ΔDAS28 > 0.6 but � 1.2 and with DAS28 > 5.1; all the remaining patients are moderate responders. Generalized linear models for ΔDAS28 and logistic regression models for EULAR response criteria were fitted. Genotypes were considered according to an additive genetic model of minor allele counts (0, 1 or 2). Therefore positive regression coefficients indicate a better response associated with minor allele additive effects. Covariates included in the models were baseline DAS28, gender, the specific TNFi, and the Spanish or Greek origin. Statistica 7.0 (Statsoft, Tulsa OK) software was used to perform these analyses. Meta-analysis of the current study with all the available non-redundant results from previous studies was done. Specifically, information from 2466 patients included in was slightly different than the used to discover new associations. All these patient sets were combined with the fixed effects model of meta-analysis, weighting each cohort by the inverse variance method, except for SNPs showing significant heterogeneity (I 2 > 50%). In this latter case, the random effects model according to DerSimonian and Laird was applied. These analyses were conducted with the R metafor package [31]. Post-hoc power analysis was conducted with G � Power 3 considering the response to treatment at any time [32].

Results
The 755 patients with RA showed characteristics of a severe disease before starting treatment with TNFi (Table 1). This was indicated by frequent erosive arthritis (70.7%), high disease activity (mean DAS28 = 5.8) and moderate to severe disability (mean HAQ = 1.5). The patients In our attempt to validate biomarkers of response to TNFi, we analyzed the 12 SNPs from GWAS in Table 2. Their genotypes passed quality control criteria, including call rate (> 99%), reproducibility (100%), fit to HWE (P > 0.05), and consistency with HapMap frequencies (all pair-wise comparisons P > 0.05). Unfortunately, linear regression did not show any significant association between the 12 SNPs and the predefined main outcome of response, ΔDAS28 at 6 months (Table 3), or with the secondary outcomes, ΔDAS28 at 3 or at 12 months (Table 3). Similarly, logistic regression did not show association of the SNPs with the non-responder vs. responder classification of the patients (S4 Table). In addition, the CD84 SNP rs6427528, which was specifically associated with the response to etanercept in a previous GWAS [14], was not associated in this subset of our patients (B = 0.16, 95% CI -0.43, 0.76, p = 0.6 for ΔDAS28 at 3 months; and OR = 1.16, 95% CI (0.68-1.97), p = 0.6 for the responder vs. nonresponder comparison; and similar results at 6 or 12 months).
Considering the convenience of summarizing the results, we performed post-hoc power analysis of our results, and combined meta-analysis with the previous data. The power analysis showed our study had only enough post-hoc power (1-β > 80%) to reproduce the rs4694890 and rs7767069 reported associations, with rs1350948 very near this level (77%). For the remaining SNPs, power was insufficient and assessment of their status relies more heavily in the combined meta-analysis than in the replication. The meta-analysis (Table 4) included a total of 3155 patients with information for the 3 SNPs identified by Plant et al. [13], and the same number, but not exactly the same patients, for the 8 SNPs drawn from Umicevic Mirkov et al. [12]. In turn, there were a total of 1178 patients treated with etanercept with information for the rs6427528 SNP identified by Cui et al. [14]. Three of the 12 associations were slightly reinforced by the meta-analysis including our results (Fig 1). The association of rs2378945 changed from p = 6.9 x10 -4 to 1.8 x10 -4 ; the corresponding to rs12142623 passed from p = 2.0 x10 -4 to 4.2 x10 -5 ; and the change for rs4651370 was from p = 1.1 x10 -4 to 5.6 x 10 −5 . In contrast, four other SNPs were no longer associated with response to treatment and the remaining five SNPs showed decreased associations (Table 4).
To complete our study, we included an exploratory analysis stratified by the different TNFi. In this exploration, we uncovered association of the NUBPL rs2378945 SNP at 3 months with ΔDAS28 only in the etanercept-treated patients (Fig 2, B = -0.50, 95% CI (-0.82, -0.17), p = 0.003). This result represented less decrease of DAS28 in the patients bearing the minor allele of rs2378945. The association with reduced response was also reflected in the analysis of the EULAR response. In effect, the minor allele of rs2378945 was associated with less good responders compared with moderate responders and non-responders considered together (OR = 0.54, 95% CI = 0.31-0.96, p = 0.035), without significant distinction between moderate responders and non-responders (p = 0.7). The circumscribed association to the 3-month evaluation could be due to the progressive improvement observed in the patients treated with etanercept (non-responders changed from 20.9% at 3 months to 8.3% and 7.2% at 6 and 12 months, respectively). An improvement that was not observed in patients treated with the other TNFi (S3 Table).

Discussion
None of the 12 GWAS-drawn SNPs showed association with the response to TNFi in our RA patients. This lack of replication is important because the SNPs were from high-quality studies and this is the first attempt of independent replication for 9 of them. Nevertheless, the combination of our results with the previous studies discriminated between 3 SNPs that showed strengthened association and the remaining 9 SNPs whose status was weakened without question. In any case, the low level of reproducibility reinforces the need for circumspect consideration and calls for more powerful studies. In this regard, the association of NUBPL rs2378945 with the specific response to etanercept could reflect the phenotype complexity.
Two of the three associations that were reinforced by our results could be considered as only one because rs12142623 and rs4651370 are near and in high linkage disequilibrium  between them (r 2 = 0.90 in subjects of European descent according to the Phase 3 of the 1000 Genomes project; https://ldlink.nci.nih.gov/). They could modify the expression of the neighboring PLA2G4A gene, which codes for the arachidonyl specific group IVA cytosolic phospholipase 2 enzyme. This enzyme regulates TNF-induced joint damaging mediators [33]. The other reinforced association corresponds to rs2378945, which is the same SNP we found associated with reduced response to etanercept in the drug-stratified analysis. This SNP is associated with the expression of the nearby NUBPL gene (GTEx, http://www.gtexportal.org.), which is a member of the ATP-binding protein family required for the assembly of the mitochondrial complex I. However, it has not any known relationship with inflammation or autoimmunity. Three of the four SNPs that were no longer associated after meta-analysis (rs7962316, rs1350948, rs4694890) had a weak standing as biomarkers of response to TNFi even before our study. They were associated in the two first phases, but not in the third, of the original GWAS highlighting them [13]. Afterward, their status was further undermined by the lack of replication in two subsequent GWAS (although the time of assessment was different) [11,12]. Therefore, our results regarding these 3 SNPs could be taken as a confirmation of their lack of association with response to TNFi. Also, the 8 SNPs drawn from the Umicevic Mirkov et al.
[12] report had been associated with response to TNFi in the first two phases of the original GWAS, but not in the additional replication collections from the same study [12]. The particularity with these 8 SNPs is that no other work had addressed them until now. Therefore, our study is the first independent validation and it differentiated between the three SNPs that were reinforced and the five SNPs that were more clearly questioned. Lastly, rs6427528 showed the strongest standing as biomarker before our study. It was associated with response to etanercept with a p-value that approached the GWAS cut-off [14]. However, not all the evidence in the original study was favorable to the association, and the meta-analysis of all results has shown a very weakened association.
The lack of reproducibility we have observed adds to the circumspection advocated in other analyses and reviews [6,23]. Currently, even the most promising biomarkers are not without doubts. For example, the first SNP that reached the GWAS level of significance, rs3794271 [34], was not replicated in a subsequent and larger study [35]. Other example is the mentioned rs6427528 in CD84 [14] that has been very undermined in our meta-analysis. However, there are reasons to expect that some SNPs will be confirmed in future studies. They could be among the SNPs that have demonstrated significant association in a meta-analysis [6], or the associations below the 5 x10 -8 threshold from two recent GWAS [18,19].
The lack of reproducibility cannot be attributed to a weak genetic component in the response to TNFi phenotype. On the contrary, the heritability of ΔDAS28 has been estimated at 0.59-0.71, which similar to other complex traits [36]. Nonetheless, the genetic structure of this component seems very complex [23,37], resembling other common diseases and biological traits as exemplified by the more than one hundred RA susceptibility loci [38]. This polygenic causality means that most loci have a small effect, requiring larger studies for their identification than the available for response to TNFi.
Other aspects that could contribute to the lack of reproducibility include variability of the response along the time and in function of the drug. The effect of time is well-known in clinical practice, where primary and secondary failure to TNFi are distinguished [39]. In this regard, there are some SNPs showing association with the response only at a given time of assessment [6], as the observed in the current study with NUBL. The effect of a specific drug, in turn, could be justified by the known differences between the TNFi [25,26]. Its potential impact on the identification of biomarkers has been revealed by the higher heritability of the response to the monoclonal antibodies (infliximab and adalimumab) than to etanercept (a soluble receptor) [23] and several genetic studies leading to the identification of drug-specific associations [14,19,21,22]. The specific association of the NUBPL SNP with the response to etanercept we have found contributes to the accumulating evidence in this direction.
Another factor likely contributing to the lack of reproducibility is the heterogeneity of the outcomes [36,40]. This issue has been poorly studied, but there are differences in heritability between the components of the DAS28 [36]. The highest estimates of heritability were obtained for the change in the swollen joint count (ΔSJC) and in the tender joint count (ΔTJC). They were followed by ΔDAS28, while heritability of change in the erythrocyte sedimentation rate (ΔESR) and in the general health scale (ΔVAS) were negligible [36]. Therefore, future studies could take advantage of this information focusing on the high heritability components. There are also reasons to prefer ΔDAS28 over the EULAR response. As a general rule, the dichotomized variables have many disadvantages in comparison with the continuous ones [41,42]. In consequence, ΔDAS28 should be more informative and less prone to false positive results than the EULAR response, but no study has compared their performance for the search of biomarkers. Finally, a large within-patient variation at successive visits has been described [40]. This variability decreased the apparent heritability of the response and reduced the power to identify genetic biomarkers. Consequently, the use of the mean of the outcomes assessed at short time intervals has been recommended [40].
Our study serves as an invitation to revise the search for genetic biomarkers of response to TNFi given the difficulties encountered for their validations. It also serves to highlight the complexity of the response phenotype, which likely involves many loci with small effects, and heterogeneity in function of the specific drug and time of assessment. Therefore, identification of reproducible biomarkers will require larger collections of samples and more aggressive patient stratification. In addition, this identification would benefit from the incorporation of alternative and frequent measures of response.
Supporting information S1