Validation study of genetic biomarkers of response to TNF inhibitors in rheumatoid arthritis

Genetic biomarkers are sought to personalize treatment of patients with rheumatoid arthritis (RA), given their variable response to TNF inhibitors (TNFi). However, no genetic biomaker is yet sufficiently validated. Here, we report a validation study of 18 previously reported genetic biomarkers, including 11 from GWAS of response to TNFi. The validation was attempted in 581 patients with RA that had not been treated with biologic antirheumatic drugs previously. Their response to TNFi was evaluated at 3, 6 and 12 months in two ways: change in the DAS28 measure of disease activity, and according to the EULAR criteria for response to antirheumatic drugs. Association of these parameters with the genotypes, obtained by PCR amplification followed by single-base extension, was tested with regression analysis. These analyses were adjusted for baseline DAS28, sex, and the specific TNFi. However, none of the proposed biomarkers was validated, as none showed association with response to TNFi in our study, even at the time of assessment and with the outcome that showed the most significant result in previous studies. These negative results are notable because this was the first independent validation study for 12 of the biomarkers, and because they indicate that prudence is needed in the interpretation of the proposed biomarkers of response to TNFi even when they are supported by very low p values. The results also emphasize the requirement of independent replication for validation, and the need to search protocols that could increase reproducibility of the biomarkers of response to TNFi.

Introduction Rheumatoid arthritis (RA) is a systemic autoimmune disease mainly characterized by inflammation of synovial joints [1]. If poorly treated, RA is very painful and incapacitating, and it can lead to joint deformities, loss of job and mobility, and premature death. Currently, the prognosis is much better than before the turn of the XXI century thanks to drugs that are specifically targeted to immune mediators [1,2]. The first drugs of this class were the TNF inhibitors (TNFi), infliximab, and adalimumab, which are monoclonal anti-TNF antibodies, and etanercept, which is a recombinant soluble TNF receptor. Now, more TNFi are available together with other drugs targeting IL6, B cells, T cells or intracellular kinases. This range of drugs is welcomed because none of them is effective in all patients. Typically about 30% of the patients fail to respond to any of the drugs, and an additional 30% of patients show only a partial response. This between-patient variability has not prevented the rheumatologists to aim for remission or low disease activity, which they seek by changing from one drug to another, and by combining them with conventional antirheumatic drugs [1,2]. The election of drug follows a trial and error approach, which is very unsatisfactory because these drugs are expensive and have potential side effects. Furthermore, delayed control of the disease worsens long term prognosis. A very attractive alternative will be to use biomarkers for personalizing the treatment [3,4].
An important effort has been directed to the discovery of genetic biomarkers of response to TNFi [5,6]. It has involved candidate gene studies and Genome-Wide Association studies (GWAs). This effort has led to some remarkable findings: the two SNPs (rs3794271 and rs284511) that have achieved association with response to TNFi below the GWAS significance threshold of p < 5 x10 -8 [7,8]; the association of the PTPRC locus in three large independent studies [9][10][11]; and two other SNPs (rs6427528 and rs113878252) with very convincing evidence of association with response to etanercept [12,13]. Regrettably, none of these results or any of the other proposed biomarkers is sufficiently validated, either because an independent validation is still pending, or because of lack of replication in other studies.

Patients
The study population consisted of 581 patients with RA according to the ACR 1987 criteria [19]. They were treated with an TNFi as the first targeted antirheumatic drug by indication of the attending rheumatologist. This treatment indication was made with independence of the study. Recruitment of samples was approved by the local ethics committees and patients gave their written informed consent. The study was conducted according to the principles of the Declaration of Helsinki (2013) and was approved by the Research Ethics Committee of Galicia (Spain, code 2014/387). The patients were of European origin (Spanish ancestry = 530 and Greek ancestry = 51). Change in Disease Activity Score 28 joints (ΔDAS28 = DAS28baseline − DAS28follow-up) was used as the primary outcome. The DAS28 composite index includes erythrocyte sedimentation rate, self-reported global patient health, and counts of swollen joints and of tender joints among 28 selected joints [20]. A score over 5.1 indicates high disease activity, whereas below 3.2 indicates low activity. In addition, response according to the European League Against Rheumatism (EULAR) criteria was assessed [21]. Non-responders are defined as patients showing 0.6 improvement from baseline DAS28, or showing modest improvement, 1.2, but remaining with high disease activity, DAS28 > 5.1. Good responders show DAS28 3.2 and improvement from baseline > 1.2. The remaining patients are considered moderate responders.

Genotyping
Eighteen SNPs were selected because of previous association with response to TNFi in RA irrespective of ethnicity or country of origin (Table A in S1 Text) [7,8,[14][15][16][17][18]. Genotyping was conducted by multiplex PCR amplification followed by single-base extension with the SNaPshot Multiplex Kit (Applied Biosystems, Foster City, California). The primers, probes and detailed protocols used for these analyses are available from the authors upon request. Replication in duplicate genotyping reactions of 10% of the samples, genotype call rate ! 94%, concordance with the Hardy-Weinberg equilibrium (p<0.01), and with SNP frequencies in the HapMap collections [22], were used for quality control (Table B in S1 Text). Two SNPs (rs10520789 and rs1679568) did not follow the Hardy-Weinberg equilibrium and were excluded. The R statistical package SNPAssoc was used for Hardy-Weinberg equilibrium assessment and allelic frequency calculation [23].

Statistical analysis
The main outcome of treatment was ΔDAS28 at 6 months. EULAR non-responders compared with responders (good + moderate) or non-responders compared with good responders (moderate responders left out) at any time, and ΔDAS28 at 3 or 12 months were considered secondary outcomes. Linear and logistic regression models for ΔDAS28 and EULAR response were fitted, respectively. Genotypes were considered in accordance with 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. Analyses were adjusted for baseline DAS28, sex and the specific TNFi. Results were interpreted considering the number of tests following the Bonferroni approach, consistence of results with the different outcomes, independence from confounding variables, and previous reports. Statistica 7.0 software (Statsoft, Tulsa, OK, USA) was used for the analyses.
Data availability: all data included in this study is available in the S1 Dataset.

Characteristics of the RA patients
Fifteen patients lacking baseline DAS28 information were excluded. The remaining 566 RA patients showed the demographic and clinical characteristics shown in Table 1. They presented a severe disease before the start of TNFi, as indicated by the frequent presence of erosions and of the specific RA autoantibodies, and by the high mean disease activity (mean DAS28 > 5.1) that was not controlled by previous treatments. The patients were evaluated during the course of treatment with infliximab (n = 381), adalimumab (n = 152) or etanercept (n = 33). Response to treatment was evaluated at 3, 6 or 12 months of follow-up with information available for a variable number of patients at each point (n = 289 at 3 months, 532 at 6 months and 365 at 12 months). All of them except 2.6% received TNFi together with conventional antirheumatic drugs. The treatment was associated with a significant decrease in RA activity, which was more marked during the first 3 months than at later times. The mean DAS28 value remained in a narrow range from the 3 rd month onward (4.1 to 3.7). In a similar way, the frequency of nonresponders according to the EULAR criteria was very uniform at the three assessment times (21.1 to 23.0%).

Analysis of association with response to TNFi
Univariate analyses identified clinical and demographic variables that were associated with ΔDAS28 at 6 months ( Table 2). The strongest association was found between ΔDAS28 and baseline DAS28. It was followed by the specific TNFi used, and by sex. These three variables retained association in the multivariate analysis including them (Table 2). Therefore, they were retained as potential confounders for all analyses. No association between ΔDAS28 at 6 months and age at the start of treatment, time since disease onset, presence of autoantibodies (rheumatoid factor or anti-cyclic citrullinated peptides) or erosions, smoking, concomitant antirheumatic drug, or baseline disability was observed ( Table 2). Two of the 18 SNPs failed quality control and were excluded from analysis. None of the remaining 16 SNPs was associated with response to TNFi in our patients when the threshold of significance was adjusted for the number of tests (p < 0.003). Specifically, none of the SNPs was associated with the primary outcome, ΔDAS28 at 6 months, or with the secondary outcomes ΔDAS28 at 3 or at 12 months, even at the level of the uncorrected 0.05 p value (Table 3). Regarding the secondary outcomes based in EULAR criteria, only a SNP showed nominal association with EULAR non-response (rs10925026 at 6 months), but it was not below the adjusted p value, and the direction of changes was not consistent at other assessment times (Table 4).

Direct comparison with previous studies
We also checked the consistency of our results with the previous studies in which the biomarkers had been proposed [7,8,[14][15][16][17][18]. This implied to consider the same time of assessment, measure of response, reference allele, and direction of change (improvement or lack of improvement) that produced the lowest p value in previous studies (Table 5). In addition, we considered the comparability of the minor allele frequencies (MAF), the proportions of Validation study of response genetic biomarkers patients on each specific TNFi and on concomitant treatment with cDMARD, and the fraction of ACPA positive patients, as all these factors could influence the association. First, we compared the direction of effect. Only 5 of the 12 SNPs, in which this information was available from the previous studies, showed consistency in direction between our results  Abbreviations: OR, odds ratio; CI, confidence interval and the previous reports, one was completely neutral (rs4925648), and the other 6 SNPs showed opposite direction (Table 5). Next, we focused in the effect sizes of the 5 SNPs showing the same direction of effect in the two studies. For the 2 SNPs assessed as ΔDAS28 (rs7070180, rs284511), the regression coefficients were markedly lower in our study than in the previous ones. Also, for the 3 SNPs assessed according with EULAR criteria (rs4612666, rs717117, rs3794271), the odds ratios were notably nearer 1.0 in our results than in previous reports. Information of direction or size of effect for the remaining 4 SNPs (rs885813, rs885814, rs11870477, rs16973982) was not available in the previous report [14], and only comparison of the p values was possible. The contrast between p values for these 4 SNPs was very notable, because none of them showed even a trend to association in our results, but showed p values 1 x10 -5 in the previous report [14]. We also found the MAF of most of the SNPs in our patients were similar to the observed in the previous studies reporting association (Table A in S1 Text). All the MAF were less than 20% different from the previously reported in Europeans, except for rs717117 that was about half as frequent in our patients as in the Danish patients from the previous report [14]. The MAF differences were larger (+30% and -33%) for the two SNPs, rs284511 and rs284515, previously reported in Japanese patients [8]. Additionally, our patients were comparable with previous ones in the percentage receiving concomitant treatment with cDMARD, which ranged from 73 to 87% except for a study with all patients on combined treatment [16], and in the fraction of ACPA positive patients, which ranged from 59 to 89% (Table B in S1 Text). Table 5. Direct comparison of results for the same outcomes in the current and previous studies. Outcomes were either ΔDAS28 at 3 or 6 months, or the EULAR criteria, which were used as a comparison of non-responders (NR) with responders (good + moderate responders (GR + M)) at 3 months, or a comparison of NR with GR either at 3 or 6 months. Only previous studies with significant associations are included. Comparison of effect size and direction, and of p values. However, there was a lot of heterogeneity in the distribution of the specific TNFi (Table B in S1 Text). Five studies (three of which were from the same patient registry) showed an even distribution in three fractions receiving either infliximab, adalimumab or etanercept [7,15,17,18,24]. Two studies, including the current one, showed about 70% of patients on infliximab and 25% on adalimumab with only 5% on etanercept [14]. Other two studies were dominated by patients on infliximab or etanercept, with only 10% on adalimumab [8,16]; and finally, a unique study included 65% of patients treated with adalimumab, with the remaining patients distributed between the other two TNFi [25].

Discussion
Our attempt to validate genetic biomarkers of response to TNFi has led to increased doubts on their status. Consequently, our results emphasize the need of independent validation of all genetic findings, even of those that reach the GWAS level of significance. In addition, the results call for a revision of the ideas justifying the search of genetic biomarkers, and to consider possible ways to improve their reproducibility. At a first glance, the most surprising lack of association in our patients corresponds to rs3794271 in PDE3A-SLCO1C1, and rs284511 near MAP3K7, because the two SNPs have reached p < 5 x10 -8 , the GWAS threshold of significance, in previous studies [7,8]. A more detailed consideration, yet, indicates that our results are not so striking. In relation with rs3794271, because its role as biomarker has already been contradicted by lack of association with response to TNFi in 1750 UK patients with RA [24]. This study was much larger than the two studies that discovered this putative biomarker (196 Danish and 315 Spanish patients, respectively) [7,14]. A fraction of the patients in the UK study differed from the patients in the two other studies, but these differences seemed irrelevant for the lack of replication [24]. Our study, which was in all respects comparable to the Danish and Spanish studies [7,14], reinforces this conclusion. In relation to the other SNP, rs284511 (and rs284515 in the same MAP3K7 locus), our study is the first independent validation attempt, but the possibility of an ethnic factor in the discordant results should be reckoned as it was discovered in Japanese patients [8]. An ethnic factor is well known in the genetic component of RA, where most loci are shared between Europeans and Asians, but 5 loci are specific to one or the other ethnic group [26]. This possibility was also suggested by the difference in MAF of the two MAP3K7 SNPs between our study and the Japanese patients, which was larger than the present for SNPs from European reports. Asian specificity of the MAP3K7 association would explain our result and that none of the previous GWAS in Europeans have mentioned these two SNPs [12][13][14]16,25,27].
Other 5 SNPs were initially discovered in the same Danish GWAS that uncovered the already discussed SNP in PDE3A-SLCO1C1 [14]. All of them showed p < 10 −4 in the 196 patients of that study, but none obtained subsequently any support. The low frequency rs717117 SNP has never been reported again. None of the other 4 SNPs was associated with response to TNFi in a GWAS of 882 RA patients from the Netherlands [25], and rs11870477 was not replicated in 315 patients of a Spanish study [7]. This latter study did not analyze the other 4 SNPs.
More intriguing is the status of rs7070180 in the GFRA1 gene, whose investigation has been prone to mishaps. This SNP showed significant association (p < 10 −4 ) with response to TNFi in stages 1 and 2 of a UK study (945 patients in total), but its assay failed in stage 3 [16]. Subsequent studies never mentioned it, until a recent GWAS in Japanese patients [8]. This study identified rs1679568, other SNP in the same GFRA1 locus, as associated with response to TNFi with p < 10 −6 . The authors took this result as concordant with the rs7070180 association in UK patients because rs7070180 was not genotyped or imputed in the Japanese patients [8]. Regrettably, our results were also incomplete. On one side, the rs1679568 assay did not pass quality control in our study. On the other, our results for rs7070180 where inconclusive: we found a smaller effect size than in Plant et al. (β = 0.07 vs. 0.30) [16], but concordant direction of effect and a trend to association (p = 0.07). Consequently, we think the status of this locus is still ambiguous.
The remaining seven SNPs had shown association with response to TNFi that was significant according with the criteria applied to their discovery studies [15,17,18], but that was less convincing from the statistical point of view (p values between 0.05 and 0.003) than the SNPs discussed above. The additional criteria leading to take these SNPs as potential biomarkers included functional studies that established a correlation between the SNPs and molecular changes that seemed important for RA pathogenesis. Six of these SNPs are located in two interesting genes from the point of view of RA pathogenesis, NLPR3 (four SNPs) and CARD8 (two SNPs). Five of them were associated with response to TNFi in a study of 1278 UK patients with RA [15]. In the same study, the authors found differences in expression of the NLPR3 and CARD8 genes between monocytes from RA patients and controls, and association linking the SNPs with the expression of NLPR3 and CARD8 in monocytes [15]. Further support for association of NLPR3 with the response to TNFi was provided by an independent study of 538 Danish patients with RA [17]. In this study, only a NLPR3 SNP (rs4612666) was tested, but it showed association. The final SNP of the seven in this group (rs3468) was selected from a very attractive study. In it, the A allele of rs3468 was identified primarily as the SNP that was most associated (p < 1.1 x10 -6 ) with increased DNA methylation at two cis-CpG sites on the LRPAP1 gene, which, in turn, was strongly associated with response to etanercept (p < 1.7 x10 -8 ). In addition, there was weak genetic association (p < 0.05) of the A allele of rs3468 with non-response to TNFi in two sets of RA patients (with 56 and 1204 patients, respectively) [18]. Therefore, the relationship was established mainly as indirect. Our results were clearly contrary to the previously reported and, together with lack of replication of the other six SNPs in this group, reinforce the well-known fact that the main determinant of reproducibility of genetic associations is the level of statistical association [28,29].
In addition to the specific comments for each SNP of the previous paragraphs, there are other factors that could have contributed to the lack of reproducibility of the results. The possible impact of differences in ethnicity has already been discussed in relation with the MAP3K7 SNPs. All other SNPs were described in studies done with European patients and ethnicity is less likely a factor for them. In addition, MAF were comparable for all the SNPs from European studies except for rs717117, whose low frequency in our patients (MAF = 5.3%) made for a lower statistical power than in the original study (MAF = 11%). Other study characteristics that could impinge in the results, the fraction of patients on concomitant treatment with cDMARD and the fraction with ACPA, can also be dismissed as differential because they were comparable between studies. Nevertheless, there was marked variability in the distribution of the specific TNFi that the patients had received. No common pattern was found among the reports, with our patients being only comparable with a previous study. The two sets of patients showed predominant treatment with infliximab and few patients on etanercept. The other studies were characterized by a variety of TNFi distributions. Considering that some patients who are resistant to a TNFi respond to a different one, it is possible that these differences in representation of the various TNFi have contributed to the lack of replication.
Our results invite to reconsider the search of genetic biomarkers of response to TNFi. The first point to ponder is the heritability of the phenotype, which was unknown until recently, when a familial component in the response to TNFi was demonstrated [30], and the heritability of ΔDAS28 was quantified (h 2 from 0.59 to 0.71 in different estimations) as sufficient for genetic studies [31]. However, change in the number of swollen joints or in the number of tender joints, which are components of DAS28, showed larger heritability than ΔDAS28 suggesting that they could be used as more discriminant outcomes of treatment [31]. Possible benefits of these alternative outcomes remain still untested, and we could not assess their performance in our study because they are not available in many of our patients. Other improvement that has been proposed consist in reducing within patient variability by using as outcome the averages of frequently repeated response assessments [32]. Another point of interest is the structure of the genetic component of response to TNFi, which has been clarified by recent studies [27,33,34]. This component is not significantly shared with the RA susceptibility loci, either taken individually, or in combination [9,33]. In addition, rare variants with strong effect do not contribute significantly to the response to TNFi [34]. Therefore, only common variants unrelated with RA susceptibility seem likely, but the effect of each of them appears to be small. This is the conclusion of a thorough study, with many analysis tools and by multiple teams, done in 2706 RA patients [27]. This conclusion emphasizes the need to continue searching with more GWAS, which are very powerful for identifying common genetic variants but need to be large to find those of small effect. Also, it emphasizes the importance of pursuing the best supported biomarkers to establish their value, as PTPRC [9][10][11], a recently discovered SNP associated with response to etanercept [13], the SNPs associated in the Japanese GWAS [8], or other SNPs with significant association in recent meta-analysis [5]. In addition, it is to be expected that biomarkers will be searched in other ethnicities as we lack information about many populations that could show particular biomarkers of value for the RA patients from the corresponding genetic background.
The available evidence, therefore, shows the need for prudence in claiming new biomarkers, the requirement for independent replication of all previous findings, and the convenience of exploring the possibility of drug-specific biomarkers among the TNFi, and the use of outcomes with high heritability, as change in the number of swollen or tender joints, and of using averages of repeated assessment of these outcomes, to increase the success of genetic studies. Progress along these lines will require new and large collections of patients.
Supporting information S1 Text. Information on selected previous studies and compatibility with the current study. It includes two tables. Table A: List of SNPs selected for this study with references and the corresponding quality control results in the current study. Table B: Characteristics of the patients included in previous studies compared with the current study. (DOCX) S1 Dataset. Raw data with all the variables and patients considered in this study. (XLSX)