Predicted expression of genes involved in the thiopurine metabolic pathway and azathioprine discontinuation due to myelotoxicity.

TPMT and NUDT15 variants explain less than 25% of azathioprine-associated myelotoxicity. There are 25 additional genes in the thiopurine pathway that could also contribute to azathioprine myelotoxicity. We hypothesized that among TPMT and NUDT15 normal metabolizers, a score combining the genetically predicted expression of other proteins in the thiopurine pathway would be associated with a higher risk for azathioprine discontinuation due to myelotoxicity. We conducted a retrospective cohort study of new users of azathioprine who were normal TPMT and NUDT15 metabolizers. In 1201 White patients receiving azathioprine for an inflammatory disease, we used relaxed Least Absolute Shrinkage and Selection Operator (LASSO) regression to select genes that built a score for discontinuing azathioprine due to myelotoxicity. The score incorporated the predicted expression of AOX1 and NME1. Patients in the highest score tertile had a higher risk of discontinuing azathioprine compared to those in the lowest tertile (hazard ratio [HR] = 2.15, 95% confidence interval [CI] = 1.11-4.19, p = 0.024). Results remained significant after adjusting for a propensity score, including sex, tertile of calendar year at initial dose, initial dose, age at baseline, indication, prior TPMT testing, and the first 10 principal components of the genetic data (HR = 2.11, 95% CI = 1.08-4.13, p = 0.030). We validated the results in a cohort (N = 517 non-White patients and those receiving azathioprine to prevent transplant rejection) that included all other patients receiving azathioprine (HR = 2.00, (95% CI = 1.09-3.65, p = 0.024). In conclusion, among patients who were TPMT and NUDT15 normal metabolizers, a score combining the predicted expression of AOX1 and NME1 was associated with an increased risk for discontinuing azathioprine due to myelotoxicity.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Azathioprine is an immunosuppressant that causes myelotoxicity in some people. Clinical Pharmacogenetics Implementation Consortium (CPIC) guidelines provide azathioprine dosing recommendations based on TPMT and NUDT15 genotype; however, these genotypes explain only 25% of azathioprine‐induced myelotoxicity.

WHAT QUESTION DID THIS STUDY ADDRESS?

The aim of this study was to determine if a risk score composed of the genetically predicted expression of genes that encode proteins in the thiopurine pathway within the liver tissue would be associated with azathioprine discontinuation attributed to myelotoxicity.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study showed that a risk score composed of genetically predicted risk expression of AOX1 and NME1 is associated with azathioprine discontinuation due to myelotoxicity.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

Having a risk score for discontinuation composed of the genetically predicted expression of NME1 and AOX1 could help discriminate patients at high risk of discontinuing azathioprine due to hematologic side effects in people who are normal TPMT and NUDT15 metabolizers.

INTRODUCTION

Azathioprine is a thiopurine used to treat several inflammatory conditions, but it is often discontinued due to dose‐dependent myelotoxicity. , Enzymes in the thiopurine metabolic pathway may have a combined effect on the azathioprine discontinuation risk. Currently, the Pharmacogenomics Knowledgebase (PharmGKB) lists 27 genes involved in thiopurine metabolism and transport, but only thiopurine S‐methyltransferase (TPMT) and nudix hydrolase 15 (NUDT15) are included in clinical guidelines for thiopurine use. Whereas certain TPMT and NUDT15 polymorphisms are associated with several‐fold increased risk for myelotoxicity among azathioprine users, , they explain fewer than 25% of bone marrow toxicity cases in routine clinical practice.

The Clinical Pharmacogenetics Implementation Consortium (CPIC) provides evidence‐based recommendations regarding pharmacogenetic tests for patient care. Their guidelines for azathioprine include dose reductions or alternative medication options based on TPMT and NUDT15 metabolizer status. Other enzymes involved in the thiopurine metabolic pathway are not included in the guidelines because studies have been small, inconclusive, or contradictory. , , , In addition, there are no data on the combined role of enzymes, nor the genes encoding them, on azathioprine toxicity.

PrediXcan uses large‐scale transcriptome datasets that are linked to genetics (e.g., Genotype‐Tissue Expression [GTEx] project ) to generate models for calculating the contribution of genetic variants to gene expression. The prediction models can then be applied in any dataset with genomewide interrogation of common variants. Genetically regulated expression accounts for part of the interindividual variability in measured transcript levels. Thus, PrediXcan yields a gene‐based test that is mechanistic by design. PrediXcan has been used previously in a pharmacogenomic study to investigate efavirenz‐related adverse events.

We hypothesized that a score combining the genetically predicted gene expression of proteins in the thiopurine metabolic pathway would be associated with a higher risk of azathioprine discontinuation due to myelotoxicity in TPMT and NUDT15 normal metabolizers.

METHODS

Data collection

This study was reviewed by the Vanderbilt University Medical Center’s (VUMC) Institutional Review Board and determined to be non‐human subjects research (IRB #180498). It was conducted in BioVU, a clinical practice‐based biobank at VUMC. In brief, BioVU includes de‐identified electronic health records (EHRs) with access to demographic characteristics, clinical notes, medical history, problem lists, medications, and diagnostic and procedure codes ; it is linked with stored DNA samples. , Within BioVU, we assembled a retrospective cohort of patients receiving azathioprine who had been genotyped using the Expanded Illumina Multi‐Ethnic Genotyping Array (MEGA)EX platform. We reviewed their medical records to confirm azathioprine use and included only new users of azathioprine who passed genotype quality control (pre‐ and post‐imputation). We defined new users as individuals that had no prior mention of azathioprine or mercaptopurine use in their EHRs. We used EHRs to collect clinical variables, including reported race, sex, age at azathioprine initiation, initial daily dose of azathioprine, indication, calendar year of initial dose, baseline white blood cell count (closest measure to initial dose within 365 days prior to and including initial dose date), reason for discontinuing, and date of last known dose before end of follow‐up. We used genotype data and classified patients based on TPMT and NUDT15 metabolizer status as per CPIC guidelines and then excluded poor, intermediate, and indeterminate metabolizers (Table S1).

Discovery cohort

The discovery cohort (N = 1201) included patients receiving azathioprine prescriptions for an inflammatory condition, such as Crohn’s disease, ulcerative colitis, vasculitis, systemic lupus erythematosus, or rheumatoid arthritis. We further limited this cohort to individuals reported as White because the GTEx project primarily includes individuals of European ancestry.

Validation cohort

The validation set included patients whose race was not reported as White (i.e., Black, Asian, other, or unknown) or who were taking azathioprine for noninflammatory indications (i.e., organ transplant, other, or unknown; N = 517).

Follow‐up and outcome

Patients entered the cohort on their first mention of azathioprine use in the EHRs, and follow‐up ended on the first of the following dates: (1) day of discontinuation; (2) last confirmed azathioprine prescription or use as per EHR +90 days; (3) lost to follow‐up; (4) day of death; or (5) end of the study (December 31, 2018). The primary study outcome was her‐confirmed azathioprine discontinuation attributed to myelotoxicity, recorded as leukopenia, neutropenia, thrombocytopenia, pancytopenia, and/or anemia. Blinded to genotype data, we reviewed clinician notes and laboratory results to make this determination.

Genetically predicted gene expression

We collected genotype information using the MEGAEX platform, which includes more than two million markers. We prepared genotyping data for imputation using McCarthy tools and imputed additional variants using Michigan Imputation Server with HRC version r1.1 reference panel and phasing with Eagle. Following standard quality control, as described previously, we estimated genetically predicted gene expression of the candidate proteins in the thiopurine metabolic pathway using PrediXcan (GTEx version 8) with MASHR version 8 expression quantitative trait locus (eQTL) weights to perform the imputation. , , , Because most drug metabolism occurs in the liver, we prespecified use of liver tissue‐specific estimations. Of the 27 genes in PharmGKB, we were able to estimate the genetically regulated expression of 19 of them. Four of the genes are involved in transport, 13 genes are involved in the metabolism of azathioprine, and two genes that can cause the cytotoxic effects of azathioprine. The transport genes are ATP‐binding cassette subfamily C member 4 (ABCC4), solute carrier family 28 member 2 (SLC28A2), solute carrier family 29 member 1 (SLC29A1), and solute carrier family 29 member 2 (SCL29A2); the metabolism genes are adenosine kinase (ADK), aldehyde oxidase 1 (AOX1), guanine deaminase (GDA), guanine monophosphate synthase (GMPS), glutathione transferases A1 (GSTA1), glutathione transferases A2 (GSTA2), glutathione transferases (GSTM1), inosine monophosphate dehydrogenase 1 (IMPDH1), inosine triphosphatase (ITPA), NME/NM23 nucleoside diphosphate kinase 1 (NME1), NME/NM23 nucleoside diphosphate kinase 2 (NME2), 5’‐nucleotidase, cytosolic II (NT5C2), and ribonucleotide reductase regulatory subunit M2 (RRM2); and the other two genes are Rac family small GTPase 1 (RAC1) and phosphoribosyl pyrophosphate amidotransferase (PPAT). We were unable to calculate estimates for ATP binding cassette subfamily C member 5 (ABCC5), hypoxanthine phosphoribosyltransferase 1 (HPRT1), phosphoribosyl pyrophosphate synthetase 1 (PRPS1), ribonucleotide reductase catalytic subunit M1 (RRM1), solute carrier family 28 member 3 (SLC28A3), and xanthine dehydrogenase (XDH) because there was insufficient expression data in the reference panel.

Statistical analysis

Demographic and clinical characteristics are presented as numbers and percentages for categorical variables and are presented as mean and standard deviation for continuous variables. We used Fisher’s exact tests to compare binary variables, Pearson’s chi‐squared tests to compare polytomous variables, and Wilcoxon’s rank sum tests to compare continuous variables.

The limited number of outcomes did not provide sufficient statistical power to use all genes in the azathioprine pathway to build a risk score and analyze the association; therefore, we used a variable selection method called relaxed Least Absolute Shrinkage and Selection Operator (LASSO) regression Relaxed LASSO is a two‐step process that employs penalized regression to select variables with a non‐zero coefficient and then uses a multivariate logistic regression model to build a risk score with those variables, avoiding over‐penalization. We developed this risk score for azathioprine discontinuation attributed to possible myelotoxicity in the discovery cohort. Once we had the risk score, we used a Cox hazard regression model to compare the risk of azathioprine discontinuation for possible myelotoxicity in patients by score tertile. Using the already estimated regression coefficients, we calculated the risk score in the validation cohort and defined the risk of azathioprine discontinuation with the Cox hazard regression model. We also adjusted by using a propensity score, which included sex, calendar year tertile of initial dose date, initial dose, age at baseline, indication, prior TPMT testing, and the first 10 principal components of the genetic data. Due to power limitations, we were unable to use the last azathioprine dose in the propensity score; therefore, in a sensitivity analysis, we exchanged the initial azathioprine dose for last azathioprine dose and reran the adjusted analysis.

We completed genotyping quality control steps using PLINK versions 1.9 and 2.0, and R version 3.6.2. , The 10 principal components were computed using PLINK (version 1.9) among individuals who passed genetic quality control after pruning single‐nucleotide polymorphisms in approximate linkage disequilibrium, removing variants with a minor allele frequency less than or equal to 0.01, and removing variants known to affect principal component calculations (HLA, Inversion 8, and Inversion 17). All analyses were conducted using STATA version 17.0.

RESULTS

The discovery cohort included 1201 White patients who were new users of azathioprine, TPMT and NUDT15 normal metabolizers, and taking the medication for inflammatory conditions; they were followed over a mean of 3.22 ± 3.85 years. Their mean age was 44.5 ± 17.6 years, and 66% were women. Forty‐seven users discontinued azathioprine due to attributed myelotoxicity (Table 1).

TABLE 1

Baseline characteristics of patients

	Cohort
	Discovery N = 1201		Validation N = 517
Female sex	795	66.2%	299	57.8%
Azathioprine indication
Systemic lupus erythematosus	115	9.6%	38	7.4%
Other connective tissue diseases	658	54.8%	92	17.8%
Inflammatory bowel disease	428	35.6%	42	8.1%
Organ transplant	0	0.0%	275	53.2%
Other	0	0.0%	69	13.3%
Unknown	0	0.0%	1	0.2%
Race: White	1201	100.0%	302	58.4%
Years of follow‐up: mean (SD)	3.22	(3.9)	3.42	(4.2)
Age at initial dose, years: mean (SD)	44.5	(17.6)	43.9	(17.3)
Initial dose, mg/day: mean (SD)	80.0 a	(48.9)	82.4 b	(47.8)

N = 1196.

N = 515.

The relaxed LASSO model selected genetically predicted expression of NME1 and AOX1 to build the score associated with discontinuation of azathioprine due to myelotoxicity among White patients. We used the coefficients from a logistic regression with these predicted expressions to build our risk score; it ranged from 0.0051253 to 0.0987391. The score for tertile 1 ranged from 0.0051253 to 0.0263788, tertile 2 ranged from 0.0263968 to 0.050482, and tertile 3 ranged from 0.0505155 to 0.0987391. Based on score tertiles, the incidence‐rate for each tertile was as follows: tertile 1 = 0.94/100 person‐years, tertile 2 = 0.67/100 person‐years, and tertile 3 = 2.00/100 person‐years. The patients in the discovery cohort with scores in the highest tertile had a higher risk of discontinuing azathioprine due to myelotoxicity compared to those in the lowest tertile (hazard ratio [HR] = 2.15, 95% confidence interval [CI] = 1.11–4.19, p = 0.024). These results remained significant after adjusting for a propensity score that included sex, calendar year tertile of initial dose, initial dose, age at baseline, indication, prior TPMT testing, and the first 10 principal components of the genetic data (HR = 2.11, 95% CI = 1.08–4.13, p = 0.030; Table 2).

TABLE 2

Risk of discontinuing azathioprine due to myelotoxicity by cohort based on tertile risk score

	At risk	Events	HR	95% CI	p value
Discovery cohort
First tertile	451	13	Reference
Unadjusted second tertile	351	8	0.75	[0.31–1.80]	0.515
Unadjusted 3rd tertile	399	26	2.15	[1.11–4.19]	0.024
PS a adjusted third tertile	397	26	2.11	[1.08–4.13]	0.030
Validation cohort
First tertile	200	17	Reference
Unadjusted second tertile	158	18	1.23	[0.63–2.39]	0.539
Unadjusted thirrd tertile	159	28	2.00	[1.09–3.65]	0.024
PS a adjusted third tertile	157	27	2.07	[1.10–3.90]	0.024

Abbreviations: CI, confidence interval; HR, hazard ratio; PS, propensity score.

PS included sex, calendar year of initial dose date, initial dose, age a baseline, indication, prior TPMT testing and the first 10 principal components.

Based on the coefficients and same score thresholds used in the discovery cohort, a similar association was observed among the 517 patients in the validation cohort, which included EHR reported non‐White patients and patients with noninflammatory indications (Table 1). There were 63 events, and the incidence‐rate for each tertile was as follows: tertile 1 = 2.57/100 person‐years, tertile 2 = 3.13/100 person‐years, and tertile 3 = 5.22/100 person‐years. Patients with scores within the highest tertile had a higher risk for discontinuing azathioprine due to myelotoxicity (HR = 2.00, 95% CI = 1.09–3.65, p = 0.024), and these results remained significant after adjustment for propensity score (HR = 2.07, 95% CI = 1.10–3.90, p = 0.024; Table 2).

Sensitivity analysis

As a sensitivity analysis we reran the adjusted analysis with last azathioprine dose rather than initial dose in the propensity score. Results in the discovery cohort (HR = 2.10, 95% CI = 1.07–4.11, p = 0.031) and the validation cohort (HR = 2.06, 95% CI = 1.11–3.85, p = 0.023) remained significant.

DISCUSSION

Our results showed that among TPMT and NUDT15 normal metabolizers, a score combining the genetically predicted gene expression of AOX1 and NME1 was associated with discontinuing azathioprine due to bone marrow toxicity.

The metabolism of azathioprine is complex and involves multiple enzymes and transporters. Whereas TPMT and NUDT15 are included in current clinical guidelines, , most patients who discontinue azathioprine due to myelotoxicity are normal TPMT and NUDT15 metabolizers, underscoring the need to examine other genes encoding enzymes involved in the thiopurine pathway. This novel proof‐of‐concept study indicates that the combined genetically predicted expression of NME1 and AOX1 among TPMT and NUDT15 normal metabolizers is associated with higher risk for discontinuation of azathioprine due to myelotoxicity.

The link between increased predicted expression of NME1 and increased risk of azathioprine myelotoxicity is consistent with the known function of the enzyme encoded by the gene. This enzyme catalyzes the biochemical reaction from thioguanosine diphosphate (TGDP) to thioguanosine triphosphate (TGTP). One of the cytotoxic effects of thiopurine drugs is through the incorporation of TGTP into RNA; thus, it is biologically plausible that increased expression of NME1 leads to increased cytotoxicity.

The association between increased expression of AOX1 and increased azathioprine myelotoxicity is more intriguing. AOX1 encodes aldehyde dehydrogenase‐1 (AO), an enzyme that metabolizes 6‐thioguanine (6TG) into 8‐hydroxy‐thioguanine (8‐OH‐TG). , Although the role of AO metabolites is currently unknown, a prior study suggested that they have immunosuppressive activity.

This study has some limitations. First, the GTEx models were developed in a population primarily comprised of individuals of European ancestry ; however, we were able to validate a score initially developed among White patients taking azathioprine for inflammatory conditions in a convenience cohort comprised of White patients taking azathioprine for other indications and patients for whom the EHR indicated non‐White race. The findings suggest that although GTEx was primarily developed using genetic information from individuals of European descent, PrediXcan‐derived results may nonetheless prove relevant for non‐Europeans. Second, the number of events and the discovery cohort were small. Third, our analysis was limited to the prespecified predicted expression of liver tissue, which we anticipated most relevant to azathioprine’s metabolism. Additional insights may be garnered by assessing the association of additional predicted tissue expression (e.g., whole blood). Finally, this study was limited to azathioprine. Future research with other thiopurines that share these genes in their metabolic pathway (e.g., mercaptopurine) may provide additional insight into the association of these genes with discontinuation attributed to myelotoxicity. Thus, further studies are needed to replicate these findings in different cohorts and test whether a polygenic risk score that includes genetic variants from unbiased approaches would further improve prediction or discrimination models. Nevertheless, this study provides a proof‐of‐concept that other genes encoding enzymes involved in the thiopurine metabolic pathway could enhance the prediction of azathioprine myelotoxicity, particularly among TPMT and NUDT15 normal metabolizers.

CONFLICT OF INTEREST

The authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

L.L.D. and C.P.C. wrote the manuscript. L.L.D., A.L.D., T.W.M., P.S.S., W.W., D.P., W.D.D., G.L., P.A., T.S.R., K.A.B., V.K.K., A.M.H., N.J.C., Q.F., C.M.S., and C.P.C. designed the research. L.L.D., A.L.D., P.A., T.S.R., and C.P.C. performed the research. L.L.D., A.L.D., J.T.Z., D.P., W.D.D., P.A., T.S.R., and C.P.C. analyzed the data. L.L.D., J.T.Z., T.W.M., P.S.S., W.W., G.L., Q.F., and C.P.C. contributed analytical tools.

Table S1

Click here for additional data file.