Replication study of polymorphisms associated with response to methotrexate in patients with rheumatoid arthritis

About 70 genetic studies have already addressed the need of biomarkers to predict the response of patients with rheumatoid arthritis (RA) to methotrexate (MTX) treatment. However, no genetic biomarker has yet been sufficiently validated. Here, we aimed to replicate a selection of 25 SNPs in the largest collection of patients up to date, which consisted of 915 patients treated with MTX. The change in disease activity (measured as ΔDAS28) from baseline was considered the primary outcome. In addition, response according to widely used criteria (EULAR) was taken as secondary outcome. We considered consistency between outcomes, P values accounting for the number of SNPs, and independence from potential confounders for interpretation of the results. Only the rs1801394 SNP in MTRR fulfilled the high association standards. Its minor allele was associated with less improvement than the major allele according to ΔDAS28 (p = 0.0016), and EULAR response (p = 0.004), with independence of sex, age, baseline DAS28, smoking, seropositivity, concomitant corticosteroid use or previous treatments. In addition, previous evidence suggests the association of this SNP with response to MTX in another autoimmune disease, juvenile idiopathic arthritis, and with high intracellular folate levels, which could contribute to poor response.

70 studies have explored polymorphisms in candidate genes including the primary target of MTX, the enzyme dihydrofolate reductase (DHFR), which converts dihydrofolate to tetrahydrofolate. Tetrahydrofolate is essential for purine synthesis, and in the biologically active form, 5-methyl-tetrahydrofolate, is an important cofactor in one-carbon metabolism. Candidate genes are involved in the transport of MTX, generation of its active metabolites, MTX polyglutamates, or in folate one-carbon metabolism, purine synthesis, or in the biosynthetic pathway for adenosine whose accumulation could be the main anti-inflammatory effect of MTX treatment 3,4 . However, no reproducible association has emerged yet from these candidate gene studies. This panorama has been completed recently, thanks to the first GWAS, which was performed in RA patients from the North of India 6 . The authors highlighted 7 SNPs with suggestive evidence of association, none of them located in candidate loci, and several SNPs in candidate genes, which were only apparent after specifically looking for them.

Materials and Methods
Patients. The study population consisted of 956 patients with RA according to the 1987 American College of Rheumatology criteria 23 . They were of European Caucasian origin, and more specifically of Spanish (269), Slovenian (103), Portuguese (230) and Polish (354) origin. The patients were recruited with written informed consent. The study was approved by the Autonomous Research Ethics Committee of Galicia (Ref. 2014/387). All protocols and methods were conducted according with the relevant guidelines (Declaration of Helsinki and the Belmont Report) and regulations (Spanish Law 14/2007 of Biomedical Research). All patients have received MTX monotherapy for 6 months. Evaluations included Disease Activity Score 28 (DAS28) at the start of treatment and at 6 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 24 . DAS28, at the start of the treatment and at 6 months, sex and age were used for all analyses. This information was incomplete in 41 patients that were excluded. Main clinical characteristics of the remaining 915 patients are summarized in Table 1.
SNP selection and testing. We selected 25 SNP for study (Supplementary Table S1 for full details and references): 14 SNPs from candidate gene studies, and 11 SNPS from a GWAS. Genotypes were determined by PCR amplification followed by single-base extension with the SNaPshot Multiplex Kit (Applied Biosystems, Foster City, California). The corresponding primers, probes and detailed protocols are available from the authors. Quality control of the results included duplicate assays of 10% of the samples, genotype call rate ≥95%, concordance with the Hardy-Weinberg equilibrium (p > 0.01) and with the SNP frequencies in the HapMap European collection (Utah residents with Northern and Western European ancestry and Iberian populations in Spain taken from http://browser.1000genomes.org/index.html). One SNP (rs1051266) was excluded due to low call rate. Here, good and moderate responders were grouped as responders (n = 604) to confront with the NR (n = 311). Linear and logistic regression models for ΔDAS28 and EULAR NR 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. Effect of each SNP was adjusted for sex, age, DAS28 at baseline, treatment with a previous DMARD and concomitant treatment with corticosteroids. In a second analysis, additional putative confounders were included: seropositivity and smoking status. Data from the four populations were combined by individual participant data metaanalysis applying two-stage methods. They involved estimation of the association parameters for each population as the first stage and the combination of the results of the four populations in the second stage. A main effects linear regression model for each population was applied to ΔDAS28. The EULAR outcome, in turn, was assessed with separate logistic regression for each population. The combination was done as fixed effects, weighting the effect corresponding to each population with the inverse variance method, and as random effects, according to DerSimonian and Laird. The fixed effects results were preferred except when heterogeneity was moderate (I 2 > 50), in which case the random effects metaanalysis was selected. Heterogeneity of effect sizes between populations was assessed with Cochran Q test and the inconsistency parameter I 2 . Results were interpreted considering the number of tests following the Bonferroni approach, concordance of results with the two treatment outcomes, independence from confounding variables, and previous evidence of association with response to MTX. Statistica 7.0 (Statsoft, Tulsa, OK, USA) was used for analyses. Post-hoc power analysis was used to assess the study capacity to replicate previous findings. The effect size to replicate was taken from previous studies, considering the EULAR response because this is the most frequently reported outcome in previous studies. The power analysis was conducted using the GPower NT software for the most extreme reported OR 26 .
Data availability. The datasets generated during the current study are available from the corresponding author on reasonable request

Results
The 915 RA patients were predominantly female (81%) and DMARD-naïve (80%) with a mean age of 53 years at the start of MTX monotherapy (Table 1). According to baseline DAS28, many patients presented a high disease activity (DAS28 > 5.1) and most of them received concomitant treatment with corticosteroids or NSAIDs (72.1% and 53.5%, respectively). After 6 months on MTX, only 27.3% of the patients showed a good response according to EULAR criteria, 38.7% showed a moderate response and 34.0% were NR. Several of the clinical and demographic variables were significantly associated with ΔDAS28 or with NR (Table 1). Specifically, men, seronegative patients, patients lacking erosions, not DMARD-experienced or not receiving concomitant corticosteroids showed larger decrease in DAS28 than the corresponding alternative patient subgroups. Associations with NR were similar (Table 1).
After genotyping and quality control, one of the 25 SNPs was discarded. Fixed effect metaanalysis of the nested models on the remaining 24 SNPs showed nominal association with ΔDAS28 (Table 2) of two SNPs from candidate genes (rs2295553 and rs1801394). Only rs1801394 in MTRR showed association withstanding Bonferroni correction (P = 0.0016). Its minor allele (A) was associated with a lower decrease in DAS28 after MTX treatment than the major allele, specifically, 0.14 less DAS28 units per each A allele according with the additive model. The magnitude of change showed a low level of heterogeneity between the four populations (I 2 = 24%) that was not significant according with the Q test (P = 0.3) and did not reflect any widely divergent result (Fig. 1).
The previous results were obtained with adjustment by five clinical variables (sex, age, baseline DAS28, previous DMARD treatment and concomitant corticosteroids), which were available from most of the 915 patients. Multivariate analysis adjusting for two other variables associated with response to MTX (smoking and seropositivity) was only possible in 643 patients (Supplementary Table S2). Therefore, the results were compared at the level of effect sizes, which are independent of the number of samples, not at the level of P values, which are very sensitive to the number of patients. The effect size of rs1801394 was almost unaltered (from β = −0.14 to −0.13; P = 0.009) after the complete adjusted analysis. No other SNP than rs1801394 showed association (P < 0.05) in this set of analyses.
Consideration of the secondary outcome of response to MTX, NR according to EULAR criteria, showed a consistent association of rs1801394 (Table 3). The A allele showed OR = 1.39 for NR (95% CI 1.11-1.75) without detectable heterogeneity between the populations (I 2 = 0). There were three other SNPs showing nominal association (rs4673990, rs2372536 and rs2650972) with little heterogeneity between the populations. As for the primary outcome, inclusion of more variables of adjustment did not modify the effect size of rs1801394 (from OR = 1.39 to 1.37) reinforcing the evidence of its independent association with MTX response (Supplementary Table S2). Two of the other SNPs showing association with NR were reinforced in the fully adjusted analysis (rs4673990 and rs2372536), but at the expense of increased heterogeneity between the populations.

Discussion
Here, we have identified a SNP, rs1801394, in MTRR associated with response to MTX with all the characteristics we had specified in advance to increase reproducibility 22 . This focus in reproducibility differs from many previous studies on the same field [3][4][5]    The association of rs1801394 (also called A66G) with response to MTX fulfilled all the requirements of our analysis: its P value was low enough to allow for correction by the number of tests, it was consistently observed with the two treatment outcomes (although the NR P value was over the corrected threshold), ΔDAS28 and EULAR criteria, and it was independent of the other variables with influence on response to MTX, including differences in disease activity, treatment and sample collection. These requisites were considered necessary to increase reproducibility, which is critical for validation of any biomarker [3][4][5]22 .
The uncritical use of P < 0.05 for claiming association has been shown to lead to many false positive results. Correction by the number of tests, although not without limitations, is widely considered to be very beneficial in genetics of complex diseases 22,27 . Other characteristic that distinguishes our analysis is the adjustment for multiple potential confounding variables. Variables like the baseline level of disease activity or the previous use of DMARDs could have a different distribution in patients with different genotypes and introduce bias in the results as they are strongly associated with MTX response. Another distinctive characteristic of our study is analysis according to the additive genetic model with values 0, 1 and 2 for the number of minor alleles. In selecting the additive model we have followed the approach of the GWAS that have chosen the additive model based in its biological plausibility and statistical reproducibility 27 .
Less consensus exists in the outcome that should be used in the search of biomarkers of response to treatment in RA. ΔDAS28 has all the advantages of quantitative traits in comparison with dichotomous traits for statistics: discrimination, robustness, increased power and reproducibility 28,29 . These considerations motivated our use of ΔDAS28 as the primary outcome for analysis, which is concordant with the large studies and joint initiatives analyzing response to biologic DMARDs [30][31][32] . The only advantage of dichotomous outcomes is its simplicity. However, this simplicity is misleading because definitions of NR vary widely between studies. In effect, some authors have used the EULAR criteria of response, others have used DAS28-based remission, or low-disease activity, or the ACR 20 and ACR 50 response criteria, and even study-specific definitions of response 4 Table 3. Metaanalysis of the association of the 24 SNPs with response to MTX according to the EULAR criteria at 6 months in the four populations of RA patients. Summary odds ratios (OR) and their confidence intervals (95% CI) obtained with fixed effects and random effects metaanalysis are presented together with the inconsistency (I 2 ) and p value of the Cochran Q test for heterogeneity. Multivariate logistic regression in each population included sex, age, baseline DAS28, previous use of DMARD and concomitant treatment with corticosteroids.
SCIENTIFIC REPORTS | (2018) 8:7342 | DOI:10.1038/s41598-018-25634-y consequence, it is suspected that heterogeneity in these definitions could contribute to lack of reproducibility. In any case, we decided to include the EULAR criteria as a secondary outcome for completeness. All the discussed characteristics were intended to increase the reproducibility of our results. However, the association of response to MTX with rs1801394 in MTRR is not yet certain. Additional support could be obtained from its credibility as a good candidate 22 . This SNPs causes the substitution of isoleucine for methionine at amino acid 22 of MTRR, 5-methyltetrahydrofolate-homocysteine methyltransferase reductase. This enzyme collaborates with MTR (5-methyltetrahydrofolate-homocysteine methyltransferase) in the transfer of a methyl group from 5-methy tetra-hydro-folic acid to homocysteine, resulting in de novo production of L-methionine and tetra-hydro-folic acid (in a reaction involving vitamin B12 as cofactor). This step in the folic acid cycle precedes the production of the active form of L-methionine, S-adenosyl-methionine, which is the methyl donor in hundreds of biologic transmethylation reactions and the donor of propylamine in polyamine synthesis. Specifically, the A allele, herein associated with poor response to MTX, has been related to a moderate decrease in total plasma levels of homocysteine, which is interpreted as increased capacity to regenerate methionine 33 . More pertinent for our study, patients with RA under MTX treatment displayed folate levels in red blood cells (RBC) that were higher in carriers of the rs1801394 A allele compared to the G/G genotype 16 . This relationship between MTRR and intracellular folate could be linked with the association we have observed because RA patients displaying a decrease in RBC folate (in the form of folate polyglutamates) concentration during MTX treatment showed a marked decrease of DAS28 at 6 months 18 . This indirect relationship seems to indicate that the A allele could decrease response to MTX by leading to increased intracellular folate, which is a factor already suspected due to the interference of folic or folinic acid supplementation with MTX treatment [34][35][36] . However, the more direct evidence linking rs1801394 with response to MTX comes from previous genetic studies. There have been three studies analyzing rs1801394 as a biomarker of response to MTX in patients with RA 18,19,37 , and two in patients with juvenile idiopathic arthritis (JIA) 38,39 . Two of these studies, one in RA 19 , and one in JIA 38 have found the same association reported here. The negative studies were either smaller (n = 48 in 18 , and 119 in 39 ) or from a different ethnic background (n = 322 Indians in 37 ). As the RA study that has found association overlaps with our study 19 , we analyzed the association between rs1801394 and treatment response after excluding the overlapping collection of patients, finding similar effect sizes (n = 664 patients; ΔDAS28 β = −0.19, SE = 0.06, P = 0.002; NR OR = 1.33, 95% CI = 1.03-1.7, P = 0.030) to those obtained in the total population. Therefore, we can conclude that previous evidence supports a functional relationship of rs1801394 with alterations in the folate pathway and provide suggestive support of its association with response to treatment, although further validation is still needed.
In this study, we analyzed other 23 SNPs from previous candidate gene studies and from a GWAS (Supplementary Table S1). However, none showed an association below the corrected threshold or was associated with the two outcomes at the nominal level of P < 0.05. In addition, three of the four that showed nominal association in any analysis, rs2295553, rs4673990 and rs2372536, were in the opposite direction to previous reports. The only concordant SNP, rs2650972, was just below the 0.05 level for the EULAR outcome (P = 0.035) and far over this threshold for ΔDAS28. This lack of reproducibility is common for candidate gene studies, as signaled by previous systematic revisions and meta-analysis [3][4][5] . The causes are many, but it is worth to highlight that the candidate gene approach responds to a concept of complex traits that has been demonstrated overly optimistic. In fact, progress in genetics has shown that many complex traits have very polygenic causality with many involved loci, each of them with very small contribution. Therefore, the candidate gene approach has been almost abandoned in favor of GWAS, which are done in much larger sample collections and with strict quality control and requirements to claim association 27 . Unfortunately, the new requirements are hard to fulfill in areas like the response to treatment in RA. This is exemplified by the study of response to TNF inhibitors, where only after multiple consortia it was possible to put together about 3000 patients 30 . The GWAS on response to MTX in comparison only included 457 patients and none of the highlighted loci reached the threshold of significance set for this type of studies 6 . Therefore, it is not surprising that we did not replicate any of the SNPs selected from that GWAS in our study.
Our study shows some limitations that include lack of complete information from all patients, possible confounding by indication common to retrospective studies without randomization, the requirement of 6 months on MTX and insufficient power to detect associations with small effect sizes, which are common in complex traits. However, statistical power of our study was sufficient to replicate previous published effect sizes. Only a SNP did not show power >0.8 to replicate previously findings with P < 0.05. This SNP was rs162040, which showed power = 0.71 for detecting association with P < 0.05. As the OR of rs162040 in our study was 1.35 without heterogeneity and the OR in the previous study was 1.45 and with similar increase in NR associated with the minor allele 6 , it is possible that lack of power has been an issue in our results for this SNP. This possibility is interesting because rs162040 is also in the MTRR locus. The next weakest power was obtained for rs2847153 that showed power = 0.82 for association with P < 0.05. However, in this case the effect size in our samples was so low (OR = 1.02) that lack of power did not seem relevant. All the other SNPs showed power ≥0.79 for replication of published effects with P < 0.01. Therefore, the association of rs1801394 needs to be replicated in other context not only to confirm its accuracy, but also to know if it extends to other patient groups, including other ethnic backgrounds, and study designs.
In conclusion, our results show only a SNP, rs1801394, in MTRR with a robust association with response to MTX monotherapy among the 24 analyzed. This association surpassed all the requirements we have imposed aiming to increase reproducibility and was supported by the study sample size, the largest to date on MTX treatment, and suggestive evidence about the SNP function and previous association data. However, it still needs confirmation. It is to be expected this finding could stimulate further collaborations and efforts to progress in the search of biomarkers predicting the important outcome of MTX treatment.