Gene Expression and Thiopurine Metabolite Profiling in Inflammatory Bowel Disease – Novel Clues to Drug Targets and Disease Mechanisms?

Background and Aims Thiopurines are effective to induce and maintain remission in inflammatory bowel disease (IBD). The methyl thioinosine monophosphate (meTIMP)/6-thioguanine nucleotide (6-TGN) concentration ratio has been associated with drug efficacy. Here we explored the molecular basis of differences in metabolite profiles and in relation to disease activity. Methods Transcriptional profiles in blood samples from an exploratory IBD-patient cohort (n = 21) with a normal thiopurine S-methyltransferase phenotype and meTIMP/6-TGN ratios >20, 10.0–14.0 and ≤4, respectively, were assessed by hybridization to microarrays. Results were further evaluated with RT qPCR in an expanded patient cohort (n = 54). Additionally, 30 purine/thiopurine related genes were analysed separately. Results Among 17 genes identified by microarray-screening, there were none with a known relationship to pathways of purines/thiopurines. For nine of them a correlation between expression level and the concentration of meTIMP, 6-TGN and/or the meTIMP/6-TGN ratio was confirmed in the expanded cohort. Nine of the purine/thiopurine related genes were identified in the expanded cohort to correlate with meTIMP, 6-TGN and/or the meTIMP/6-TGN ratio. However, only small differences in gene expression levels were noticed over the three different metabolite profiles. The expression levels of four genes identified by microarray screening (PLCB2, HVCN1, CTSS, and DEF8) and one purine/thiopurine related gene (NME6) correlated significantly with the clinical activity of Crohn’s disease. Additionally, 16 of the genes from the expanded patient cohort interacted in networks with candidate IBD susceptibility genes. Conclusions Seventeen of the 18 genes which correlated with thiopurine metabolite levels also correlated with disease activity or participated in networks with candidate IBD susceptibility genes involved in processes such as purine metabolism, cytokine signaling, and functioning of invariant natural killer T cells, T cells and B cells. Therefore, we conclude that the identified genes to a large extent are related to drug targets and disease mechanisms of IBD.


Introduction
Ulcerative colitis (UC) and Crohn's disease (CD) are chronic remitting and progressive inflammatory bowel diseases (IBD). Primarily, 5-aminosalicylic acid (5-ASA) and glucocorticosteroids are used in the treatment. Glucocorticosteroid dependent or refractory patients are eligible for immunomodulatory therapy with the purine analogues azathioprine or 6-mercaptopurine, methotrexate and/or anti-TNF-antibodies [1,2,3]. Azathioprine and 6-mercaptopurine are prodrugs, converted in vivo to active metabolites via a complex metabolism [4] ( Figure S1). Two main metabolite groups are produced; the phosporylated thioguanine nucleotides (6-TGNs) which comprise thioguanosine mono-, diand triphosphates, and methylated thioinosine phosphates measured as meTIMP [5]. Both metabolite groups contribute to the immunomodulatory effects in different ways [4,6,7,8,9]. Up to 30% of IBD-patients discontinue thiopurine therapy due to adverse events or refractoriness [1,10,11]. A cut off at .230-260 pmol 6-TGN/8610 8 red blood cells (RBC) has been proposed as a lower limit for clinical efficacy [12], but controversy exists regarding its utility, since an important overlap exists between patients in remission and those with active disease. High concentration of the other major metabolite, meTIMP, has mainly been associated with hepatotoxicity [13] and myelotoxicity [11]. Thus, a high meTIMP/6-TGN concentration ratio indicates an increased risk of both therapy failure and adverse events [14,15,16].
The xanthine oxidase inhibitor allopurinol in combination with a reduced dose of thiopurine (,25-33% of original dose) may be considered in patients with this unfavourable metabolite profile. The combination therapy switches the metabolism towards enhanced 6-TGN production and has been shown safe and effective in IBD [15,17]. It also reduces glucocorticosteroid requirements and hepatic as well as nonhepatic side effects [15,17,18]. The use of 5-ASA has also been suggested as an alternative to manipulate the metabolite profile by inhibition of TPMT [19,20,21]. However, the effect of 5-ASA on 6-TGN in vivo varies between studies [11,13,20,22,23,24,25,26] and it does not seem to change the concentration of meTIMP [23,24]. The use of 5-ASA to modulate the metabolite profile has not been implemented in clinical practice in a similar way as allopurinol.
The underlying mechanism to why a proportion of patients preferentially metabolize AZA and 6-MP to meTIMP is currently unknown. Interindividual variation in metabolite profiles and drug response may be explained by differences in activities of drugmetabolizing enzymes, and their correlation with corresponding transcript or protein levels. A well-known cause of adverse reactions to thiopurines is reduced or absent thiopurine Smethyltransferase (TPMT) activity, where patients with low TPMT activity accumulate myelotoxic concentrations of 6-TGNs if treated with standard doses [27]. However, not all interindividual differences in thiopurine metabolism and response can be attributed to variations in TPMT [14,15,28,29,30], since a large number of enzymes are involved.
Blockage of inosine 59-monophosphate dehydrogenase (IMPDH) activity could, based on its position in the metabolic pathway of thiopurines ( Figure S1), restrict the formation of 6-TGNs and therefore explain a high meTIMP/6-TGN concentration ratio. Indeed, the IMPDH activity was inversely correlated with the concentration of meTIMP in our previous work. However, no correlation with the concentration of 6-TGN was observed [25,31].
Inosine triphosphatase (ITPase) is involved in a metabolic loop where 6-TIMP is reconverted from 6-thioinosine triphosphate ( Figure S1). ITPase deficiency may contribute to the metabolite profile by increasing the concentration of methylated metabolites (methyl thioinosine triphosphate) [32]. Furthermore, most nucleoside analogues enter and exit cells via nucleoside transporters [33]. Impaired function, as well as up-or down-regulation of transport proteins with different specificities for thiopurine metabolites probably affect the metabolite profiles, as may variations in activities of intracellular nucleotidases and kinases.
The aim of this study was to explore the molecular basis of differences in metabolite profiles. We performed a whole genome expression analysis in blood samples from thiopurine treated patients with IBD and related the findings to metabolite concentrations and clinical patient characteristics.

Ethics Statement
The study was approved by the Ethics Committee at Linköping University, Sweden, dnr 03-260 and M58-06. Written informed consent was obtained from all patients.

Patients
An explorative cohort of IBD-patients (n = 21) with normal TPMT activity (.8.9 U/mL pRBC; units per mL packed RBC) was selected based on differences in their metabolite profiles. We defined a high metabolite concentration ratio as meTIMP/6-TGN .20 based on metabolite determinations in our laboratory and on studies of deviant metabolism by others [14,15,34]. In our experience, this ratio corresponds to the 75 th percentile of the meTIMP/6-TGN concentration ratios (median 12) amongst 1220 patients with IBD on thiopurine therapy and TPMT activity in the normal range. Ten patients with meTIMP/6-TGN concentration ratios .20 (R20) were included in a microarray analysis, as were four patients with metabolite ratios corresponding to the median metabolite ratio (range 10.0-14.0, Median) and seven patients displaying a profile with a metabolite ratio #4 and 6-TGN $100 pmol/8610 8 RBC (R4). The cut-off $100 pmol 6-TGN/ 8610 8 RBC was employed to ensure acceptable analytical precision.
All patients had been on long term medication with an unchanged thiopurine dose for at least 2 weeks prior to blood sampling. Metabolite profiles were stable as judged from historical records, with at least two observations of the same kind available in 19/21 patients. Patients who had received blood transfusion within 4 months were not included. Disease activity (score $5 in active disease) was noted based on the Walmsley index for UC (n = 10) and the Harvey-Bradshaw index for CD (n = 10) and patient characteristics were noted (Table S1).
The data from the microarray analysis was further validated by reverse transcription quantitative PCR (RT qPCR) in an expanded patient cohort (n = 54), including the initial study population with the exception of one sample with insufficient amount of RNA (Table S2).

Isolation of RNA from Peripheral Blood Samples
Blood samples were collected in PreAnalytiX Paxgene TM blood RNA tubes (Becton Dickinson, Franklin Lakes, NJ). RNA was isolated using the PreAnalytiX Paxgene TM blood RNA kit (Qiagen, Hilden, Germany) on the day of sampling, according to the manufacturers instructions. RNA concentration was assessed with NanodropH ND-1000 spectrophotometer (Nanodrop Technologies, Wilmington, DE) and RNA integrity with 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA).

Microarray Analysis
RNA samples were treated with the SnX TM globin depletion reagent [37] and analysed by AROS Applied Biotechnology (Aarhus, Denmark) with the Affymetrix GeneChip Human Genome U133 Plus 2.0 array (Affymetrix, Santa Clara, CA), representing the entire human genome with more than 38 500 well characterized genes.

RT qPCR
Real-time PCR was performed with the FAST 7500 real-time PCR system and reagents from Applied Biosystems (Foster City, CA) with 5-10 ng cDNA per reaction in a final volume of 10 mL. Thirty-two potential reference genes (TaqManH Express Human Endogenous Control Fast Plate, Applied Biosystems) were evaluated for low sample-to-sample variation using the Normfinder algorithm [38] and cDNA from six patients. The mRNA expression [C T (threshold cycle)] of each target gene (Table S3) was normalized against the expression level of the selected reference genes (GUSB, YWHAZ, and MRPL19) with Genex Professional software version 4.3.8 (MultiD Analysis AB, Göteborg, Sweden) to obtain a delta-C T (dC T ). The relative expression was determined for each gene in relation to the sample with the lowest expression (highest C T ).
By exploring the Pharmacogenomics knowledge base (http:// www.pharmgkb.org/), the KEGG pathway of purine metabolism (http://www.genome.jp/kegg/pathway.html#nucleotide) and the thiopurine literature, 30 genes with a proven or potential relationship to the mechanisms, metabolism or transmembrane transport of purine/thiopurines were selected and analysed separately by RT qPCR (Table S3).

Data Analysis
Analysis of microarray data. The image files (cel format) were imported to the GeneSpring GX 11 software (Agilent Technologies, Santa Clara, CA). Data was background corrected and normalized by the robust multiarray analysis algorithm [39].
In order to identify new candidate genes differently expressed over metabolite profiles, low-intensity signals were removed by a stringent filtering, retaining all data with a signal intensity greater than the 75 th percentile of all intensities. Thereafter, only intensities above the 20 th percentile in 100% of at least 1 of the 3 metabolite profiles were retained, leaving 7325 probe sets for further analysis. Genes differently expressed between metabolite profiles were identified by analysis of variance (ANOVA) without correction for multiple testing and with the Student-Newman-Keuls post hoc test at a level of statistical significance (P-value) of 0.001.
The 7325 probe sets were further included in a Spearman rank order correlation analysis against the individual dose-normalized metabolite concentrations and the meTIMP/6-TGN concentration ratio. The results were considered statistically significant if P,0.001.
Pathway analyses. In order to detect differences between metabolite profiles, gene set enrichment analysis (GSEA) was performed on the normalized data set, using the Broad Institute GSEA software, version 3.0 [40,41]. Patients with meTIMP/6-TGN concentration ratio .20 were compared with patients with ratios #4 using the C2 molecular signatures database (MSigDB) of pathways of the Kyoto Encyclopedia of Genes and Genomes (KEGG), containing 186 gene set (414 pathways).
The interactions between gene products identified in this study and genes present at susceptibility loci identified for CD, UC and IBD [42] were evaluated with the Search Tool for the Retrieval of Interacting Genes/Proteins database, STRING, version 8.3 [43]. All prioritized candidate susceptibility genes were extracted from Table S2 of Jostins et al. [42], whereas for susceptibility loci with no prioritized genes, all genes were extracted. Candidate susceptibility genes with erroneous interactions with genes identified in this work were manually removed from the network.
Statistics. Dose-normalized metabolite concentrations (pmol metabolite per mg azathioprine) were used when investigating relationships between gene expression and metabolite concentrations. 6-mercaptopurine doses were converted to azathioprine doses, assuming a conversion factor of 2.08 [44].
Correlations between variables were evaluated using the Spearman rank order correlation coefficient, R s . Median (range) values are given. For group comparisons of continuous variables, the Mann-Whitney U-test, or the Kruskal Wallis test were used. For categorical variables, Fishers exact test was used. Two-sided testing was used and considered statistically significant if P,0.05. Multiple linear regression analyses, using backward stepwise removal or inclusion of variables, were applied to assess the relationship between log transformed dose-normalized metabolite concentrations or meTIMP/6-TGN concentration ratios as dependent variables and gene expression levels expressed as dC T values and the use of concomitant drugs as independent variables. A P to enter of 0.15 and a P to exit the models of 0.15 was adopted. Log transformation was necessary to normalise the distribution of the dependent variables.

Genes Identified by Microarray Screening -ANOVA
Four genes that were differently expressed over metabolite profiles were identified in the explorative patient cohort using the stringent filtering of microarray data (P,0.001); FAM46A, SLX1A, TGOLN2, UBE2A and included in the analyses with RT qPCR (Table S3).

Genes Identified by Microarray Screening -Spearman Rank Order Correlation Analyses with Metabolite Concentrations
Among the significant genes identified by means of Spearman rank order correlation analyses of the 7325 probe sets, the top five or six genes with the most significant correlations with the individual metabolite concentrations, or the meTIMP/6-TGN concentration ratio, were selected for further evaluation.
In total, fourteen candidate genes were identified by Spearman rank correlation analyses, one of which overlapped with the four genes identified by ANOVA (UBE2A). No genes with known association with the metabolic pathways of thiopurine drugs or purines were identified in the microarray screen (employing a pvalue ,0.001).
RT qPCR Data vs. Microarray Data -General Screening Altogether, 17 genes, identified by microarray screening, were taken to RT qPCR in the expanded patient cohort (n = 54) ( Table  S3). The relative gene expression levels of nine genes, correlated with the concentration of meTIMP, 6-TGN or the meTIMP/6-TGN concentration ratio (Table 1, Table S5, and exemplified in Figure 1).
Four genes (CD1D, DEF8, HVCN1, and TUSC2) correlated positively with 6-TGN and all except HVCN1 correlated negatively with the meTIMP/6-TGN concentration ratio. In the microarray data, DEF8 and HVCN1 displayed a significant positive correlation with 6-TGN, whereas the other two genes displayed a significant negative correlation with the meTIMP/6-TGN concentration ratio.
CTSS, FAM156A, GNB4, and PLCB2 were negatively correlated with both the concentration of meTIMP and the meTIMP/6-TGN concentration ratio. In the microarray data, FAM156A and GNB4 displayed a significant inverse correlation with meTIMP, whereas the other two genes displayed a significant inverse correlation with the meTIMP/6-TGN concentration ratio.
LAP3 correlated negatively only with the meTIMP/6-TGN concentration ratio, both in the expanded cohort and in the microarray screen.
Substantial interindividual variation in the gene expression levels was observed, with a considerable overlap in expression levels between the three metabolite profiles ( Table 2). Nevertheless, the gene expression levels of CD1D, CTSS, FAM156A, GNB4, and PLCB2 were lower in patients with meTIMP/6-TGN concentration ratios .20 as compared with those with a median metabolite ratio (Median). Of these five differences, four were also noticed when comparing patients with high metabolite ratios (.20) with those with low metabolite ratios (#4) ( Table 2).

Genes Related to the Metabolic Pathway of Thiopurines
Thirty genes potentially associated with the mechanism, metabolism or transmembrane transport of thiopurines or purines, were analysed with RT qPCR (Table S3). When the relative gene expression levels were evaluated, nine genes correlated with the concentration of meTIMP, 6-TGN and/or the meTIMP/6-TGN concentration ratio (Table 3, Table S6, and exemplified in Figure 2).
A positive correlation was observed between the gene expression level RAC1 and the concentration of 6-TGN, whereas HPRT1 correlated negatively with this metabolite (both P = 0.04).
Two genes, XDH and NT5C1B, were not expressed in peripheral blood as judged by the RT qPCR data.
Patients with meTIMP/6-TGN ratios .20 showed lower expression levels of MGST2, TPMT and NME6, but higher expression levels of IMPDH2 and NT5E than patients with metabolite ratios #4. However, a large overlap in gene expression levels over the three metabolite profiles was observed ( Table 4).
The TPMT activity in RBC did not correlate with the mRNA expression of TPMT (R s = 0.10, P = 0.45). The GSEA showed no significant enrichment of genes of any pathway in any phenotype (different metabolite profiles) when a false discovery rate of 0.25 and 1000 permutations of each phenotype was applied (data not shown). However, RAC1, PLCB2 and GNB4 overlapped with the KEGG chemokine signaling pathway (P = 4.6610 25 ) and RAC1 and PLCB2 with the KEGG Wnt signaling pathway (P = 0.002).
Sixteen out of 18 genes (89%) which correlated significantly with meTIMP and/or 6-TGN or the meTIMP/6-TGN concentration ratio in the expanded patient cohort (n = 54) interacted in networks with 41 genes present at 6, 4 or 21 susceptibility loci identified for CD, UC or IBD [42], respectively, as judged by the STRING analysis (Figure 3; FAM156A and HVCN1 were not present in any network). All the network-identified candidate susceptibility genes belonged to genes prioritized by Jostins et al. [42] with the exception of four genes at two susceptibility loci without prioritized genes (CTH, ADK, CAMK2G and PLAU). Genes selected for a potential association with the metabolism (IMPDH2, HPRT1, TPMT, NT5E, ENTPD1 and NME6) or transmembrane transport (SLC29A2) of thiopurines or purines connected with candidate susceptibility genes linked to the KEGG pathway of purine metabolism, except MGST2 which linked to candidate susceptibility genes involved in the KEGG pathway of glutathione metabolism ( Figure 3, and data not shown). Genes identified by microarray screening were linked to the following KEGG pathways through interactions with candidate susceptibility genes: an inter-pathway connection between 'cysteine and methionine metabolism' and 'glutathione metabolism' (LAP3), inositol phosphate metabolism and signaling through phosphatidylinositol, calcium, chemokine or Wnt (PLCB2), chemokine and Wnt signaling (RAC1), chemokine signaling (GNB4) (Figure 3, and data not shown). For four of the genes (CD1D, CTSS, TUSC2, and DEF8), the STRING analysis identified interactions solely based on textmining and indicated involvement in functioning of invariant natural killer T cells, T cells and B cells, cytokine and Genes that showed a significant correlation (P,0.05) between relative gene expression levels (RT qPCR data) and the concentration of 6-TGN, meTIMP or the meTIMP/ 6-TGN concentration ratio in the expanded patient cohort (n = 54). b Significant (P,0.001) positive (P) or negative (N) correlations observed in the microarray data of the explorative patient cohort (n = 21) are indicated. doi:10.1371/journal.pone.0056989.t001 chemokine signaling, anti-proliferative and apoptotic effects, and inositol phosphate metabolism ( Figure 3, and data not shown).

Active Disease vs. Remission
The gene expression levels of CTSS, DEF8, HVCN1, NME6, and PLCB2 were higher in remission than in active disease (P#0.04), whereas a decreased expression of IMPDH2 was associated with remission (P = 0.008). Eight of 9 patients with active disease had CD. Considering CD patients only, the results were essentially the same with inverse relationships between the gene expression levels of CTSS (R s = 20. The measured concentrations of meTIMP and 6-TGN were no different in patients in remission (n = 43) compared with those with active disease (n = 9, P$0.13).

Regression Analyses
All genes that displayed an individual correlation with one of the metabolites or the meTIMP/6-TGN concentration ratio in the expanded patient cohort, were, together with concomitant therapy with 5-ASA or glucocorticosteroids (distribution between metabolite profiles, P = 0.11 and 0.006, respectively) assessed using multiple linear regression analyses.
The models were essentially the same using untransformed data (data not shown).

5-ASA and Glucocorticosteroid Therapy
The gene expression level of TPMT was higher among patients   Genes that showed a significant correlation (P,0.05) between relative gene expression levels (RT qPCR data) and the concentration of 6-TGN, meTIMP or the meTIMP/6-TGN concentration ratio in the expanded patient cohort (n = 54). b None of the genes selected by means of their potential relationship with the mechanism, metabolism or transport of purines/thiopurines fulfilled the established significance level (P,0.001) in the microarray data. doi:10.1371/journal.pone.0056989.t003  However, there was a large overlap in gene expression levels between the two groups (data not shown).

Discussion
An aberrant thiopurine metabolism with preferential formation of the methylated metabolites has been related to lesser efficacy and an increased risk of adverse events [11,14,15,16]. At present it is unknown why up to 25% of patients display such a metabolite profile. In order to explore the molecular basis of differences in thiopurine metabolite concentrations we adopted a broad approach using whole genome expression analysis in blood samples from patients with distinct metabolite profiles.
Based on the microarray screening we selected seventeen candidate genes for further analysis using RT qPCR in the expanded patient cohort, and nine of these genes demonstrated significant correlations with meTIMP and/or 6-TGN concentrations. However, among the genes there was none with a known or suspected relationship to the purine/thiopurine metabolism, transport or drug effects. Furthermore, although large interindividual differences in gene expression levels were observed (Table 2  and Table 4), only small differences were noticed between the three metabolite profiles (R20, Median and R4). Thus, based on gene expression, we did not find a strong influence on the three distinct thiopurine metabolite profiles.
Thirty genes potentially associated with the mechanism, metabolism or transmembrane transport of thiopurines or purines were explored by means of RT qPCR in the expanded patient cohort (irrespective of the microarray result). The expression levels of nine genes correlated significantly with meTIMP and/or 6-TGN concentrations. All of these genes, except RAC1, putatively affect metabolism or transport of thiopurines or purines. RAC1 encodes a small GTPase and is involved in the signaling from the T-cell receptor in activated CD4+ cells. Its downstream targets promote cellular survival. The 6-TGN metabolite 6-TGTP has the potential to induce apoptosis in these cells by blocking the Rac1 protein [7]. Rac1 also facilitates the interaction between antigen presenting cells and effector cells, a process disturbed by 6-TGTP [45].
Multiple regression analyses, using gene expression levels and concomitant therapy with 5-ASA or corticosteroids as predictors, explained at most 46% of the variation in the dependent metabolite variables. Four of the five genes included in the regression models had no established prior relationship with the metabolic scheme of thiopurines, although NT5E has been suggested to facilitate cellular uptake of thiopurine metabolites by extracellular de-phosphorylation [46]. Possibly, the thiopurine metabolite profile depends on several interacting genes, each with an individual small effect. However, it is also possible that the metabolite profile mainly reflects a post-transcriptional regulation. In clinical practise, thiopurine metabolites are measured in RBC, whereas gene expression levels were measured in the target cells of therapy. The use of RBC as a surrogate compartment for the target cells of therapy (the mononuclear cells) may obscure the ability to establish relationships between gene expression levels and metabolite concentrations.