Pharmacogenetic Sex-Specific Effects of Methotrexate Response in Patients with Rheumatoid Arthritis

Methotrexate (MTX) is a commonly used drug for the treatment of rheumatoid arthritis (RA), but its effectiveness can vary greatly among patients. Pharmacogenetics, the study of how genetic variations can affect drug response, has the potential to improve the personalized treatment of RA by identifying genetic markers that can predict a patient’s response to MTX. However, the field of MTX pharmacogenetics is still in its early stages and there is a lack of consistency among studies. This study aimed to identify genetic markers associated with MTX efficacy and toxicity in a large sample of RA patients, and to investigate the role of clinical covariates and sex-specific effects. Our results have identified an association of ITPA rs1127354 and ABCB1 rs1045642 with response to MTX, polymorphisms of FPGS rs1544105, GGH rs1800909, and MTHFR genes with disease remission, GGH rs1800909 and MTHFR rs1801131 polymorphisms with all adverse events, and ADA rs244076 and MTHFR rs1801131 and rs1801133, However, clinical covariates were more important factors to consider when building predictive models. These findings highlight the potential of pharmacogenetics to improve personalized treatment of RA, but also emphasize the need for further research to fully understand the complex mechanisms involved.


Introduction
Rheumatoid arthritis (RA) is a chronic autoimmune disease that causes inflammation of the synovial tissue, leading to destruction of bones and cartilage [1]. It affects between 0.5% and 1.0% of the population, with a higher prevalence in women and urban areas in the northern hemisphere [2]. The European League Against Rheumatism (EULAR) has been using, since the 1980s, methotrexate (MTX), a conventional synthetic Disease-modifying antirheumatic drug (DMARD), as the first-line drug for RA treatment. However, MTX has a significant inter-patient variability, with 30-50% of patients failing to achieve remission, and up to 30% experiencing adverse effects that exclude them from treatment due to toxicity [3,4]. MTX is also used in the treatment of other autoimmune diseases (such as psoriasis, 10 to 25 mg/week) and certain types of cancer (such as lymphoblastic leukemia, 200 mg/m 2 intravenous over two hours up to high doses of 1-3 g/m 2 as a 24 h continuous infusion), usually a somewhat higher dose than for RA (7.5 to 20 mg/week) [5].
MTX has been studied to explain the variability in clinical response to the drug [4]. To date, over 120 SNPs in 34 genes have been identified as potentially contributing to

Participant Recruitment
For this study, we used the Merida cohort, which consists of RA patients recruited since March 1990 and followed up in monographic medical consultations followed between 2012 and 2015 at the Hospital General de Mérida. These patients were recruited into a prospective observational cohort study using the following inclusion and exclusion criteria. Inclusion criteria: patients over the age of 18 years who were followed between 2012 and 2015, diagnosed with RA according to American College Rheumatology criteria [28], protocolized follow-up in the RA monographic consultation (this information is on the patient's evaluation sheet), receiving, or has ever received, treatment with MTX monotherapy. During monotherapy, patients could also take concomitant nonsteroidal anti-inflammatory drugs, oral prednisone, and folic acid, but not other DMARDs from baseline. Exclusion criteria: high probability of patient loss, inability to provide accurate data, and patients enrolled in clinical trials. Furthermore, patients who had received biologics, cyclophosphamide, or a combination of MTX with other DMARDs prior to monotherapy were excluded.
Informed consent was obtained from all RA patients, and the study was approved by the Bioethics and Biosafety Commission of the University of Extremadura (No. 42/2011).
MTX was administered under a rapid escalation pattern, initially administered orally in a single dose up to 10 mg/day, in two doses starting at 10 mg/day. Patients were switched to the subcutaneous route if they did not tolerate oral administration well, the amount of 15 mg per day was exceeded, or the desired efficacy was not achieved.

DNA Extraction and Genotyping
Eight polymorphisms in seven genes encoding proteins involved in MTX metabolism (pharmacokinetics and pharmacodynamics) and MTX toxicity were selected for genotyping in our patients according to a literature review ( Figure 1). The analyzed polymorphisms were:

Variables of Interest
Two indicators of MTX efficacy were obtained based on the Disease Activity Score (DAS) 28 with C-reactive protein (DAS28-CRP) levels: response to MTX (DAS28 < 3.2) [29] and RA remission (DAS28 < 2.6). The DAS is a scoring system to assess the activity of rheumatoid disease, having been recommended by the EULAR for this purpose both in clinical studies and in daily clinical practice. The DAS index for measuring disease activity was developed during the 1980s [30] and was published by Van der Heijde et al. (1990) [31]. DAS uses 44 joints in its count. The DAS index combines information regarding the number of swollen, tender, and acute phase reactant joints, and a global measure of health status. DAS28 is a simplified version using only 28 joints for counting and was published by Prevoo et al. (1995) [32]. Therefore, DAS28-CRP is a variant of DAS28 that substitutes the erythrocyte sedimentation rate for the C-reactive protein count and takes 28 specific joints.
To establish response criteria, patients had to have received MTX monotherapy at the highest tolerated dose for at least 6 months and responders who achieved and maintained a DAS28-CRP < 2.6 were in clinical remission, regardless of whether they continued MTX monotherapy or no treatment.
Forty-four cofactors and comorbidities were obtained to complete the statistical model ( Table 2). Demographic and clinical variables with more than 5% of cases were added to the analysis. Included confounding variables can be organized into different groups: Demographic variables (n = 8): age, sex at birth, educational level, dedication, consumption of alcohol, tobacco, coffee, and tea.
Pharmacological variables (n = 10): use of DMARDs previously, number of previous DMARDs, time between diagnosis and treatment with DMARD, time between first symptoms and treatment with DMARD, initial doses of prednisone, doses of folic acid, maximum doses of MTX, maximum tolerated doses of MTX, time between diagnosis and treatment with MTX, MTX monotherapy duration.

Statistical Analysis
Multivariable logistic models were used to assess the association between genetic variants and the dependent variables: response to MTX RA remission and adverse effects, as dichotomous variables. The 44 confounding variables gathered in this study were added to each of the linear models. In this first approach, we were not interested in obtaining the best explanatory model but rather in removing any possible confounding effect these variables may have on the dependent variables. For each genotyped SNP, we tested four different models of inheritance: dominant, recessive, codominant, and overdominant.
Once the different associations between genetic variants and the dependent variables were assessed, we used multivariable lasso logistic regression and receiver operating characteristic (ROC) curves to weight the contribution of the significant genetic variants in a prediction model. Lasso logistic regression was developed by randomly reallocating individuals into two groups, training, and test, with equal sizes. K-fold cross-validation was used with 150 folds.
All statistical analyses were performed using R v4.2. Significance was assessed using the False Discovery Rate, and the STROBE (strengthening the reporting of observational studies in epidemiology) methodology was followed [33].

Clinical Characteristics
The aim of this study was to assess the effect of different genetic polymorphisms on the response to MTX in a group of RA patients treated with this drug. Table 2 shows details about the clinical characteristics of cases and controls.
The cohort was composed of 301 patients (67.1% women) with a median age of 49 years. Figures 2-4 show the distribution of response to MTX, RA remission, and the different adverse events for each SNP, respectively.
The aim of this study was to assess the effect of different genetic polymorphisms on the response to MTX in a group of RA patients treated with this drug. Table 2 shows details about the clinical characteristics of cases and controls.
The cohort was composed of 301 patients (67.1% women) with a median age of 49 years. Figures 2-4 show the distribution of response to MTX, RA remission, and the different adverse events for each SNP, respectively. Due to the lack of cases of hematological adverse events and hepatic adverse events (4.6% and 8.9% of cases, respectively; Table 3), the analysis was developed for the total number of adverse events (50.5%) and for gastrointestinal adverse events (28.9%). This is a retrospective study with no previous sample size calculation. Due to the lack of cases of hematological adverse events and hepatic adverse events (4.6% and 8.9% of cases, respectively; Table 3), the analysis was developed for the total number of adverse events (50.5%) and for gastrointestinal adverse events (28.9%). This is a retrospective study with no previous sample size calculation.  Figure 1 shows the allele frequencies and presents a relative comparison between genotypes in a pie chart. In general, the study has a good representation of both alleles for each SNP, but this is not the case for ITPA and AMPD1, where very few alternative alleles are present (9.4% and 11.5%, respectively), and in most cases, they are in heterozygous form. As a result, we have a small representation of homozygous genotypes for these SNPs (seven individuals for ITPA and two individuals for AMPD1), which prevents the testing of certain models of inheritance.

Association between Genetic Variants and MTX Efficacy
The efficacy of the MTX treatment was measured through two variables: response to MTX (DAS28 < 3.2) and AR remission (DAS28 < 2.6).
In Figure 3, forest plots of the different analyses are shown, all the 44 confounding variables were added to the multivariable linear methods. This figure shows that ITPA rs1127354 homozygous C/C genotype is associated with having a higher probability of being a responder to MTX (logOR = 2.7; 95% CI = 0.8, 4.6) ( Figure 3). Thus, chances of being a responder decrease when comparing codominant inheritance C/A vs. C/C or overdominant C/A vs. A/A + C/C and increase in a dominant model (C/C vs. C/A + A/A). ITPA effect was not significant when AR remission was considered.  FPGS rs1544105 and MTHFR rs1801131 polymorphisms showed association with AR remission. Dominant inheritance of FPGS homozygous G/G genotype is also associated with a higher probability of AR remission under a codominant (logOR = −1.9; 95% CI = −3.2, −0.6) and dominant models (logOR = 1.3; 95% CI = 0.1, 2.4). In contrast, FPGS heterozygote G/A genotype decreases the probability of AR remission under an overdominant inheritance (OR = −1.9; 95% CI = −3.1-0.75). Finally, MTHFR rs1801131 T/T homozygous genotype is associated with a higher probability of AR remission in a recessive (logOR = 1.6; 95% CI = 0.3, 3.0) or codominant inheritance (logOR = 1.8; 95% CI = 0.3, 3.2) (Figure 3).
In order to analyze the sex-specific pharmacogenetic effects, the interaction between sex at birth and each genotype was assessed. For those SNPs showing a significant interaction effect with sex at birth, a separate analysis by sex was developed. In Figure 4, panels a and b, it is possible to find some exploratory figures with GGH, FPGS and ABCB1 genotypes, response to MTX and AR remission distributions by sex, while panels c and d show the sex-specific effects in response to MTX and AR remission, for different genes and SNPs ( Figure 4).  Regarding response to MTX (Figure 4c), the effects were male specific, no significant effects were found in women. Thus, GGH rs1800909 heterozygous T/C was found to be associated with an increase prevalence of responders to MTX under an overdominant (logOR = 3.3; 95% CI = 0.1, 6.5) and codominant (logOR = 4.6; 95% CI = 0.1, 9.1) inheritance . Sex-specific effects in AR remission (DAS28 < 2.6) and response to MTX (DAS28 < 3.2). (a) Bar plots with the percentage of responders to MTX (DAS28 < 3.2) in each SNP genotype for GGH rs1800909 and ABCB1 rs1045642. (b) Bar plots with the percentage of individuals with AR remission (DAS28 < 2.6) in each SNP genotype for GGH rs1800909 and FPGS rs1544105. (c) GGH rs1800909 and ABCB1 rs17602729 effect estimates (logOdds ratio) over response to MTX, and 95% CI, for men (x axis) and women (y axis) are shown. (d) GGH rs1800909 and FPGS rs1544105 effect estimates (logOdds ratio) over response to MTX, and 95% CI, for men (x axis) and women (y axis) are shown.

Association between Genetic Variants and MTX Adverse Events
Due to our sample size and the number of adverse events, analyses were only possible using the total number of adverse events and, specifically, the gastrointestinal adverse events, as can be seen in Figure 5. In Figure 6, forest plots of the different analyses are shown, all 44 confounding variables were added to the multivariable linear methods.

Assessing Genomic Effects to Predict RA Remission and Adverse Events
Besides the estimation of the effect of different genes on RA remission and the presence of adverse events, it is also interesting to assess the contribution of these genes to a predictive model. In order to achieve this, we first chose the most relevant covariables using a lasso logistic regression. Among the 44 covariables used in this study, different sets were flagged as relevant. The lasso approach highlighted ten variables for response to MTX and AR remission: radiological erosions, tuberculosis, cardiovascular risk factors, maximum dose of MTX, duration of the MTX treatment, time between first symptoms and treatment, DAS28 levels at basal, dose of folic acid, MTX monotherapy duration, and being a smoker. Eleven variables were highlighted for the presence of any adverse events: DAS28 levels at basal, being a smoker, age, alcohol intake, obesity, tuberculosis, baldness, hypertransaminasemia, number of risk factors, maximum dose of MTX, and MTX monotherapy duration. Figure 8 shows the ROC curves and AUCs for different logistic models for AR remission and the presence of any adverse event. For each of the three dependent variables, the model with the highest AUC is the one that includes all the 44 covariables and the significant genes found in this study. Genes, by themselves, do not have predictive power (AUC = 0.56, 0.58, and 0.57 for response to MTX, AR remission, and all adverse events, respectively). Even more, it seems that adding significant genes to the lasso-flagged list of covariables does not significantly improve the AUC.
The analysis (all patients and by sex) of the response to MTX, the association with the remission of RA, the associations with the presence of any adverse event, are included in eight different complementary tables (see Supplementary Material).

Assessing Genomic Effects to Predict RA Remission and Adverse Events
Besides the estimation of the effect of different genes on RA remission and the presence of adverse events, it is also interesting to assess the contribution of these genes to a predictive model. In order to achieve this, we first chose the most relevant covariables using a lasso logistic regression. Among the 44 covariables used in this study, different sets were flagged as relevant. The lasso approach highlighted ten variables for response to MTX and AR remission: radiological erosions, tuberculosis, cardiovascular risk factors, maximum dose of MTX, duration of the MTX treatment, time between first symptoms and treatment, DAS28 levels at basal, dose of folic acid, MTX monotherapy duration, and being a smoker. Eleven variables were highlighted for the presence of any adverse events: DAS28 levels at basal, being a smoker, age, alcohol intake, obesity, tuberculosis, baldness, hypertransaminasemia, number of risk factors, maximum dose of MTX, and MTX monotherapy duration. Figure 8 shows the ROC curves and AUCs for different logistic models for AR remission and the presence of any adverse event. For each of the three dependent variables, the model with the highest AUC is the one that includes all the 44 covariables and the significant genes found in this study. Genes, by themselves, do not have predictive power (AUC = 0.56, 0.58, and 0.57 for response to MTX, AR remission, and all adverse events, respectively). Even more, it seems that adding significant genes to the lassoflagged list of covariables does not significantly improve the AUC.  (2) and (3) were: DAS28 levels at basal, being a smoker, age, alcohol intake, obesity, tuberculosis, bald-  (2) and (3) were: DAS28 levels at basal, being a smoker, age, alcohol intake, obesity, tuberculosis, baldness, hypertransaminasemia, number of risk factors, maximum dose of MTX, and MTX monotherapy duration. Selected genes and inheritances for models (1), (3), and (4) were ADA rs244076 (codominant inheritance), FPGS rs1544105 (dominant and codominant inheritances), GGH rs1800909 (recessive inheritance), and MTHFR rs1801131 (dominant, codominant, and overdominant inheritances).

Discussion
In the present study, we examined the impact of eight SNPs in seven genes on the effectiveness and toxicity of MTX treatment. We also included clinical records, which allowed us to consider a larger number of factors in our analysis. This may be the reason why the present results differ from those previously reported. The present results suggest that ITPA (rs1127354) and ABCB1 (rs1045642) may have a pharmacogenetic effect on the efficacy of MTX measured through DAS28; ADA (rs244076) and MTHFR (rs1801133) may affects its toxicity, while GGH (rs1800909), MTHFR (rs1801131), and FPGS (rs1544105) may have an effect in both the efficacy and toxicity of MTX treatment.
ITPA gene encodes the inosine triphosphate pyrophosphatase adenosine pathway enzyme which hydrolyzes inosine triphosphate and deoxyinosine triphosphate to the monophosphate nucleotide ( Figure 1). We showed that homozygote ITPA (rs1127354) C/C is associated with an increased probability of being a responder to MTX. This outcome was not observed in most of the previous works [7,19,20]. Also, the heterozygous (C/A) was found to be associated with non-responsiveness [26].
FPGS homozygous G/G genotype (rs1544105), which encodes the folylpolyglutamate synthase mitochondrial enzyme (Figure 1), has been found by this study to be associated with AR remission in males but not in females. Homozygous A/A and heterozygous A/G genotypes were also associated with an increase in gastrointestinal adverse effects in the whole dataset and in females, respectively. This outcome, in the general population, was also found in the bibliography [21]. Nevertheless, most previously published studies did not find a significant association [7,14,19,24,25,34].
GGH gene encodes for the gamma-glutamyl hydrolase which removes gamma-linked glutamate from MTX. Dominant homozygous T/T SNP rs1800909 was found to be associated with both levels of DAS28 (through the response to MTX and AR remission) and adverse events in men, this genotype being related to a better response to medication and fewer adverse events. Effects of GGH are highly heterogenous in the bibliography, with published papers which do find significant association [29] and others that find no effect whatsoever [7,28,30].
MTHFR gene, which encodes for the methylenetetrahydrifolatereductase enzyme, is one of the most studied genes in the context of MTX treatment. Our study found that both SNPs rs1801133 and rs1801131 are associated with MTX toxicity and rs1801131 AR remission. As can be seen in Table 1, the effects of these SNPs are not clear since it is possible to find original articles and meta-analyses with and without [13,20,24,34] significant effects reported for both efficacy and toxicity.
Lack of consistency was also found when analyzing the effects of ABCB1 (rs1045642), which encodes the cell membrane protein P-glycoprotein. We found a significant effect of the heterozygous rs1045642 (C/T), under an overdominant model, in male responders to MTX. However, is possible to find any kind of result in the bibliography as can be observed in Table 1. Finally, heterogeneity of outcomes could also be found in ADA (rs244076) when comparing our outcomes with the ones at Table 1. Heterozygous G/A seems to provide the individual with protection against adverse events (Figures 6 and 7); however, this is the first time ADA (rs244076) was found to have a significant effect on the prevalence of adverse events.
The study suffers certain limitations that need to be addressed. First, even though the sample size used in this study is among the largest for this type of study, it is still relatively limited and does not provide a representative picture of the population. This means that it may not have sufficient representation of certain genotypes, and it is possible that weaker but significant effects are not being detected. A limited sample size also prevents us from testing interactions between genes or even the development of models to assess prediction models. Furthermore, this study focused its attention on eight SNPs from seven genes, but there are more than forty genes and almost a hundred SNPs that may be involved in the pharmacokinetics of MTX. Nevertheless, our study, as far as we know, is the one that includes the greatest number of covariables and comorbidities. Our results suggest that these covariables are more relevant than genomic information when building predictive models, as shown by the ROC analysis ( Figure 8).

Conclusions
In conclusion, this study has demonstrated that the pharmacogenetic analysis of MTX treatment is more complex than previously thought and that we are still far from having a complete understanding of the process. Our findings have identified an association of ITPA rs1127354 and ABCB1 rs1045642 with response to MTX, polymorphisms of the FPGS rs1544105, GGH rs1800909, and MTHFR genes with disease remission, GGH rs1800909 and MTHFR rs1801131 polymorphisms with all adverse events, and ADA rs244076 and MTHFR rs1801131, rs1801133 polymorphisms with gastrointestinal adverse events. Findings also have highlighted the importance of considering sex-specific effects. Additionally, we have shown that the inclusion of clinical covariates is crucial for obtaining meaningful pharmacogenetic models. Further studies should take these factors into account to better understand the pharmacogenetics of MTX treatment.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/pharmaceutics15061661/s1, Table S1: Response to MTX association. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S2: Response to MTX association by sex. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S3: AR remission association. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S4: AR remission association by sex. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S5: Association with presence of any adverse event. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S6: Association with presence of any adverse event by sex. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S7: Association with presence of gastrointestinal adverse event. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables; Table S8: Association with presence of gastrointestinal adverse event by sex. Log OR, 95% CI and FDR of the logistic regression in shown for all the SNPs analyzed in this study. Analysis were developed using a multivariable logistic regression with 44 covariables. Funding: This research was supported by Instituto de Salud Carlos III, Ministerio de Ciencia e Innovación and Unión Europea-European Regional Development Fund (Proyecto IMP/00009).

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki and approved by the Bioethics and Biosafety Commission of Extremadura (No. 42/2011).

Informed Consent Statement:
Informed consent was obtained from all patients.

Data Availability Statement:
The data presented in this study are available in deidentified form on request from the corresponding author. The data are not publicly available due to privacy restrictions.