Pharmacogenetic model for predicting adverse effects of methotrexate in patients with rheumatoid arthritis

Introduction

According to Russian and European clinical guidelines, methotrexate (MTX) is used as first-line therapy for rheumatoid arthritis (RA) and, if necessary, is combined with other conventional synthetic disease-modifying antirheumatic drugs (csDMARDs) and/or symptomatic anti-inflammatory drugs. Emphasis is placed on the importance of regular laboratory monitoring to avoid adverse reactions (AEs) to MTX. The switch to biologic or targeted synthetic DMARDs is made in cases of MTX intolerance or inefficacy [1-4].

Studies report varying frequencies and patterns of MTX AEs. According to Pincus T et al., at least 72.9% of patients experienced at least one AE [5]. Gastrointestinal reactions (nausea, vomiting, stomatitis, inflammatory and erosive mucosal lesions) were noted in 20-30% of patients, making them the most common limitation in therapy. Hepatotoxicity with elevated liver enzymes (alanine aminotransferase and aspartate aminotransferase, ALT/AST), steatosis, and fibrosis was diagnosed in approximately 10-15% of patients. The development of severe fibrosis and liver cirrhosis was much rarer but required long-term monitoring of liver function. In this context, the recommendations of the Drug-Induced Liver Injury Network (DILIN), emphasizing the need to exclude other causes of liver injury and continuously monitor organ status, are useful [6]. Hematological disorders (leukopenia, thrombocytopenia, anemia) were observed in 5-10% of patients and required regular complete blood count monitoring as well as folic acid supplementation to reduce the risk of severe AEs. Rare but dangerous pulmonary complications, including interstitial pneumonitis, were identified in 1-2% of patients, especially during the first months of MTX administration, and required immediate treatment discontinuation and specialized care [5, 7, 8].

Currently, there is no reliable and clear understanding of all the mechanisms underlying the development of MTX AEs. Particular attention is paid to single nucleotide polymorphisms (SNPs) of proteins involved in MTX metabolism, which affect the pharmacokinetics and pharmacodynamics of the drug. Polymorphisms previously studied by us that affect MTX therapy outcomes are presented in Table 1; the investigation of the association of AEs with SNPs of DHFR (rs70991109, rs70991108) is currently ongoing.

Functionally, the proteins involved in MTX metabolism can be divided into groups: transporters, metabolic enzymes, and regulators of purine and adenosine metabolism. Transporters include ABCB1, a protein that mediates MTX efflux from the cell, and SLC19A1, responsible for folate and MTX influx into the cell. Polymorphisms of these genes affect both drug efficacy and tolerability. Enzymes involved in intracellular MTX metabolism include FPGS, which catalyzes MTX polyglutamation by adding glutamate residues, thereby activating the drug, and GGH, which performs deconjugation by removing glutamate residues. Their ratio determines both the duration of action and the accumulation of active MTX forms, influencing pharmacological effect and treatment safety. A number of enzymes play key roles in the metabolism of adenosine, which possesses anti-inflammatory properties: 1) adenosine deaminase (ADA), which removes an amino group and converts adenosine to inosine, reducing extracellular adenosine levels; 2) adenosine monophosphate deaminase 1 (AMPD1), which transforms adenosine monophosphate (AMP) into inosine monophosphate (IMP) with ammonia release; 3) inosine triphosphatase (ITPA), which protects the cell from purine metabolite toxicity by removing a phosphate group, preventing the accumulation of deoxyadenosine triphosphate; and 4) 5-aminoimidazole-4-carboxamide ribonucleotide transformylase (ATIC), involved in de novo purine synthesis and serving as a key enzyme in adenosine synthesis. MTX exerts an indirect inhibitory effect on enzymes of the adenosine pathway, promoting the accumulation of adenosine inside and outside the cell, thus realizing its anti-inflammatory effect [9]. Key enzymes of folic acid metabolism, essential for DNA synthesis, repair, methylation, and amino acid metabolism, are MTHFR and DHFR. MTHFR converts folate into the active form 5-methyltetrahydrofolate, which participates in the conversion of homocysteine to methionine, while DHFR reduces folic acid to tetrahydrofolate, necessary for MTHFR function. MTX directly inhibits DHFR and indirectly reduces MTHFR activity, thereby disrupting DNA synthesis and cell proliferation, providing its cytotoxic effect [10]. Polymorphisms in the genes encoding the above proteins alter MTX transport, metabolic transformation, and elimination, determining the individual risk of AEs. Consideration of these features allows personalization of therapy, reducing the risk of drug discontinuation due to intolerance. The search for reliable genetic markers is currently ongoing.

Table 1. Key proteins of MTX metabolism and their encoding SNPs

Note: * — according to the old nomenclature.

Materials and Methods

The study was approved by the Local Ethics Committee of Chelyabinsk State Medical Academy of the Ministry of Health of the Russian Federation (Protocol No. 10 dated November 25, 2012). The study included 294 patients with a confirmed diagnosis of RA. The work was conducted as a prospective cohort model with participant enrollment over ten years. Monitoring for MTX AEs in the monotherapy regimen was carried out for 6 months.

MTX AEs were assessed at each visit based on patient complaints, clinical examination (including mucosal status), and laboratory parameters using the Naranjo scale to establish causality. Hepatotoxicity was diagnosed upon persistent elevation of ALT and AST levels (De Ritis ratio 1.33±0.42) with subsequent normalization after drug discontinuation. Leukopenia was defined as a white blood cell count below 3.0×10⁹/L.

Pharmacogenetic analysis focused on SNPs of genes encoding MTX metabolism and transport: ABCB1, ADA, AMPD1, ATIC, FPGS, GGH, ITPA, MTHFR, SLC19A1. These genes cover processes of MTX influx (1 gene) and efflux (1 gene), polyglutamation (1 gene), deconjugation (1 gene), as well as the folate pathway (1 gene) and the adenosine pathway (4 genes), including ATIC, involved in de novo purine synthesis. SNP selection was based on dbSNP data, with confirmed association with MTX efficacy and/or AEs in PubMed or PharmGKB. The minor allele frequency of each SNP in the population was at least 5%. No studies of these SNPs in the context of RA treatment with MTX had previously been conducted in the Russian Federation. The candidate approach allowed the selection of 12 SNPs from 9 genes: one SNP for 6 genes and two for 3 genes. Genomic DNA extraction from peripheral venous blood samples was performed using the commercial kit "Protrans DNA Box 500" (Germany). Genotyping was performed using polymerase chain reaction (PCR). For the analysis of 10 single nucleotide polymorphisms — SLC19A1 rs1051266 (G80A), ABCB1 rs1128503 (C1236T) and rs2032582 (A2677C), GGH rs3758149 (C-401T), FPGS rs84451422 (A>C) and rs1544105 (C105T), ATIC rs2372536 (C347G), ADA rs244076 (T>C), AMPD1 rs17602729 (C34T), ITPA rs1127354 (C94A) — primers domestically produced by TestGen LLC were developed. For MTHFR rs1801131 (A1298C) and rs1801133 (C677T) polymorphisms, commercially available reagents were used. Detection of amplification products was performed by endpoint FLASH analysis on QuantStudio (Applied Biosystems) thermal cyclers. Data processing and interpretation were carried out using QuantStudio Design and Analysis Software (version 1.5.2).

Statistical analysis to identify complex complementary gene relationships, considering possible non-additive interactions and predicting the risk of AEs, was performed using Multifactor Dimensionality Reduction (MDR) to obtain decision rules of the "if–then" type. An exhaustive search algorithm for the best models was employed, along with entropy graphs based on information analysis. The diagnostic efficiency of the obtained models was assessed using sensitivity and specificity indicators, reliability was evaluated based on 10-fold cross-validation results, and statistical significance was assessed using the chi-square criterion. Calculations and graphical constructions were performed using the mdr package (version 3.0.2).

Results

AEs were recorded in 82 patients (27.89% of 294 participants). Some patients experienced a combination of several AEs. The frequency of AEs was higher in non-responders to MTX therapy (n=51; 17.3%) compared to responders (n=31; 10.5%). Three main categories of AEs were identified: hepatotoxicity (n=50; 17.0%), gastrointestinal reactions (n=29; 9.9%), including nausea and vomiting (n=26; 8.8%) and stomatitis (n=3; 1%), as well as leukopenia (n=3; 1%). A combination of hepatotoxicity and stomatitis was observed in two patients.

In the first stage of pharmacogenetic analysis, the frequencies of SNPs of the selected genes were studied in relation to the occurrence of MTX AEs. It was found that the TT genotype of the MTHFR rs1801133 polymorphism was significantly more frequent in patients with AEs: 13.4% (11 cases) vs. 6.1% (13 cases) in the group without AEs (p=0.041, odds ratio (OR) = 2.37, 95% confidence interval (CI) [1.02; 5.54]). For different types of AEs, such as hepatotoxicity, gastrointestinal disorders, and hematological changes, no significant correlations with SNPs were identified. Analysis of allele combinations showed that their frequency ranged from 0.0099 to 0.0569 in the group without AEs and from 0.0089 to 0.0349 in the group with AEs. These data did not allow the identification of a reliable marker for predicting AEs of therapy; therefore, the correlation of SNP combinations with MTX AEs was investigated, and prediction models were generated [11].

Building pharmacogenetic models for predicting AEs requires the selection of markers that provide high accuracy with a minimum number of variables. Model efficiency is determined by its ability to accurately predict outcomes in different samples from the general population. Optimization is achieved by using a limited set of predictors while maintaining high predictive power. To identify genes and their interactions affecting the risk of MTX AEs in patients with RA, the MDR method was applied. The model-building process included sequential stages of selecting genetic factors and assessing their contribution to AE development.

Automated analysis was the first stage in searching for optimal pharmacogenetic models, covering one to three genes, using an exhaustive search algorithm without considering the biochemical roles of the proteins encoded by the studied polymorphisms. As a result, the three best models were identified, including one, two, and three genes respectively: MTHFR rs1801133; FPGS rs1544105 + GGH rs3758149; ABCB1 rs2032582 + MTHFR rs1801133 + FPGS rs4451422 (Table 2). Data analysis revealed that none of these models were suitable for diagnostic application: they were characterized by low support in cross-validation (2/10, 3/10, 6/10, respectively), indicating insufficient reliability, and also had limited diagnostic efficiency not exceeding 64%. Furthermore, the simple model based on the MTHFR rs1801133 polymorphism did not reach statistical significance (p <0.24). It should be noted that three of the five polymorphisms included in the automatic selection models are related to the polyglutamation system (GGH rs3758149, FPGS rs1544105 and rs4451422) and will be included in models created based on information search in the next stage of analysis.

Table 2. Characteristics of models for predicting AEs to therapy based on the results of automated MDR analysis (n~AEs~ = 82, n~comparison group~ = 212)

In the second stage, information analysis was performed, considering the biochemical function of the proteins encoded by the studied genes. To improve the quality of predictive models, a graph illustrating the contribution of genes to the development of AEs was constructed, based on information analysis using Shannon entropy. This measure of uncertainty allowed the determination of information gain as the difference between the probability distributions of the system with and without considering individual genetic elements. The graph is represented by vertices (genes) and edges (their interactions), where the information gain values in percentages reflect the contribution of both individual genes and their pairwise interactions to total entropy. Edge thickness corresponds to the magnitude of the gain, and color visualizes the nature of the interaction: orange and red — synergistic, non-additive (epistatic), gene-effect-enhancing interactions; green and blue — additive interactions with redundant information; brown indicates a weak or independent influence (Fig. 1).

Figure 1. Entropy graph for the contribution of genes and their interactions to the risk of AEs during MTX therapy

The predominance of orange and red shades in the circular graph indicates the dominance of non-additive epistatic gene interactions. Redundancy was weak and limited to interactions with the MTHFR rs1801133 gene (green edges). The contribution of individual polymorphisms to AE risk was minimal (ranging from 0.01% for GGH rs3758149 to 1.02% for MTHFR rs1801133), whereas the influence of gene interactions was significantly higher, reaching 1.73%. The maximal interactions concerned both the transport system (SLC19A1 rs1051266 + ABCB1 rs2032582) at 1.56% and the polyglutamation system (GGH rs3758149 + FPGS rs1544105) at 1.53%, (GGH rs3758149 + FPGS rs4451422) at 1.73%, which is also confirmed by the short branches on the dendrogram (Fig. 2).

Figure 2. Similarity dendrogram for the contribution of genes and their interactions to the risk of AEs during MTX therapy

In predictive models for MTX AEs, the polyglutamation system genes demonstrated significance due to strong interactions, dictating the necessity of their inclusion in all models. It should be emphasized that when modeling therapy inefficacy, these genes were also active; however, the corresponding model was inferior to the final model, which included genes responsible for the adenosine and folate metabolic pathways, as well as intracellular methotrexate transport.

The first model based on information analysis, including three genes of the polyglutamation system (FPGS rs4451422 + FPGS rs1544105 + GGH rs3758149), demonstrated maximal adverse effects in carriers of the CT genotype of the FPGS rs1544105 polymorphism, except for heterozygotes for GGH rs3758149 and FPGS rs4451422 (Fig. 3). An increased risk was also observed in carriers of the combination of ancestral alleles of all three genes (indicated by the dark gray cell in the upper left corner of the figure).

The model with high reliability (10/10), high significance (p=0.003), and good sensitivity (74.7%) had low specificity (44.8%); therefore, the number of polymorphisms was expanded to 5 and 6.

Figure 3. Bar charts of the number of patients in cells of genotype combinations for the three-gene polyglutamation model

Columns on the left — number of patients with AEs, columns on the right — number of patients without AEs; dark gray cells — high-risk combinations, light gray cells — low-risk combinations, white cells — no combination of genotypes present.

The five-gene model, combining the transport and polyglutamation systems (SLC19A1 rs1051266 + ABCB1 rs2032582 + FPGS rs4451422 + FPGS rs1544105 + GGH rs3758149), retained statistical significance (p <0.001) and cross-validation reproducibility (10/10), demonstrating higher sensitivity (79.3%) and specificity (63.2%) compared to the polyglutamation model (Table 3).

Table 3. Characteristics of models for predicting AEs for therapy based on the results of MDR information analysis (n~AEs~ = 82, n~comparison group~ = 212)

The expanded six-gene model, supplemented with the ATIC rs2372536 gene, retained a high degree of validity (p <0.001) and reliability (10/10) with qualitative indicators of diagnostic efficiency (80.4%) and particularly high sensitivity — 91.5% (Table 2). The circular graph reflects the interaction of the ATIC gene with FPGS (rs4451422 — 0.86% and rs1544105 — 1.05%) (Fig. 1), and the dendrogram places the ATIC allele among the three pairs with the greatest mutual synergistic effect (short red branches, see Fig. 2).

Thus, using the MDR method, the contribution of 12 allelic polymorphisms of 9 genes to the risk of developing MTX AEs in the treatment of RA was investigated. The greatest diagnostic value was demonstrated by two models combining genes of the transport system (SLC19A1 rs1051266, ABCB1 rs2032582), polyglutamation (GGH rs3758149, FPGS rs4451422, rs1544105), and ATIC rs2372536 in the six-gene model. The models are distinguished by high statistical significance (p <0.001), cross-validation reliability, and sensitivity/specificity of 79.3% and 91.5% / 63.2% and 69.3%, respectively. The "if–then" rule of these models opens up possibilities for predicting the risk of MTX AEs and implementing them into clinical practice. It is assumed that the results may be extrapolated to methotrexate therapy for undifferentiated arthritis. The six-gene model is patented and ready for practical application [12].

Discussion

In the Russian Federation, the direction associated with the transition to personalized medicine, the development of high-tech healthcare methods, and health-saving measures, including the rational use of drugs, is recognized as a priority [13]. These goals relate to the treatment of RA as the most common and socially significant disease in rheumatology. The study of predictors of therapeutic response to MTX as a first-line treatment for this disease is of primary importance. Predicting its inefficacy and AEs predetermines the choice in favor of other DMARDs. The effectiveness of MTX therapy is partially reflected by clinical markers, but their significance for the risk of developing AEs is low [14]. Pharmacogenetic studies of MTX intolerance predisposition have the greatest predictive value, as demonstrated in our earlier work [15-17].

Recognition of the complexity of the pathogenetic processes of AEs stimulated further research and the generation of a computer program for pharmacogenetic testing in clinical practice. For these purposes, a list of key genetic markers influencing the development of MTX AEs was defined, and a domestic test system was developed and validated by TestGen LLC.

The practical model was based on pharmacogenetic models that consider the combined influence of genes on MTX adverse effects. Models including genes encoding enzymes of transport, conjugation/deconjugation, and the adenosine pathway of MTX action demonstrated the greatest statistical significance. The five-gene model, combining the transport and polyglutamation systems (SLC19A1 rs1051266 + ABCB1 rs2032582 + FPGS rs4451422 + FPGS rs1544105 + GGH rs3758149), reflects the biological mechanisms of methotrexate metabolism: limited drug entry into the cell and impaired conversion to active forms. The influence of adenosine pathway polymorphism on the development of MTX AEs is represented in the six-gene model (SLC19A1 rs1051266 + ABCB1 rs2032582 + FPGS rs4451422 + FPGS rs1544105 + GGH rs3758149 + ATIC rs2372536). Thus, MTX AEs are determined by allelic variants of genes encoding the entire process of MTX metabolism.

Graphical representations of the final models, including 5 and 6 genes, illustrate the nature of their interactions (Fig. 4). Among the transport system genes, a strong interaction between SLC19A1 and ABCB1 is observed. In the polyglutamation system, the GGH polymorphism interacts with both FPGS variants, with no inter-pair links within FPGS. The ATIC alleles (adenosine pathway) exhibit non-additive epistatic interaction with both FPGS variants (glutamation system) and SLC19A1 (transport system). Inter-system links between transport and polyglutamation genes are generally about half as strong compared to intra-system links, suggesting the existence of two relatively independent components of the predictive model. Each component identifies only a portion of at-risk patients, and their joint inclusion in the model allows achieving high diagnostic efficiency.

Five-gene model	Six-gene model

Figure 4. Entropy graph for the final models of gene interactions in the risk of treatment-related AEs.

Graph construction method: Fruchterman–Reingold algorithm.

Practical implementation of this work is possible using a computer program titled "Prediction of Methotrexate Adverse Effects in Rheumatoid Arthritis Based on Patient Genotyping Results." Since the data were obtained from a cohort of the European population, there is a high probability of reproducing the results in most regions of the Russian Federation [12].

Conclusion

The informativeness of a combined analysis of pharmacogenetic markers for predicting tolerability of MTX therapy in patients with RA exceeds the predictive power of individual allelic variants. The diagnostic efficiency of models based on biological data on the role of proteins encoded by the studied genes surpasses that of automated models. Two genetic models for predicting MTX tolerability have been developed, based on polymorphisms encoding MTX transport and conjugation. The model "SLC19A1 rs1051266 + ABCB1 rs2032582 + GGH rs3758149 + FPGS rs1544105 + FPGS rs4451422 + ATIC rs2372536" possesses the greatest diagnostic significance, with a sensitivity of 91.5% and specificity of 69.3%. In clinical practice, it is proposed to apply the model in the form of a formal "if–then" rule to predict the development of MTX AEs.