Genetic variations in methotrexate metabolic pathway genes influence methotrexate responses in rheumatoid arthritis patients in Malaysia

Methotrexate (MTX) is the most widely used disease-modifying anti-rheumatic drug (DMARD) for rheumatoid arthritis (RA). Many studies have attempted to understand the genetic risk factors that affect the therapeutic outcomes in RA patients treated with MTX. Unlike other studies that focus on the populations of Caucasians, Indian and east Asian countries, this study investigated the impacts of six single nucleotide polymorphisms (SNPs) that are hypothesized to affect the outcomes of MTX treatment in Malaysian RA patients. A total of 647 RA patients from three ethnicities (NMalay = 153; NChinese = 326; NIndian = 168) who received MTX monotherapy (minimum 15 mg per week) were sampled from three hospitals in Malaysia. SNPs were genotyped in patients using TaqMan real-time PCR assay. Data obtained were statistically analysed for the association between SNPs and MTX efficacy and toxicity. Analysis of all 647 RA patients indicated that none of the SNPs has influence on either MTX efficacy or MTX toxicity according to the Chi-square test and binary logistic regression. However, stratification by self-identified ancestries revealed that two out of six SNPs, ATIC C347G (rs2372536) (OR 0.5478, 95% CI 0.3396–0.8835, p = 0.01321) and ATIC T675C (rs4673993) (OR 0.5247, 95% CI 0.3248–0.8478, p = 0.008111), were significantly associated with MTX adequate response in RA patients with Malay ancestry (p < 0.05). As for the MTX toxicity, no significant association was identified for any SNPs selected in this study. Taken all together, ATIC C347G and ATIC T675C can be further evaluated on their impact in MTX efficacy using larger ancestry-specific cohort, and also incorporating high-order gene–gene and gene–environment interactions.

www.nature.com/scientificreports/ have been in use since the 1950s 5 . Since methotrexate (MTX) was readapted in the late 1980s, it has become the most widely used csDMARD 6 . MTX is a folate anti-metabolite that suppresses disease activity and reduces joint pain. A low dose (15-25 mg) of MTX per week is prescribed to patients through either subcutaneous or oral administration for at least three months and this drug has been proven to be an effective DMARD for RA 7,8 . However, about 30-50% of RA patients do not respond to MTX and thus ruling out MTX as a treatment option for these non-responders [9][10][11][12][13] . Moreover, up to 35% RA patients are forced to discontinue MTX due to the adverse drug effects including stomatitis, gastrointestinal upset, headache or minor central nervous system disturbance, hair loss, ulcers, liver toxicity, pancytopenia and pneumonitis 6,[9][10][11][12][13] . Despite these limitations, MTX is still the gold standard for the treatment of RA and in Malaysia, the usage of MTX has increased sixfold from 1997 to 2007 14 . Malaysia Clinical Practice Guidelines on RA recommends that patients are to take a combination of MTX with non-biologic DMARDs or folic acid (minimum 5 mg/week) when MTX monotherapy shows signs of failure in efficacy or side effects, respectively 15 . Hence, anti-RA medicine can become unexpectedly lengthy, costly and ranges from not effective to partially effective in meeting treatment expectation.
Available literature has converged to hypothesize that the variability in MTX efficacy and toxicity is due to a dysregulation in the MTX pathway ( Supplementary Fig. S1) [16][17][18] . As multiple enzymes mediate the metabolism of MTX, it is conceivable that the alterations of enzymes' availability and activity have a direct impact on MTX treatment. Our study focused on the potential dysregulation of 4 key enzymes: folylpolyglutamate synthase (FPGS), γ-glutamyl hydrolase (GGH), aminoimidazole-4-carboxamide ribonucleotide (AICAR) transformylase, and inosine triphosphate pyrophosphatase (ITPA). FPGS is an important enzyme responsible for converting MTX into a range of polyglutamate forms-methotrexate polyglutamate (MTX PG) and genetic variants of FPGS likely affect the retention of MTX bioavailable in the cell. The therapeutic effects of MTX in RA patients rely on its conversion to MTX PG 18 . The conversion of MTX to MTX PG by FPGS can be reversed by GGH 16 ( Supplementary Fig. S1). GGH removes the polyglutamates from MTX PG by performing serial trimming on the long chain MTX PG to yield MTX which is able to be exported from the cell. Therefore, variants of GGH such as rs11545078 (C452T) causes defective GGH and leads to an increase of intracellular concentration of MTX PG 19 , rs3758149 (C-401 T), on the other hand, promotes enzymatic activity of GGH that results in the hydrolysis of MTX PG 20 . AICAR transformylase (encoded by ATIC gene) converts AICAR to formyl-AICAR (FAICAR), playing a role in the de novo purine synthesis ( Supplementary Fig. S1). MTX is able to block AICAR transformylase and hence causes an increase of intracellular level of AICAR which has a role in the activation of adenosine signaling pathway, resulting in a series of anti-inflammatory activities 21 . Several studies shown two variants of ATIC, C347G (rs2372536) 20,22 and T675C (rs4673993) 23,24 achieve good clinical response to MTX. ITPA is also played an important role in the adenosine signaling pathway. The presence of C94A (rs1127354; chr20:3213196) mutation may modulate the ability of ITPA to convert inosine triphosphate (ITP) to inosine monophosphate (IMP) required for de novo purine synthesis which will affect the amount of intracellular adenosines which can be exported for binding to adenosine receptors 21 , leading to multiple anti-inflammatory mechanisms.
Since FPGS, GGH, AICAR and ITPA are known to determine the rate limiting steps in MTX metabolism, the effects of genetic variations such as single nucleotide polymorphisms (SNPs) could modulate the pharmacokinetic properties of these enzymes 17,23,[25][26][27] . In our study, six SNPs were selected for the genotyping analysis as based on their previously reported association with MTX treatment outcomes. These SNPs include FPGS A1994G (rs10106; chr9:127813796), GGH C452T (rs11545078; chr8:63026205), GGH C401T (rs3758149; chr8:63039169), ATIC C347G (rs2372536; chr2:215325297), ATIC T675C (rs4673993; chr2:215347616) and ITPA C94A (rs1127354; chr20:3213196). In addition, these candidate gene association studies were previously and mainly conducted in the Caucasian population and no similar study has been carried out for Malaysian RA patients. Therefore, the novelty of this study is to reassess the possible association of the SNPs with MTX treatment outcome in a Malaysian population. We hypothesize that specific SNP changes that can alter gene function are able to explain the variability observed in MTX efficacy and toxicity. Thus, this study aimed to genotype six SNPs of the candidate genes (FPGS, GGH, ATIC and ITPA) in MTX metabolic pathway and determine their association with MTX therapeutic outcomes in Malaysian RA patients.

Results
Characterization of the studied population. Demographics. This study recruited 647 RA patients from Sunway Medical Centre (n = 268), Hospital Tuanku Ja'afar Seremban (n = 284) and Hospital Selayang (n = 95) ( Table 1). The number of Chinese RA patients, most of which were recruited from Sunway Medical Centre, was twice the amount compared to the number of Malay or Indian RA patients. This account for the skewed distribution of the ancestry-specific groups in our patient sampling with over 80% of patients recruited from the private hospital (Sunway Medical Centre) were Chinese (N = 326), compared to the number for Malay (N = 153) or Indian (N = 168) patients. Shahrir et al. 28 was the first to report on the number of RA registries in Malaysia, showed majority of RA patients in Malaysia are Indian (54.5%), followed by Malay (31.4%), Chinese (11.6%), Indigenous (1.2%) and others (1.3%). However, the authors also admitted the possibility of incomplete data because they have not included RA patients from private and university hospitals in their study. In 2008, the population sizes of Malay, Chinese and Indian were 53.3%, 26.0% and 7.7%, respectively. Due to the lack of information of the prevalence of RA in different ancestry-specific groups in Malaysia, and due to limitation in the sampling of patients, we cannot conclusively determine the prevalence of RA in the three main ancestryspecific groups of Malays, Chinese and Indians in Malaysia. The female RA patients (88.7%) outnumbered male RA patients (11.3%) in this study. This sex-imbalance is consistent with the current literature which have shown a higher number of female RA patients.  Table 2). In addition, genotype counts for the six SNPs were compared among the three ancestry-specific groups using chi-square test. The results revealed that except GGH C452T (rs11545078) and ITPA C94A (rs1127354), other four SNPs significantly differ in genotype frequencies among the Malay, Chinese and Indian RA patients (dashed line boxes in Table 2).
Association of six metabolic SNPs with MTX efficacy and toxicity in three ancestry-specific RA patients. When the association study of SNPs with the MTX treatment was carried out using the entire cohort (n = 647), there was no significant difference between ARs and IRs as well as Non-ADR groups and ADR groups in the allelic association tests. Logistic regression was then performed to test the standard models of disease penetrance (dominant, recessive, additive) for the interaction of six SNPs with MTX efficacy and toxicity www.nature.com/scientificreports/ in the cohort of 647 RA patients. The forest plot ( Fig. 1) for the association between SNPs and the MTX efficacy and toxicity was performed using R package ggplot2 29,30 indicated no significant association between the SNPs with either MTX efficacy or MTX toxicity (Table 3). We then stratified the study cohort into three ancestry-specific groups, Malay, Chinese and Indian. The numbers of AR and IR of the stratified groups with their respective genotypes are presented in Table 4. Using these data, significant differences were observed in ATIC C347G (rs2372536) (OR 0.5478, 95% CI 0.3396-0.8835, p = 0.01321) and ATIC T675C (rs4673993) (OR 0.5247, 95% CI 0.3248-0.8478, p = 0.008111) ( Table 5). Based on the effect sizes obtained from our analyses, the risk of patients becoming IR was reduced by ATIC C347G (rs2372536) and ATIC T675C (rs4673993) for approximately 55% and 57%, respectively. When the inheritance models were applied to the ancestry-specific stratification, it could be inferred that: (i) ATIC C347G (rs2372536) was associated with AR in Malay RA patients under dominant and additive models; (ii) the minor allele of ATIC T675C (rs4673993) under three genetic models (dominant, recessive and additive) may predict a higher success rate in MTX treatment among Malay RA patients ( Table 5). All positive results for the association between SNPs and MTX efficacy are as shown in the forest plot (Fig. 2). All the six SNPs were not significantly associated with the MTX efficacy in either Chinese or Indian RA patients. Furthermore, there was no significant association of all six SNPs with MTX toxicity in the three ancestry-specific groups (Tables 4, 6). It is worth mentioning a recent paper concluded an association of MTX side effects with FPGS A1994G (rs10106) of the south Indian Tamil with RA 31 . More than 80% of Malaysian Indian study subjects are also of Tamil origin, despite this the association between MTX toxicity and the FPGS A1994G (rs10106) in this ancestry-specific group is not evident in our study.

Discussion
Majority of the available drugs used for the treatment of RA were clinically evaluated in European ancestries, this raises a concern about their efficacy and toxicity in other ancestry groups globally 32 . Asian populations especially those in South East Asia were considerably under-represented in pharmacogenomic and pharmacogenetic studies of RA 32,33 . Hence, the present study evaluated the outcomes of MTX treatment in three major ancestry groups in Malaysia and their association with 6 SNPs from the enzymes involved in MTX metabolism. Comparing with the studies conducted by geographical locations, our study attempted to delineate ancestry specific risk factors that would increase the precision of the proposed association.
The MAF of ATIC T675C (rs4673993) recorded in Genome Aggregation Database (gnomAD) and 1000 Genomes is 0.3251 and 0.2855, respectively 34,35 . In our study, the allele frequency of ATIC T675C (rs4673993) for the overall cohort is 0.4 and it is 0.44, 0.31 and 0.54, respectively, in Malay, Chinese and Indian populations in Malaysia (Table 2). By comparing the allele frequency between our study and the public database, we noticed Table 2. Allele and genotype frequencies of the three ancestry-specific groups in Malaysia. *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001; ****P ≤ 0.0001. a Data are presented in number (percentage) unless otherwise indicated. b Data are presented in frequency (number) unless otherwise indicated. www.nature.com/scientificreports/ that our population is carrying a higher allele frequency of ATIC T675C (rs4673993). When comparing allele frequencies by ethnicity within our study cohort, there was a significant difference between Malay, Chinses and Indian for this SNP. Apparently, the allele frequency of ATIC T675C (rs4673993) in Indian and Malay subjects were significantly higher than that observed in Chinese subject and in the public databases. Interestingly, our study suggested that the Malay RA patients with minor allele of ATIC T675C (rs4673993) have a better treatment outcome upon MTX monotherapy. In other words, this minor allele was associated with an increased remission rate in Malay RA patients following the treatment of MTX. A few studies have also demonstrated the impact of ATIC T675C (rs4673993) on MTX treatment outcome. Prospective studies conducted by Lee et al. and Iannaccone et al. have shown that the minor allele of ATIC T675C (rs4673993) was significantly associated with low disease activity in RA patients with MTX monotherapy 23,24 . These two studies were conducted in the USA with 120 and 262 RA patients, respectively, as the subsets of Brigham and Women's Hospital Rheumatoid Arthritis Sequential Study (BRASS). Moreover, a meta-analysis performed by Chen (2017) indicated that ATIC T675C (rs4673993) could predict the responsiveness of MTX treatment in recessive model (OR 2.54; 95% CI 1.23-5.26) 36 . The authors combined two studies to yield a total of 698 Caucasians and observed a significant favouritism of ATIC T675C (rs4673993) in RA patients having response to MTX treatment. On the other hand, a retrospective study by Lima et al. (2014) gave a totally different conclusion 37 , whereby more than fourfold increase in risk of MTX inefficacy was associated with ATIC T675C (rs4673993) in a population of 233 adults ( ≥ 18 y.o.) of Portuguese Caucasian RA patients. The result discrepancy may be due to different ancestral lineages of RA patients enrolled in the respective studies. In our case, the association of polymorphism ATIC T675C (rs4673993) with the responsiveness to MTX treatment could only be observed in Malay but not in Chinese and Indian RA patients.
Current RA literature consistently highlights the hypothesis that the anti-inflammatory action of MTX is achieved through the indirect inhibition of AICAR. The ATIC T675C (rs4673993) SNP is positioned in the intronic region of ATIC. To our knowledge, there are no functional studies on this particular SNP in ATIC activity. Nevertheless, the intronic SNP either interferes the transcriptional regulation of the coding-enzyme or is in linkage disequilibrium (LD) with another coding SNP 23,[38][39][40] . In the present study, since similar effect size of ATIC T675C (rs4673993) and ATIC C347G (rs2372536) was observed, both SNPs can be in LD. Nevertheless, this observation need further validation, since the current sample size is too small to perform a LD test and the lack of a reference panel for ATIC T675C (rs4673993) and ATIC C347G (rs2372536) in Malay patients.    22 studied 422 Caucasian RA patients in Poland who were treated with MTX therapy and found that GG minor genotype significantly exhibited a good response to MTX. However, the lack of association between rs2372536 polymorphism and the clinical response to MTX was also reported in some studies 42,43 . Recently, two meta-analyses were performed to investigate the association between ATIC C347G (rs2372536) and MTX response 44,45 . The first meta-analysis was based on five studies of 1056 RA patients in which 722 were MTX responders and 334 were non-responders. This analysis found the difference of ATIC C347G (rs2372536) between Caucasions (Spain, Slovenia and Netherlands) and Asians (India), being associated with www.nature.com/scientificreports/ non-responsiveness to MTX treatment in Caucasians but not associated in Asians 44 . The second meta-analysis combined two European (Spain and Netherlands), one East Asian (Japan) and two South Asian (India) studies with 458 MTX responders and 398 non-responders in total 45 . When combining five studies, ATIC C347G (rs2372536) demonstrated a significant association with non-responsiveness of MTX under the dominant and codominant models. Yet, geographical stratification showed that the association of ATIC C347G with MTX response was still observed in Europeans in pre-allele, dominant and codominant models but not in South Asian populations 44 . Despite all studies above demonstrated a significant association between ATIC variants and MTX efficacy, the results were raher inconsistent. Common factors for inconsistency such as small sample size and insufficient statistical power, study design, medication dosage, grouping criteria, and patient condition could cause limitations in the association study. Moreover, gene-gene interactions within folate and adenosine biosynthesis pathways may complicate the association study between SNPs and MTX treatment outcomes 46 . In fact, RA has complex inheritance patterns and no single genetic variant has a decisive role in MTX efficacy or MTX toxicity in the treatment of RA. By using the MDR (Multifactor Dimensionality Reduction) method, a cohort of 255 RA patients treated with MTX in the USA was evaluated with the efficacy of MTX treatment, and the results showed that 53% MTX responders was associated with high-order interactions among SNPs in ITPA (C94A), RFC1 (G80A), and ATIC (C347G) genes 46 . Upon excluding the predisposing genotype combinations, a 3.8-fold lower efficacy was observed 46 . Later, the same researchers extended their study of gene-gene interactions using ITPA (C94A), RFC1 (G80A), and ATIC (C347G) to another 3 RA cohorts (USA, Dutch and Swedish) 47 . Both USA and Dutch cohorts (n = 435) confirmed a predisposing genetic attribute significantly associated with methotrexate response [odds ratio (OR) = 2.9, 95% confidence interval (CI) 1.9-4.2; P < 0.001]. Although the association of combined SNPs with MTX responsiveness in the Swedish cohort (n = 530) could not be determined, the association was observed after the non-genetic factors, age, sex and anti-citrullinated protein antibody (ACPA) status were included in MDR analysis 47 . Thus, individual variants of ATIC may not play a direct role in MTX efficacy, future studies shall map the ATIC variants to drug response as based on the detection of nonlinear multigene interactions, this may improve the accuracy of predicting the MTX efficacy. In addition, other non-genetic covariates should be considered because the association study between genetic variants and MTX efficacy sometimes seems oversimplified understanding the MTX response in RA.
AICAR transformylase contains two domains which are MGS (methylglyoxal synthetase) like domain and AICAR binding domain 48 . ATIC C347G (rs2372536) causes the substitution of threonine (Thr) with serine (Ser) at position 116. Thr116 lies in the binding pocket of MGS-like domain and is the first residue of α8 helix which Table 5. The asscoacition between SNPs and MTX efficacy in three different ancestry-specific groups. a OR: odd ratio; b 95% CI: 95% confidence intervals frequency; c p-value: probability value; d No association analysis is possible. www.nature.com/scientificreports/ likely serves as a N-cap residue stabilizing the helix by interacting with the amide groups from the main chain. We proposed that the side-chain hydroxyl group of Thr116 forms hydrogen bonds with the amide groups of Val117 and Glu118 (green arrowhead in Supplementary Fig. S2), while its main-chain carboxyl group forms hydrogen bond with the amide group from Glu119 (blue arrowhead in Supplementary Fig. S2). The methyl group of Thr116 might stabilize the hydrogen bond between Thr116 and the main chain. As Thr116 is substituted with serine, the methyl group can be removed and this results in a more flexible C-N rotation. In other words, Ser116 causes the rearrangement of the protein structure at N-cap and thus, potentially affects the substrate-binding affinity and AICAR transformylase enzyme activity. This explains why the RA patients with the minor allele of ATIC C347G (rs2372536) might have a phenotypic change in response to MTX.

Conclusion
Present study suggested that minor allele of ATIC C347G (rs2372536) and ATIC T675C (rs4673993) could influence the response to MTX monotherapy in Malay patients with RA, while the other four SNPs failed to demonstrate their associations with the reduction of disease activity following the MTX monotherapy. ATIC C347G (rs2372536) and ATIC T675C (rs4673993) are not the only ancestry-specific SNPs since any variations appear in the genes of MTX metabolic pathway are potentially able to affect the effectiveness of MTX treatment. As more Malay-specific SNPs can be revealed, the prediction of poor response would enable patient to be placed on alternative drugs, while those with predicted good response could proceed with MTX treatment. As for the future recommendation, ancestry specific signal of ATIC should be validated in a larger replication cohort of a similar ancestry group profile to reduce the Type II error rate of MTX treatment response. Both ATIC C347G (rs2372536) and ATIC T675C (rs4673993) warrant an in-depth investigation, especially in the , which is also known as SLC19A1, transports MTX into the cell while ABCB1 (ATP Binding Cassette B1 or P-glycoprotein) is ATP-dependent pump that moves MTX out of the cell. Gene variants of SLC19A1/RFC-1 (rs1051266 and rs2838956) and ABCB1 (rs1045642 and rs1128503) have been shown to influence MTX influx and efflux in and out of the cell, respectively, leading to the therapeutic effects [49][50][51][52] . MTHFR (methylenetetrahydrofolate reductase) is another well studied metabolic gene of MTX and it has been associated with the MTX inadequate responders. Nucleotide substitutions in MTHFR such as C677T (rs1801133) and A1298C (rs1801131) which result in single amino acid substitutions can greatly reduce the production of functional reductase and its enzymatic activity 53 . Last but not least, any variations in enzymes involved in the intracellular MTX metabolic pathway may affect functional properties of the enzyme and hence the responses to MTX efficacy and toxicity.  www.nature.com/scientificreports/ 28 (DAS28CRP) for at least 6 months. On the other hand, inadequate responders have the same RA treatment as adequate responders but failed to achieve clinical remission or low disease activity as defined by DAS28CRP but at present been treated with other DMARDs, mono, duo (excluding MTX and HCQ group) or triple therapy, either csDMARDs or tsDMARDs or bDMARDs. Hence, patients may or may not be on MTX at the point of recruitment.

Methods
MTX toxicity. Patients were categorized into two groups for potential toxicity and side effect: non-adverse drug reaction (Non-ADR) group and adverse drug reaction (ADR) group. The categorization was marked based on whether they have experienced drug intolerance during the MTX treatment. All side effects were recorded from the start of the MTX treatment until the withdrawal due to adverse drug reactions.
Blood sample and clinical data collections. A total of 647 RA patients were involved in this study after their informed consent was obtained. A total of 5 ml of the whole blood sample was collected in ethylenediaminetetraacetic acid (EDTA) tube from individual patients by venepuncture during patients' regular visits at the hospitals. Besides, clinicopathological and demographic data were also extracted from patients' clinical records and linked to deidentified patient blood samples collected for this study.
TaqMan® SNP genotyping assays. Genomic DNA from patients was isolated from patients' whole blood samples. Briefly, buffy coat was obtained from the blood sample by centrifugation. Genomic DNA was then extracted from the buffy coat by using QIAamp DNA Blood Mini Kit (Qiagen, Germany) according to the manufacturer's instructions. SNP genotyping of FPGS A1994G (rs10106), GGH C452T (rs11545078), GGH C401T (rs3758149), ATIC C347G (rs2372536), ATIC T675C (rs4673993) and ITPA C94A (rs1127354) were performed by using TaqMan® SNP Genotyping Assays (Thermo Fisher, USA) according to manufacturer's instructions. The genotype data of each participant were analyzed using an online software named "Genotyping V4.2" (Thermo Fisher Connect™). A total of 5% of the samples (n = 33) for each respective SNP were randomly selected for PCR amplification (AmpliTaq Gold™ 360 Master Mix, Thermo Fisher Scientific, USA) (Supplementary Table S1) and subsequently for the genotype verification by Sanger sequencing (1st BASE Pt Ltd). Sequencing results were curated with SnapGene V4.3.10 (from GSL Biotech; available at https:// snapg ene. com).

Statistical analysis.
Genotype and allele frequency of all the selected six SNPs were calculated. A chisquare independence test was performed to test the association between SNPs and ethnic groups. Chi-square test and binary logistic regression were performed to investigate the association between SNPs and MTX efficacy and MTX toxicity (PLINK V1.09) 54 . Effect sizes of potential associations were calculated as odds ratio (OR) and 95% confidence intervals (CI) as a measure of the association between the categorical variables. A p-value of < 0.05 was considered to be statistically significant.

Data availability
The data underlying this article cannot be shared publicly due to the privacy of individuals that participated in the study. The data will be shared on reasonable request to the corresponding author.