In our drugs' optimization results, Paclitaxel, Sorafenib, and Vincristine were found to have potential therapy effect in the basal-like subtype of breast cancer. Consistent with clinical findings, Paclitaxel and Vincristine were the optimal adjuvant therapy for triple-negative breast cancer [35,36]. Sorafenib is a multiple targeted agent which can inhibit tumor cell proliferation and angiogenesis by inhibiting the activation of multiple different kinases [37], and our results indicate that Sorafenib plays a therapeutic role in the basal-like subtype of breast cancer mainly through affecting specific risk pathway crosstalk mediated by hsa-miR-30a, hsa-miR-222, and hsa-miR-193a. Some studies have confirmed that hsa-miR-30a, hsa-miR-222, and hsa-miR-193a play key roles in breast cancer [38–40]. Irinotecan, an antitumor enzyme inhibitor mainly used for the treatment of colorectal cancer [41], is suitable for the her2-enriched subtype, which mediates the specific crosstalk among the risk pathways of the her2-enriched subtype through regulating hsa-miR-23a and hsa-miR-324. In accordance with the result of WT Kuo and Eissa [42,43], hsa-miR-324 and hsa-miR-23a have distinct biological functions in breast cancer. Vemurafenib has long been approved for the treatment of metastatic melanoma with BRAF mutation [44], and our results showed that this drug had a damaging effect on the specific crosstalk of risk pathway of the luminal A subtype through action on hsa-miR-145. Just as some researches have shown that miR-145 is a potential cancer biomarker and serves as a novel target for cancer therapy, including breast cancer [45]. Vorinostat as an anticancer agent that inhibits histone deacetylases, approved for cutaneous T-cell lymphoma [46], and plays a key role in the epigenetic regulation of gene expression. Vorinostat could act on the specific risk pathways crosstalk of the luminal B subtype via 14 miRNAs (Figure 6), which have been found to play important roles in the occurrence and development of breast cancer, such as hsa-miR-155, hsa-miR-34a, hsa-miR-17, hsa-miR-22, and hsa-miR-140 [47–51].

**Figure 6.** The mechanism of optimal therapeutic drugs in each subtype of breast cancer. Sorafenib, Paclitaxel, and Vincristine were applicable for the basal-like subtype treatment, Irinotecan was optimum for the her2-enriched subtype treatment, Vemurafenib was suitable for the luminal A subtype treatment, and Vorinostat was applied to the luminal B subtype treatment.

#### **4. Discussion**

Breast cancer is a complex disease with high heterogeneity in terms of the underlying molecular alterations, the cellular composition of tumors, and even the clinical outcomes. Different subtypes exhibit distinct biological behavior, prognosis, and usually different responses to drug treatment [52], yet identifing applicable drugs for each subtype still largely remains limited. Therefore, it is urgently needed to develop a systematic pipeline to identify medications for different subtypes of breast cancer.

The occurrence and development of tumors is a complex process involving many steps, links, and factors. It is mostly the action of a single molecule (gene or miRNA) that leads to poor therapeutic effect among many chemotherapeutic regimens [53]. In recent years, many researches have revealed that the occurrence of tumors is closely related to the abnormality of biological pathways, and crosstalk of abnormal pathways is one of the prime reasons for the poor outcomes of tumor treatment [54]. Studies have shown that regulatory molecules such as non-coding RNA participate in the anomaly of biological pathways through the regulation of genes, adding to the difficulty of cancer treatment [55]. In order to actually reflect the intricate crosstalk of pathways, we have developed a new method based on biological pathways—that is, reconstruction of biological pathways which include both genes and miRNAs. We have also identified the optimal drugs by quantifying the effect of candidate drugs on miRNA-mediated crosstalk of pathways. We have successfully identified the specific crosstalk of pathways in each subtype of breast cancer and revealed their pathogenesis respectively by applying this method. Moreover, we also screened applicable drugs for each subtype of breast cancer. We successfully screened the most suitable drugs for each subtype of breast cancer, including Paclitaxel

and Vincristine, which are breast cancer treatment drugs in clinical application. On the basis of the original application, we accurately identified their applications in each subtype, such that Paclitaxel and Vincristine were best for basal-like, Irinotecan was suitable for her2-enriched, and Vorinostat was the optimal drug for luminal B subtypes. We also identified other anticancer drugs application in each subtypes of breast cancer. The results show that our approach could help doctors to further improve treatment strategies with the current menu of chemotherapy options.

Currently, several methods have been proposed to optimize drugs for human cancers. For example, Lamb et al. provided a computational method to connect diseases and their potential therapeutic small molecules based on gene expression profiles form disease and cultured human cells treated with bioactive small molecules respectively [56]. Gottlieb et al. predicted novel drug indications based on multiple drug–drug and disease–disease similarity measures [57]. Furthermore, Malas et al. prioritized drugs using the semantic information between drug and disease concepts [58]. Comparing with these methods, our study has some unique features. First, we considered the role of non-coding RNAs in our approach. Second, our study optimized anticancer drugs by measuring their effects for mediating the crosstalk between risk pathways, which was an important molecular mechanism in the initiation and progression of human cancers. Finally, we optimized candidate drugs for different breast cancer subtypes, which may further promote the precise use of drugs for human cancer.

There are also several limitations in our study. First of all, drugs targeting miRNAs for therapeutic purposes are limited, and there are many drugs without miRNA targets. Secondly, miRNAs affected by the drugs are required for further study. We believe that more and more drugs that regulate miRNAs and drug-regulated miRNAs will be discovered with the development of in-depth study on the interaction of drugs and miRNAs, and our method can identify the optimal therapeutic agent for complex diseases more accurately and comprehensively. In summary, the results in this study highlight that dissecting subtype-specific risk pathway crosstalk could provide novel insights into the underlying molecular mechanisms and thus promote the drug discovery for various breast cancer subtype. Moreover, we focused on breast cancer in this study, but the method proposed here could also be applied to many other complex diseases, as pathway crosstalk is widespread in biological systems and the dysregulation of which play a critical role in the occurrence of disease.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/10/9/657/s1, Table S1: Detailed information of the 36 anticancer drugs and their targets.

**Author Contributions:** Conceptualization, Y.X. and Y.X.; data curation, J.W. and C.Z.; formal analysis, Y.X., H.Z., and C.H.; investigation, S.L., H.Z., and J.W.; methodology, Y.X.; supervision, Y.X. and L.W.; validation, H.Z. and Q.W.; visualization, S.L. and C.H.; writing—original draft, Y.X. and S.L; writing—review and editing, C.Z., Q.D., D.S., Y.X., and L.W.

**Funding:** This work was supported by the Fundamental Research Funds for the Provincial Universities (2017JCZX51,2017JCZX53, 2017JCZX54); the Heilongjiang Postdoctoral Science Foundation (LBHZ16123, LBH-Z17218, LBH-Z17110); the China Postdoctoral Science Foundation (2018M631943, 2016M600260); the China Postdoctoral Science Special Foundation (2019T120280); the National Natural Science Foundation of China (61903105, 31801115, 31801107, 61803129); and the Wu lien-teh Youth Science Fund Project of Harbin Medical University (WLD-QN1407)

**Acknowledgments:** Not applicable.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Impact of** *CYP2C9* **and** *VKORC1* **Polymorphisms on Warfarin Sensitivity and Responsiveness in Jordanian Cardiovascular Patients during the Initiation Therapy**

#### **Laith N. AL-Eitan 1,2,\*, Ayah Y. Almasri 1 and Rame H. Khasawneh 3**

Received: 15 October 2018; Accepted: 21 November 2018; Published: 27 November 2018 -

**Abstract:** Warfarin is an oral anticoagulant frequently used in the treatment of different cardiovascular diseases. Genetic polymorphisms in the *CYP2C9* and *VKORC1* genes have produced variants with altered catalytic properties. A total of 212 cardiovascular patients were genotyped for 17 Single Nucleotide Polymorphisms (SNPs) within the *CYP2C9* and *VKORC1* genes. This study confirmed a genetic association of the *CYP2C9*\*3 and *VKORC1* rs10871454, rs8050894, rs9934438, and rs17708472 SNPs with warfarin sensitivity. This study also found an association between *CYP2C9* and *VKORC1* genetic haplotype blocks and warfarin sensitivity. The initial warfarin dose was significantly related to the *CYP2C9*\*3 polymorphism and the four *VKORC1* SNPs (*p* < 0.001). There were significant associations between rs4086116 SNP and TAT haplotype within *CYP2C9* gene and rs17708472 SNP and CCGG haplotype within *VKORC1* gene and warfarin responsiveness. However, possessing a *VKORC1* variant allele was found to affect the international normalized ratio (INR) outcomes during initiation of warfarin therapy. In contrast, there was a loose association between the *CYP2C9* variant and INR measurements. These findings can enhance the current understanding of the great variability in response to warfarin treatment in Arabs.

**Keywords:** *CYP2C9*; *VKORC1*; warfarin; warfarin initiation phase of therapy; INR; pharmacogenetics study

#### **1. Introduction**

Warfarin is a commonly prescribed oral anticoagulant that is employed for the treatment of venous and arterial thromboembolic disorders and cardiac valve replacements [1]. However, interindividual genetic variation causes great variability in dosage requirements, making the latter a problematic issue for physicians. A higher or lower dose than needed could lead to bleeding and thrombotic risk, respectively [2,3]. Two-thirds of warfarin dose variation was due to environmental factors like age, body mass index, smoking status, gender, and diet, among others, while the remaining one-third is caused by genetic factors such as the *CYP2C9* and *VKORC1* genes [4–6].

Belonging to the cytochrome P450 superfamily, the *CYP2C9* gene is involved in the metabolism and clearance of S-warfarin, the latter of which is a racemic form of warfarin together with R-warfarin [7,8]. *CYP2C9* is located on the long arm of chromosome 10, and like other members of *CYP2C*, *CYP2C9* is highly polymorphic [9–11]. Although there are over 50 single nucleotide polymorphisms (SNPs) located in the regulatory and coding region of the *CYP2C9* sequence, the most studied *CYP2C9* polymorphisms are *CYP2C9\*2* (R144C) and *CYP2C9\*3* (I359L) [12,13]. It was found

that individuals carrying the *CYP2C9* \*2 \*3 alleles reduced the elimination of S-warfarin and, therefore, plasma concentrations of the latter increased significantly compared to the wild-type allele and the individuals with these variant alleles need a lower warfarin maintenance dose [14].

Similarly, the *VKORC1* gene encodes the vitamin K epoxide reductase complex subunit 1 and it is a warfarin target [5]. The vitamin K epoxide reductase complex subunit 1 normally catalyzes the carboxylation reaction of the vitamin K-dependent protein glutamic acid residues in order to activate it, the latter of which are responsible for catalyzing the clotting factor pathway [15]. Several genetic studies conducted in different populations suggest that the G3673A (rs9923231), C6484T (rs9934438), and G9041A (rs7294) polymorphisms of the *VKORC1* gene are the most common and well-studied [7,16].

Analysis of *CYP2C9* and *VKORC1* gene polymorphisms revealed that they were responsible for 10% and 40% of warfarin dose requirement variance, respectively [17]. Combined with clinical data, both of the aforementioned genes can explain up to 60% of the warfarin variance [18,19]. More than 50 years ago, the dose of warfarin was determined by trial and error, with an initial dose (2–10 mg/day) dependent on the indication of warfarin and clinical factors, regardless of the effect of the genetic factor [20]. Patients are treated with warfarin in two stages, the first called the initiation phase of treatment, which is considered as the stage in which the initial international normalized ratio (INR) value of the patient is unstable and fluctuates up and down. With the second (maintenance) stage of the therapy, the patient is within the therapeutic INR for at least two consecutive visits [21]. However, pharmacogenetics for specific populations including Jordanian Arabs is necessary. Therefore, the objective of this study was to recognize genetic variations within the *CYP2C9* and *VKORC1* genes that are involved in warfarin sensitivity and responsiveness in Jordanian cardiovascular patients of Arab descent during initiation of treatment.

#### **2. Materials and Methods**

#### *2.1. Patient Population and Study Design*

The study population consisted of 212 unrelated warfarin intake patients selected from the Jordanian-Arab population, from the Anticoagulation Clinic at the Queen Alia Heart Institute (QAHI) in Amman, Jordan. Informed consent was obtained from all subjects. The study protocol was approved by the Human Ethics Committee at Jordan University of Science and Technology in Irbid, Jordan, and the Royal Medical Services in Amman, Jordan. Ethical approval code: 13/78/2014.

In this study, patients have cardiovascular diseases and are prescribed warfarin as an anticoagulant therapy. Inclusion criteria involved patients being 18 years or older and having received warfarin for at least three months. Patients who did not provide informed written consent, did not visit the anticoagulation clinic regularly, took *CYP2C9* inducers or inhibitors, or did not have a complete data set were excluded.

Initially, 350 patients were screened and, based on the aforementioned inclusion and exclusion criteria, 300 patients were approached to participate in this study (Figure 1). 80 patients were subsequently excluded because of inability to complete the treatment program or refusal. Of the remaining patients, 220 accepted to be part of this study, after which an additional eight patients were excluded from the final analysis because of a failure in genotyping. In total, whole sets of data were obtained from 212 patients with cardiovascular disease who were being treated with warfarin. Data was collected on demographic (age, gender, and body mass index) and lifestyle (smoking status and diet) characteristics as well as medical history (diseases and clinical features of warfarin therapy). Clinical features included the target INR, mean weekly warfarin doses required to reach the target INR, and the use of concomitant medication, all of which are required to be known in order to adjust warfarin doses. All data were blinded and obtained through semi-standardized interviews and medical records.

**Figure 1.** Flow chart depicting study design. INR: international normalized ratio.

All subjects in this study received warfarin anticoagulant therapy according to RMS anticoagulation protocol, which began with 2.5 to 10 mg nightly doses. INR monitoring is required at least once a week for the first three to four weeks after the initiation of therapy. After three consecutive visits, a patient with stable INR reaches the maintenance dose and is monitored for warfarin treatment administration by clinic protocol.

#### *2.2. Outcome Measurement*

Oral anticoagulant therapy was mandated by the prothrombin time (PT) that is, evaluated using an automated method over STAGO coagulometric unit in the QAHI laboratory. To calculate INR, there was a blood coagulation (clotting) test. INR monitoring for dose adjustments was determined by the physician and pharmacist. INR values between January 2014 and November 2015 were obtained from medical records, and these values were then used to divide the patients according to warfarin responsiveness into: (1) good responders, with an INR value within the target range; (2) poor responders, with an INR value under the target range; and (3) extensive responders, with an INR above the target range.

Further, based on their warfarin sensitivity, patients were divided into resistant, normal, or sensitive to warfarin as follows: (1) warfarin resistance (or Poor metabolizer), largest daily doses were required to keep a patient's INR within a therapeutic range (dose required > 49 mg/week); (2) warfarin normal patient (or Intermediate metabolizer), intermediate doses required (doses between 21–49 mg/week); and (3) warfarin sensitive patient (or Extensive metabolizer) lowest doses required (required dose <21 mg/week) [22].

#### *2.3. SNP Selection, DNA Extraction, and Genotyping*

In this study, 17 SNPs in the *CYP2C9* and *VKORC1* genes were selected from public databases and genotypes. Information about the aforementioned SNPs is shown in Table S1. Genomic DNA was extracted within one week of blood collection using the commercially available Wizard Genomic DNA Purification Kit (Promega Corporation, Madison, WI, USA) according to the manufacturer's instructions. After extraction, the DNA was diluted in 96-well plates using an automated robotic system to achieve concentrations of 20 ng/μL (50–500 μL). Concentrations were confirmed with the Nano-Drop ND-100 (Thermo Scientific, Wilmington, DE, USA). Genotyping was carried out by means of the MassARRAY®system (iPLEX GOLD) (Sequenom, San Diego, CA, USA), which was carried out at the Australian Genome Research Facility (AGRF) (Sequenom).

#### *2.4. Statistical Analysis*

Discrepancy and call rates were calculated using Microsoft Excel, and the deviation from Hardy Weinberg Equilibrium (HWE) was assessed using the Pearson X2 test. Minor allele frequencies (MAF) and HWE *p*-values for genotypic distribution were calculated via the Court lab-HW calculator. To test which of the chosen SNPs is associated with warfarin response, various statistical genetic association analyses were conducted, such as the chi-square, nonparametric correlation tests (Kruskal-Wallis and Tukey Pairwise comparison) and haplotype genetic analysis test. The Statistical Package for the Social Sciences (SPSS) version 21.0 and the SNPStat Web Tool (https://www.snpstats.net/start.htm) were used to perform all analyses.

#### **3. Results**

#### *3.1. Study Group*

The study group comprised of 212 unrelated Jordanian-Arabs patients treated with warfarin, with a mean age (±SD) of 56.03 (±17.68) years, a median age of 60, and an age range of 18 to 85 years. There were 34 poor responders (16%), 146 moderate responders (68.9%), and 32 extensive responders (15.1%). Table 1 summarizes the demographic, lifestyle, and medical characteristics of each of the three groups.

In total, 17 SNPs (100%) passed the quality control measures for throughput genotyping and were analyzed by the MassARRAY®system (iPLEX GOLD) with high accuracy and a 97% average success rate. The genotypic discrepancy average (±SD) rate over the 17 loci was only 0.06% (±0.0004%) out of the entire cohort (212 subjects). Genotypic and allelic frequencies are shown in Table S2.

For the 17 SNPs examined in this study, all were in accordance with the HWE. Ten polymorphisms (rs104894539, rs104894540, rs104894541, rs104894542, and rs61742245 in *VKORC1*, and rs28371685, rs28371686, rs72558191, rs9332131, and rs9332239 in *CYP2C9*) were non-polymorphic. In contrast, the seven remaining SNPs (rs10871454, rs8050894, rs9934438, and rs17708472 located in *VKORC1*, and rs1799853, rs4086116, and rs1057910 located in *CYP2C9*) were polymorphic and thus included in the study. The minor alleles and their frequencies for the successful genotyped SNPs are shown in Table S3.

#### *3.2. Effect of* CYP2C9 *and* VKORC1 *Polymorphisms on Warfarin Sensitivity during Initiation Phase of Therapy*

Regarding the association of *VKORC1* and *CYP2C9* SNPs with warfarin sensitivity among the three inclusion groups, significant differences in proportions among genotypes were observed at all tested *VKORC1* SNPs (*p* < 0.001) (Table 2). Significant differences were also observed between two SNPs of the *CYP2C9* gene (rs4086116 (*p* = 0.012) and rs1057910 (*p* < 0.001)), as shown in Table 2. Moreover, there was a significant association observed between *VKORC1* and *CYP2C9* haplotypes and warfarin sensitivity (*p* < 0.0001) (Table 3).

**Table 1.** Descriptive analysis of demographics and clinical characteristics of 212 cardiovascular patients treated with warfarin at the Queen Alia Heart Institute.

a Mean Standard deviation in square brackets. b CHD: Chronic heart disease. c MVR:Mitral valve replacement.

d AVR: Aortic valve replacement. e AF: Atrial Fibrillation. f DVR: Double valve replacement.

\* Chi-Square Test with *p*-value < 0.05 is considered significant.

**Table 3.** Frequencies of the haplotypes of *VKORC1* and *CYP2C9* genes among the 212 warfarin sensitive patients.

\* Genetic haplotype frequency of 212 warfarin intake patients, \*\* *p*-value < 0.05 is considered significant.

*3.3. Effect of* CYP2C9 *and* VKORC1 *Polymorphisms on Warfarin Required Dose during Initiation Phase of Therapy*

Carriers of *CYP2C9* and *VKORC1* polymorphisms had a significantly increased required dose compared with wild-type subjects or carriers of only one polymorphism of *CYP2C9* or *VKORC1* (Table 4).

\* Kurskal Wallis test with *p*-value < 0.05 is considered significant, Mean Standard deviation in square brackets.

#### *3.4. Effect of* CYP2C9 *and* VKORC1 *Polymorphisms on Warfarin Responsiveness during Initiation of Therapy*

There were no significant differences in patient responder groups regarding the *VKORC1* and *CYP2C9* SNPs except for the *VKORC1* rs17708472 (*p* = 0.042) and the *CYP2C9* rs4086116 (*p* = 0.005) SNPs (Table 5). However, significant associations were found between genetic haplotypes of CCGG *VKORC1* and TAT *CYP2C9* and warfarin sensitivity, with *p* = 0.02 and *p* = 0.018, respectively (Table 6).

#### *3.5. Effect of* CYP2C9 *and* VKORC1 *Polymorphisms on INR Treatment Outcome*

There were no significant differences observed between the *CYP2C9* SNP genotypes and INR values measured at start of treatment for 212 cardiovascular patients treated with warfarin. In contrast, significant differences were observed between INR values measured at the initiation phase of therapy and certain *VKORC1* SNPs, namely rs10871454 (*p* = 0.006), rs8050894 (*p* = 0.007), and rs9934438 (*p* = 0.009), as shown in Table 4.

**Table 5.** Association of *VKORC1* and *CYP2C9* SNPs with response to warfarin during the initiation phase of therapy of 212 cardiovascular patients.

\* Chi-Square Test with *p*-value < 0.05 is considered significant.

**Table 6.** Frequencies of the haplotypes of *VKORC1* and *CYP2C9* genes among the 212 warfarin responsiveness patients.

\* Genetic haplotype frequency of 212 warfarin intake patients, \*\* *p*-value < 0.05 is considered significant.

#### *3.6. Correlation Between Warfarin Dose and Clinical Data*

Finally, there was no significant correlation between warfarin dose and body mass index, age, gender, co-morbidities, or the treatment indication (*p* = 0.505).

#### **4. Discussion**

Earlier studies on warfarin pharmacogenetics provide evidence that common *VKORC1* and *CYP2C9* polymorphisms with clinical and environmental factors are responsible for over half of the variability in warfarin required dose [23,24]. Genotyping patients who are carriers of *VKORC1* and *CYP2C9* variant alleles has been proven to reduce the risk of over-anticoagulation compared to the traditional initial dose approach [25,26].

In this study, our goal was to identify genetic factors associated with sensitivity and responsiveness to warfarin treatment during the initiation phase of treatment in Jordanian-Arab patients with cardiovascular disease. The results of the current pharmacogenetic study strongly suggest that there is a significant association of the *VKORC1* rs8050894, rs10871454, rs9934438, and rs17708472 SNPs and the *CYP2C9* rs4086116 and rs1057910 SNPs and their haplotypes with the required warfarin dosage and warfarin sensitivity. This study also reported that there is a genetic association between *VKORC1* rs17708472 SNP and CCGG genetic haplotype block and *CYP2C9* rs4086116 SNP and TAT genetic haplotype block with warfarin responsiveness during the initiation phase of therapy.

The allelic frequencies of *CYP2C9* and *VKORC1* SNPs in our population were similar to those found in other ethnic groups, as in the case of *CYP2C9\*2*, with 10% in our population, 6% in American and European populations, and 1% in Africans [27]. However, allelic frequencies of *VKORC1* SNPs were found to differ drastically from other populations. For example, rs10871454 was 52% in our population compared to 41% in Americans, 39% in Europeans, and 6% in African populations [27]. With regard to the association of *CYP2C9* polymorphisms with warfarin sensitivity, our results are consistent with the study by Takahashi et al. (2001) which shows that *CYP2C9* \*2 and \*3 polymorphisms reduce warfarin clearance [28] as *CYP2C9* is the major metabolizing enzyme of warfarin, therefore, reduction of activity results in lower required doses needed to achieve the therapeutic INR. We found a strong association of *CYP2C9\*3* (rs1057910 A>C) and *CYP2C9* (rs4086116 C>T) genotypes with warfarin sensitivity during the initiation stage of treatment with *p* < 0.001 and *p* = 0.012, respectively (Table 2). This study reported that individuals with one variant allele were associated with an increased risk of warfarin sensitivity. For example, 41.5% of the patients who carried the rs1057910 A>C variant allele were sensitive to warfarin, compared to 8.8% of the wild-type patients were sensitive. Moreover, carrying a *CYP2C9* TCC genetic haplotype block was significantly associated with warfarin sensitivity with *p* < 0.0001 (Table 3).

In disagreement with other studies, the *CYP2C9\*2* variant did not show a significant association with warfarin sensitivity (*p* = 0.744). This can possibly be explained by the genotypic frequency having an impact on the association; in our population only two patients were homozygous for the variant allele (TT), 32 were heterozygous patients (CT), and 105 patients were homozygous for the wild-type allele (CC) (Table 2). Moreover, it has been proposed that patients who carry one *CYP2C9\* 2* allele results in a dose reduction compared with the wild-type dose [14,29,30]. Although the allelic frequency of *CYP2C9\*2* (rs1799853) in our samples was in agreement with the American and European populations [27], this study did not find significant differences between this SNP and variability in required doses in the initiation phase of therapy (*p* = 0.366) as shown in Table 4 and Figure 2.

Clinical pharmacogenetic studies suggested that patients who carry the *CYP2C9\*3* (rs1057910) C allele leads to a dose reduction of 28–41% [12,29,30] in the Caucasian American population and from 12–38% of the Asian population compared to the wild-type [31–34]. In alignment with the aforementioned studies, we revealed that there is a 34.3% reduction in warfarin dose. Patients carrying one variant allele (C) required 41.75 mg/week, in comparison with the wild-type allele which required 27.43 mg/week (*p* < 0.001), as shown in Table 4 and Figure 2. Therefore, in our study the *CYP2C9\*3* allele has a greater effect on variation in warfarin dose during the initiation phase of therapy compared with *CYP2C9\*2*. In the case of individuals carrying rs4086116 C>T variant allele, this resulted in 23.9% and 37.9% reduction on warfarin dose compared to wild-type with *p* = 0.016, as shown in Table 4 and Figure 2.

**Figure 2.** The distribution of warfarin dose by CYP2C9 genotypes during the initiation phase of therapy for 212 Jordanian cardiovascular patients: X axis represents different CYP2C9 genotypes, Y axis represents the proportion of patients across each genotype, Blue column represents a sensitive group who required the lowest warfarin dose (<21 mg/week), Purple column represents the intermediate group who required moderate warfarin dose ((21–49) mg/week), Yellow column represents a resistant group who required the highest warfarin dose (>49 mg/week). (**A**) Distribution of warfarin dose by rs1799853 variant. (**B**) Distribution of warfarin dose by rs1057910 variant. (**C**) Distribution of warfarin dose by rs4086116 variant.

For the four studied *VKORC1* SNPs, we observed a strong association of *VKORC1* SNPs with warfarin sensitivity (*p* < 0.001). For example, patients with the T allele for rs10871454 C>T showed a high risk of warfarin sensitivity with 25.5% (CT) and 48.3% (TT) reduction of the required dose, respectively. In this case, the drug target enzyme could be expressed in smaller amounts and, therefore, low doses of the drug can obtain a therapeutic INR in an initial phase of therapy (Table 2). Furthermore, *VKORC1* genetic haplotype analysis showed a significant association between three *VKORC1* genetic haplotype blocks and sensitivity to warfarin with *p* < 0.0001 (Table 3). Moreover, Schelleman et al. (2007) and Wadelius et al. (2005) reported that a correlation exists between *VKORC1* SNP 1173 C> T (rs9934438), and the variation in warfarin dose patients carrying the variant allele of this SNP is related to a reduction in the required dose compared to the wild-type [5,35]. In alignment with these results, we found that patients who carry one variant allele required an average dose of 39.47mg/week and two variant allele carriers needed an average of 26.82 mg/week, while patients carrying the wild type CC needed an average dose of 53.02 mg/week with *p* < 0.001 as shown in Figure 3. Limdi et al. (2007) and Shrif et al. (2011) showed that *VKORC1* rs8050894 (1542G>C) were associated with lower warfarin doses in European Americans and Sudanese patients, respectively [36,37]. Accordingly, our results found a significant association between lower required warfarin dose and this SNP at an initial phase of therapy with *p* < 0.001 as shown in Figure 3.

**Figure 3.** The distribution of warfarin dose by *VKORC1* genotypes during the initiation phase of therapy for 212 Jordanian cardiovascular patients: X axis represents different *VKORC1* genotypes, Y axis represents the proportion of patients across each genotype, Blue column represents a sensitive group who required the lowest warfarin dose (<21 mg/week), Purple column represents the intermediate group who required moderate warfarin dose ((21–49) mg/week), Yellow column represents a resistant group who required the highest warfarin dose (>49 mg/week). (**A**) Distribution of warfarin dose by rs10871454 variant. (**B**) Distribution of warfarin dose by rs9934438 variant. (**C**) Distribution of warfarin dose by rs8050894 variant. (**D**) Distribution of warfarin dose by rs17708472 variant.

For the last *VKORC1* SNP *VKORC1*\*4 C<T (rs17708472), our study is similar to the study by Haug et al. (2008), which reported that this SNP was associated with higher dose requirements [38]. Our results showed that this variant was associated with significant differences in initial warfarin required dose; patients who were homozygous for the variant allele (TT) genotype required an average dose of 57.8 mg/week, heterozygous (CT) patients required an average dose of 41.63 mg/week, while wild-type (CC) patients required an average dose of 33.65 with *p* < 0.001 (Table 4).

With regard to the correlation of *CYP2C9* and *VKORC1* SNPs and warfarin responsiveness, we compared SNP genotypes with the warfarin responder groups (poor, good, and extensive responders). Significant differences were found between *VKORC1* rs17708472 (C>T) genotypes and the three different responder groups (Table 5); 33.3% of the patients carrying the variant allele (TT) were within the extensive responder group (meaning this variant allele was associated with increased risk of over-anticoagulation), compared to 11.3% of the wild-type (CC) patients who were extensive responders (*p* = 0.042). Accordingly, Kringen et al. (2011) have also shown that patients who carry this SNP are associated with an increased risk of the existence of therapeutic INR (over-anticoagulation) [39]. Otherwise, *VKORC1* rs10871454, rs8050894, and rs9934438 alleles show no significant differences between SNP genotypes within the three responder groups in our population, with *p* = 0.171, 0.235, and 0.226, respectively (Table 5). In contrast, *VKORC1* CCGG genetic haplotype block showed a significant association with warfarin responsiveness with *p* = 0.02 (Table 6).

Moreover, significant differences were observed between *VKORC1* rs10871454, rs8050894, and rs9934438 SNPs and INR value during the initiation phase of therapy (*p* = 0.006, 0.007, and 0.009, respectively), while rs17708472 SNP showed no significant differences (*p* = 0.493). Therefore, in our population the *VKORC1* SNP genotypes are associated with the generation of a high or low INR during the initiation phase of therapy (Table 4).

Taube et al. (2000) reported that an individual carrying an allelic variant of *CYP2C9* was not associated with an increased incidence of severe over-coagulation during long-term treatment [40]. Correspondingly, our results show no significant differences between the *CYP2C9\*2* and \*3 genotypes within the three responder groups with *p* = 0.076 and 0.910, respectively (Table 5). Conversely, our study reported significant differences between the proportion of *CYP2C9* rs4086116 (C>T) genotype and the three different responder groups; 50% of TT carriers were within the extensive responder group compared with 10.6% of wild type CC carriers (*p* = 0.005), which means that TT carriers are associated with increased risk of over-anticoagulation (Table 5). Moreover, significant association was observed between *CYP2C9* TAT genetic haplotype block and warfarin responsiveness with *p* = 0.018 (Table 6).

In addition, we did not observe significant differences in the three studied SNPs (*CYP2C9\*2, \*3,* and *CYP2C9* (C> T) rs4086116) and the INR value during the initiation phase of therapy. Therefore, in our population *CYP2C9* is not associated with the generation of a high or low INR as shown in Table 4.

Confirmation of our results and ongoing research including additional factors will be accomplished in a larger patient cohort, including genetic factors such as *OATP* transporters (mediates the uptake of warfarin into hepatocytes), *CYP3A4*, *CYP1A1*, and *CYP1A2* enzymes (metabolizing of R-warfarin), or *GGCX* encoded *gamma-glutamyl carboxylase* (the reduced vitamin K–form to activate coagulation factors) [5,41]. Application for individualized warfarin treatment will be both beneficial and efficient for cardiovascular patients in the future. Finally, the majority of the population included in this study is elderly, with only 15% of the subjects under 40 years of age. Therefore, additional study is needed in children and young adults.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/9/12/578/s1. Table S1: SNPs ID, their position and genotyping data based on whole cohort (*N* = 220). Table S2: List of SNPs, their minor allele frequencies, and HWE *p* values for genotypic distribution at each locus based on 220 patients. Table S3: The frequency of the allele and genotype for all polymorphisms in cardiovascular patients treated with warfarin.

**Author Contributions:** L.N.A.-E. designed the study. L.N.A.-E., A.Y.A. and R.H.K. were responsible for clinical data and blood samples collection. L.N.A.-E. and A.Y.A. analyzed and interpreted the data. L.N.A.-E. and A.Y.A. prepared the manuscript. All authors helped in reviewing the manuscript.

**Funding:** This work was supported by the Deanship of Research at Jordan University of Science and Technology under grant number 203/2014.

**Acknowledgments:** We are grateful to Pharmacist Nadia Al-Omary, and engineer of Nutrition Miss Hanna from anticoagulation clinic for cooperation and providing excellent samples for this study and for coordinating the Anticoagulant Service, and Miss Diana, laboratory technician, for the blood collection and INR measurements.

**Conflicts of Interest:** The authors declare no conflict of interest, financial or otherwise.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **Azathioprine Biotransformation in Young Patients with Inflammatory Bowel Disease: Contribution of Glutathione-S Transferase M1 and A1 Variants**

**Marianna Lucafò 1,2, Gabriele Stocco 3,\*, Stefano Martelossi 2, Diego Favretto 2, Ra**ff**aella Franca 4, Noelia Malusà 5, Angela Lora 2,4, Matteo Bramuzzo 2, Samuele Naviglio 2, Erika Cecchin 1, Giuseppe To**ff**oli 1, Alessandro Ventura 2,4 and Giuliana Decorti 2,4**

Received: 19 February 2019; Accepted: 1 April 2019; Published: 4 April 2019

**Abstract:** The contribution of candidate genetic variants involved in azathioprine biotransformation on azathioprine efficacy and pharmacokinetics in 111 young patients with inflammatory bowel disease was evaluated. Azathioprine doses, metabolites thioguanine-nucleotides (TGN) and methylmercaptopurine-nucleotides (MMPN) and clinical effects were assessed after at least 3 months of therapy. Clinical efficacy was defined as disease activity score below 10. Candidate genetic variants (*TPMT* rs1142345, rs1800460, rs1800462, *GSTA1* rs3957357, *GSTM1*, and *GSTT1* deletion) were determined by polymerase chain reaction (PCR) assays and pyrosequencing. Statistical analysis was performed using linear mixed effects models for the association between the candidate variants and the pharmacological variables (azathioprine doses and metabolites). Azathioprine metabolites were measured in 257 samples (median 2 per patient, inter-quartile range IQR 1-3). Clinical efficacy at the first evaluation available resulted better in ulcerative colitis than in Crohn's disease patients (88.0% versus 52.5% responders, *p* = 0.0003, linear mixed effect model, LME). TGN concentration and the ratio TGN/dose at the first evaluation were significantly higher in responder. *TPMT* rs1142345 variant (4.8% of patients) was associated with increased TGN (LME *p* = 0.0042), TGN/dose ratio (LME *p* < 0.0001), decreased azathioprine dose (LME *p* = 0.0087), and MMPN (LME *p* = 0.0011). *GSTM1* deletion (58.1% of patients) was associated with a 18.5% decrease in TGN/dose ratio (LME *p* = 0.041) and 30% decrease in clinical efficacy (LME *p* = 0.0031). *GSTA1* variant (12.8% of patients) showed a trend (*p* = 0.046, LME) for an association with decreased clinical efficacy; however, no significant effect on azathioprine pharmacokinetics could be detected. In conclusion, GSTs variants are associated with azathioprine efficacy and pharmacokinetics.

**Keywords:** azathioprine; inflammatory bowel disease; glutathione-S transferase; pharmacogenetics; pharmacokinetics

#### **1. Introduction**

Inflammatory bowel disease (IBD) is a chronic, relapsing and remitting disease of the gastrointestinal tract that comprises two main entities, Crohn's disease (CD) and ulcerative colitis (UC). The disease has a peak onset in subjects 15 to 30 years old, and its incidence is rising in the pediatric population [1]. Despite the recent introduction in therapy of biologicals, thiopurines continue to be widely used in this disease; indeed, these are cheap drugs, and maintain at least 20% of patients in a state of stable long term steroid free clinical remission [2]. Among thiopurines, azathioprine is mainly used as an immunosuppressant in IBD, and, although it has a well described risk benefit profile, adverse drug reactions are relatively common, occurring in 15–18% of patients, and can be severe enough to require the withdrawal of therapy [3,4] In addition, a significant proportion of patients does not respond to therapy with this agent [5–7]. The reasons of this high heterogeneity in clinical response is not clear yet [8]; however, variability in azathioprine metabolism can be important; indeed, azathioprine is a prodrug that requires metabolic conversion to its active form. The first step in this metabolic conversion is mediated by conjugation with glutathione, resulting in the formation of mercaptopurine. This reaction is in part nonenzymatic but it is even controlled, as demonstrated by recent publications, by the enzyme glutathione S-transferase (GST in particular by the isoforms A and M [9]. The latter is also inactive and is converted by the enzymes of the purine salvage pathway to the active thioguanine nucleotides, which are responsible of the cytotoxic and apoptotic effect of these drugs [10]. Mercaptopurine is metabolized by the enzyme xanthine oxidase in the liver, and by thiopurine methyltransferase (TPMT) and inosine triphosphatase (ITPA), mainly in extra hepatic tissues [10]. Polymorphisms in genes involved in azathioprine metabolism can hence influence the efficacy and toxicity of this drug [6].

In this study, we aimed to evaluate the contribution of candidate genetic variants involved in azathioprine biotransformation on azathioprine efficacy and pharmacokinetics in young patients with IBD.

#### **2. Materials and Methods**

#### *2.1. Patient Characteristics*

111 patients with IBD were enrolled by the Gastroenterology Unit of the Pediatric Hospital "Burlo Garofolo" in Trieste, Italy between March 2004 and February 2015. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Institutional Ethics Committee (Projects RC 23/2005 and 12/2013). All subjects and parents gave their informed consent for inclusion before they participated in the study. The inclusion criteria were age less than 30 years, previous diagnosis of IBD and treatment with azathioprine for at least 3 months. The patients enrolled are all the patients taking azathioprine at "Burlo Garofolo" in Trieste in the time-frame of the study. Blood samples for azathioprine metabolites measurement and for genotyping were taken at the appropriate clinic visit. Timing of metabolite level measurement was determined by the clinical setting of azathioprine administration at the hospital: generally, azathioprine metabolites levels were measured after 3, 6, and 12 months of treatment and then every year. Patients were treated with a dose-escalating strategy to reduce the risk of adverse events starting, however, from a relatively high dose (median of 2 mg/kg). At subsequent follow-up visits (every 3 months), the dose was increased or reduced so as to obtain the optimal clinical response; the criteria used to increase or reduce the dose of azathioprine were the level of disease activity and laboratory parameters used to monitor azathioprine toxicity (in particular leukocyte, erythrocyte and platelet counts, hemoglobin concentration, mean corpuscular volume, liver enzymes alanine aminotransferase, aspartate aminotransferase and γ-glutamyltransferase, and amylase levels). According to current guidelines, genotyping information was shared with the clinicians only for patients presenting *TPMT* variant alleles, in order to allow increased monitoring of adverse events.

Clinical response was assessed using Pediatric Crohn's Disease Activity Index and Pediatric Ulcerative Colitis Activity Index, respectively [11], for CD and UC patients at the time of blood sample collection for the first metabolites' measurement, which occurred at least 3 months since the beginning of therapy. The disease was considered inactive if the disease activity index was <10 at the time of sample collection.

#### *2.2. Measurement of Azathioprine Metabolites*

Metabolites (TGN andMMPN) were measured in patients' erythrocytes using the high performance liquid chromatography assay by Dervieux and Boulieu [12]. The ratio between TGN and the dose of azathioprine was calculated considering, for each individual measurement of the metabolites, the dose the patients took the day the blood sample was recorded.

#### *2.3. Genotypes*

Genomic DNA was extracted from peripheral blood samples using a commercial kit (Sigma, Milan, Italy), to characterize genetic polymorphisms in the candidate genes *TPMT* (rs1142345, rs1800460 and rs1800462), *GSTA1* (rs3957357), *GSTM1* (deletion), and *GSTT1* (deletion). The considered genotypes and method of analysis are described in Table 1. Genotypes for *TPMT* rs1800462 was determined by polymerase chain reaction (PCR) with allele specific oligonucleotides (ASO). Primers used were, for the wild-type allele, as forward P2W 5 -GTATGATTTTATGCAGGTTTG-3 and as reverse P2C 5 -TAAATAGGAACCATCGGACAC-3 ; primers. For the variant allele, a second tube was used with P2M 5 -GTATGATTTTATGCAGGTTTC-3 as forward primer and the above-mentioned P2C 5 -TAAATAGGAACCATCGGACAC-3 as reverse. PCR protocol for these primers were: initial denaturation 5 min at 94 ◦C, followed by 37 cycles with 30 s at 94 ◦C, 30 s at 57 ◦C, and 2 min at 72 ◦C, with a final extension for 10 min at 72 ◦C. PCR product was visualized on a 2% agarose gel. In case of a patient carrying the wild type allele, the product (254 bp) was present with the P2W and P2C primers; in case of patients carrying the variant allele, with the P2M and P2C primers for *TPMT* rs1800460 and rs1142345, PCR- restriction fragment length polymorphism (RFLP) was used. For rs1800460, primers used were: forward 5 -AGGCAGCTAGGGAAAAAGAAAGGTG-3 and reverse 5 -CAAGCCTTATAGCCTTACACCCAGG-3 . PCR protocol for these primers was: initial denaturation 5 min at 94 ◦C, followed by 37 cycles with 30 s at 94 ◦C, 30 s at 55 ◦C, and 2 min at 72 ◦C, with a final extension for 10 min at 72 ◦C. The DNA amplification produces an amplicon of 694 bp, which is subsequently digested enzymatically with the enzyme Mwol (concentration of 1 U/10 μl) incubated for 90 min at 60 ◦C. The enzyme recognizes the wild-type site and cuts the DNA strand into two fragments of 443 bp and 251 bp, while it does not cut the variant fragment. A 2% agarose gel was prepared for visualization. For rs1142345, primers used were forward 5 -AATCCCTGCTGTCATTCTTCATAGTATTT-3 and reverse 5 -CACATCATAATCTCCTCTCC-3 . PCR protocol was the same as *TPMT* rs1800460. PCR produces a 401 bp amplicon, which is subsequently digested enzymatically with the Accl enzyme (concentration 1 U/10 μl) and incubated for 90 min at 37 ◦C. The enzyme recognizes the variant site and cuts the DNA strand into two 252 pb and 149 pb fragments while the wild-type strand is not cut. A 2% agarose gel was prepared for visualization. *GSTM1* and *GSTT1* genotypes were determined by MULTIPLEX-PCR-ASO as previously described [13], in which three pairs of primers were used simultaneously: a specific pair for the T isoform, one for the M isoform and one for the β-globin gene, which acts as an internal positive control in order to verify the amplification. The three pairs of primers lead to three fragments of different sizes: 480 bp (*GSTT*), 286 bp (β-globin), and 219 bp (*GSTM*). The primers used have the following sequence: *GSTM* Forward: 5 -GAACTCCCTGAAAAGCTAAAGC-3 ; *GSTM* Reverse: 5 -GTTGGGCTCAAATATACGGTGC-3 ; β*-Globin*Forward: 5 -GAAGAGCCAAGGACAGGT-3 ; β*-Globin* Reverse: 5 -CAACTTCATCCACGTTCACC-3 ; *GSTT* Forward: 5 -TTCCTTACTGGTCCTCACATCTC-3 ; *GSTT* Reverse: 5 -TCACCGGATCATGGCCAGCA-3 . PCR protocol for these primers were: initial denaturation 5 min at 94 ◦C, followed by 37 cycles with 30 s at 94 ◦C, 30 s at 57 ◦C, and 2 min at 72 ◦C, with a final extension for 7 min at 72 ◦C. All PCR reactions described were carried out using RedTaq

polymerase (Sigma, Milan, Italy), with the addition of dNTPs 0,25 nM and with a primer concentration of 1 mM.

**Table 1.** Assay used for genotyping of the considered variants.

ASO: allele specific oligonucleotides; *GST*: glutathione-S-transferase, *TPMT*: thiopurine-S-methyl transferase; RFLP: restriction fragment length polymorphism. PCR: polymerase chain reaction.

For *GSTA1*, pyrosequencing was employed (Table 1), since this genotyping method was already validated in the laboratory. The primers used for the pyrosequencing were: forward 5 -ATCCAGTAGGTGGCCCCTTG-3 , reverse 5 -ACCGTCCTGGCTCGACAA-3 (biotinylated). Sequencing primer was: 5 -GCTTTTCCCTAACTTGAC-3 . PCR protocol for these primers were: initial denaturation 10 min at 95 ◦C, followed by 40 cycles with 30 s at 95 ◦C, 30 s at 66 ◦C, and 30 s at 72 ◦C and with a final extension for 10 min at 72 ◦C. PCR produces a 148 bp amplicon. For pyrosequencing, we used PSQ96MA (Qiagen, Hilden, Germany). PCR amplifications were performed in an Eppendorf Mastercycler gradient, with TaqGold DNA Polymerase (AB Applied Biosystems, Foster City, CA, USA).

#### *2.4. Statistical Analysis*

Statistical analysis was performed using the software R (version 2.15). The association between pharmacological phenotypes of interest (i.e., clinical efficacy of treatment, dose of azathioprine, TGN metabolites concentrations, MMPN metabolites concentrations, ratio TGN/dose) and the considered demographic variables, IBD type, cotreatment or genotypes in a univariate analysis, was evaluated using linear mixed effects model built using the phenotype as the dependent variable, each covariate as the fixed effect and the patients as the random effect in the model. For clinical efficacy, the first available measurement was used, while for other pharmacological variables, all available measurements were used.

Multivariate analysis was carried out to test the independence of the effects of the genotypes significant in the univariate analysis on the phenotypes considered (i.e., TGN or MMPN concentrations, dose of azathioprine, ratio TGN/dose); for this multivariate analysis generalized linear models of the Gaussian family were used considering individually each phenotype from the univariate analysis as the dependent variable and the covariates significant in the univariate analysis as the independent variables. Normality of the phenotype was tested by the Shapiro test and log10 transformation was applied if needed, in order to adjust the normality of the distribution.

#### **3. Results**

#### *3.1. Patients Enrolled and Samples Collected*

111 patients were enrolled from March 2004 to February 2015; median age was 15.05 years (IQR 12.28–16-82), 52 (46.8%) were females. Clinical and demographic characteristics of the enrolled patients are reported in Table 2.

**Table 2.** Demographic and clinical characteristics of the patients enrolled in the study.

For continuous variables, median, 1st–3rd quartiles values are reported. To report age and length of treatment median and interquartile range are provided; for patients with more than one measurement of azathioprine metabolites, median age and length of treatment were used.

Azathioprine metabolites were measured in 257 samples (median 2 per patient, IQR 1-3). Among these, 89 were obtained during treatment with azathioprine alone and 161 during treatment with azathioprine and other medications and in particular: 93 with an aminosalicylate, 18 with an aminosalicylate and a glucocorticoid, 15 with infliximab, 10 with an antibiotic, 4 with an aminosalicylate and an antibiotic, 4 with an antibiotic and a glucocorticoid, 3 with an aminosalicylate, an antibiotic and a glucocorticoid, 2 with an infliximab and a glucocorticoid, and 1 with an aminosalicylate and infliximab; for 7 patients, information about concomitant treatment could not be retrieved.

*3.2. Measurement of Azathioprine Metabolites: Association with Demographic and Clinical Covariates*

Results of measurements together with azathioprine dose are shown in Table 3.

**Table 3.** Summary of azathioprine's dose and metabolites' concentrations.

MMPN indicated methylated nucleotides, TGN indicates thioguanine nucleotides.

Concentration of TGN metabolites were associated with IBD type, with UC patients showing slightly increased concentrations (Supplementary Figure S1, LME *p* = 0.047), but not with gender or treatment length. Azathioprine dose was strongly associated with age, with younger patients taking higher doses (Supplementary Figure S2, LME *p* = 0.0001), but not with gender, IBD type or treatment length. Concentration of MMPN metabolites or the ratio between TGN concentration and azathioprine dose were not associated with IBD type, gender, or treatment length. Interestingly, the ratio between TGN concentration and azathioprine dose was strongly associated with azathioprine dose when the analysis was limited to pediatric patients (i.e., with age less than 18, Supplementary Figure S3, LME *p* = 0.0043). Clinical efficacy, defined as disease activity score below 10 at the time of first sample collection for measurement of azathioprine metabolites, was assessed in all patients. Azathioprine was more effective in UC than in CD patients (88.0% versus 52.5% responders, LME *p* = 0.0003), while gender, age, and duration of azathioprine treatment were not associated with azathioprine efficacy. A higher concentration of TGN metabolites at the first evaluation was observed in patients in remission (Figure 1, LME *p* = 0.0099), similarly a positive correlation was observed with TGN/dose ratio (Figure 1, LME *p* = 0.0023).

**Figure 1.** Response to azathioprine (AZA) and thioguanine nucleotides (TGN) concentration, as pmol/8 × 10ˆ8 erythrocytes (left panel) or ratio between TGN concentration/daily azathioprine dose as pmol/8 × 10ˆ8 erythrocytes/mg/kg/day (right panel). *p*-values are from linear mixed effect model (LME).

On the contrary, azathioprine dose and the concentration of MMPN metabolites were not associated with a clinical response (data not shown).

#### *3.3. Genotyping*

Results of genotyping are reported in Table 4.

**Table 4.** Genotype distribution in the 111 patients enrolled in the study.

GST indicates glutathione-S-transferase, *TPMT* indicates thiopurine-S-methyl transferase.

All polymorphisms evaluated respected Hardy-Weinberg equilibrium and their distribution was in accordance with literature data for subjects of Caucasian ethnicity. For the association between genetic variants and azathioprine pharmacokinetics, *TPMT* rs1142345 variant (4.8% of patients) was associated with increased TGN (LME *p* = 0.0042), TGN/dose ratio (LME *p* < 0.0001), decreased azathioprine dose (LME *p* = 0.0087) and MMPN (LME *p* = 0.0011; Figure 2), as well established [2]. Interestingly, all patients with variant *TPMT* were in remission at the first evaluation of thiopurine metabolites, in comparison to 65% of patients with wild-type *TPMT* (LME *p* = 0.041, Figure 2).

**Figure 2.** *TPMT* rs1142345 and azathioprine (AZA) dose, thioguanine (TGN) and methylmercaptopurine (MMPN) metabolites and efficacy. Concentration of azathioprine metabolites is expressed as pmol/8 × 10ˆ8 erythrocytes (U). *p*-values are from linear mixed effect model (LME).

*GSTM1* deletion (58.1% of patients) was associated with a 18.5% decrease in TGN/dose ratio (LME *p* = 0.041, Figure 3) and 30% decrease in clinical efficacy (LME *p* = 0.0031; Figure 3). Additionally, MMPN was reduced in patients with deletion of *GSTM1* (LME *p* = 0.039; Figure 3).

**Figure 3.** *GSTM1* deletion and azathioprine (AZA) dose, thioguanine (TGN), and methylmercaptopurine (MMPN) metabolites and efficacy. Concentration of azathioprine metabolites is expressed as pmol/8 × 10ˆ8 erythrocytes (U). *p*-values are from linear mixed effect model (LME).

*GSTA1* variant (12.8% of patients) showed a trend for an association with decreased clinical efficacy (LME *p* = 0.046, Figure 4); however, no significant effect on azathioprine pharmacokinetics could be detected (Figure 4).

**Figure 4.** *GSTA1* rs3957357 variant and azathioprine (AZA) dose, thioguanine (TGN) and methylmercaptopurine (MMPN) metabolites and efficacy. Concentration of azathioprine metabolites is expressed as pmol/8 × 10ˆ8 erythrocytes (U). *p*-values are from linear mixed effect model (LME).

*GSTT1* deletion was not associated with azathioprine pharmacokinetics and efficacy (data not shown). Multivariate analysis supported the results of the univariate analysis (Table 5).

**Table 5.** Multivariate analysis considering for each pharmacological dependent variable covariate significant in the univariate analysis.

*GST*: glutathione-S-transferase, IBD: inflammatory bowel disease, UC: ulcerative colitis, CD: Crohn's disease, MMPN: methylated nucleotides, TGN: thioguanine nucleotides, *TPMT*: thiopurine-S-methyl transferase. The effect size represents the increase (positive value) or decrease (negative value) in the value of the dependent variable for each independent variable listed. *p*-values are from linear mixed-effect models.

#### **4. Discussion**

Despite the introduction of new and effective biologics in the therapy of IBD, the thiopurine drugs azathioprine and mercaptopurine continue to be frequently used for maintaining remission in these diseases. The problem with these drugs is that they are ineffective in a significant percentage of patients, and also induce side effects that can be severe [2]. The reasons for this variability are not clear; however, a number of studies have suggested that variations in enzymes involved in their metabolism can be involved.

For azathioprine, this agent is the prodrug of mercaptopurine, and has to be converted to produce its pharmacological activity. This conversion can occur spontaneously, but is also catalyzed by the enzymes GST, in particular the A and M isoforms [16]. In rat liver homogenates, Kaplowitz et al. have demonstrated that, while at high pH (pH = 8.0) the nonenzymatic and enzymatic reactions occur at similar levels, at pH closer to physiological values, the enzymatic reaction prevails [17]. In addition, in homogenates of human livers obtained from transplant donors, treatment with furosemide, an inhibitor of soluble GSTs [18], inhibited the conversion of azathioprine to mercaptopurine [19]. Additional evidence of a role of GSTs in azathioprine metabolism has been obtained in animal models; indeed, pretreatment of rats with the GST inhibitor probenecid increased the proportion of azathioprine in rat liver and reduced GSH consumption. Similarly, less hepatic GSH depletion was observed after azathioprine treatment in Gunn rats, a model of hyperbilirubinemic rat [17]. Of interest, bilirubin is also a GST inhibitor, with some studies indicating a stronger inhibitory effect of bile acids on GSTM1 in comparison to other isoforms [20].

We previously showed that the frequency of *GSTM1* deletion was significantly lower in patients that developed an adverse event in comparison to patients that tolerated azathioprine treatment, in agreement with a model in which patients with *GSTM1* deletion are less sensitive to the effects of azathioprine, putatively because of the contribution of this enzyme on the conversion of azathioprine to mercaptopurine [13]. Moreover, in a recent previous study [21], we evaluated the effects of GST polymorphisms on azathioprine biotransformation in a cohort of young patients with IBD, tolerant to azathioprine therapy and taking the drug for more than 3 months. Patients with the deletion of *GSTM1* tolerated a dose of azathioprine significantly higher in comparison to patients with normal *GSTM1*. Moreover, the amount of active TGNs generated in patients with the deletion of *GSTM1* was significantly decreased in comparison to patients with a normal genotype. Multivariate analysis confirmed that this effect was independent from that of other genes with a significant effect, such as *TPMT*, the main gene known to influence mercaptopurine metabolism [22].

The present study is the first report of an association between azathioprine efficacy and *GSTM1* and *GSTA1* variants in young patients with IBD. Moreover, we confirmed the reduced TGN/azathioprine dose ratio in patients with *GSTM1* deletion we previously reported, which may be associated with the described lower efficacy of azathioprine in patients with this genotype. This could support the need for genotype adjusted tailored therapy, possibly testing the efficacy of strategies leading to higher TGN concentration in patients with *GSTM1* deletion, such as increased azathioprine dose or co-treatment with an aminosalycilate [23], even if prospective studies are needed to further support these strategies.

Therefore, all these studies support a role of *GSTM1* on azathioprine efficacy, mediated by an increased conversion of azathioprine to mercaptopurine. The reaction catalyzed by GSTM1 likely occurs after oral administration mainly in the intestine and the liver, modulating the amount of mercaptopurine and TGNs that are released in the main circulation [24].

Azathioprine dose is strongly associated with patients' age in the present study, an observation consistent with our previous results in children with IBD, showing that these patients require higher doses of azathioprine to achieve similar therapeutic efficacy and TGN concentration [25]. *TPMT* activity indeed is significantly higher in children than in adults [26]; interestingly, in pediatric patients (age less than 18 years), we could observe a lower ratio of TGN/azathioprine dose, as in our previous report. However, when the analysis was extended to young adults (age less than 30 years), the correlation between age and the TGN/dose ratio was lost. This may be related to environmental factors, including epigenetic determinants, even if more studies are needed to elucidate the mechanisms underlying these observations [27].

Considering studies by other groups, our results are in agreement with a recently published paper describing a lower efficacy of azathioprine in patients with *GSTM1* deletion, even if the results were not fully significant [28]. Moreover, in our study, we observed reduced concentration of MMPN nucleotides during azathioprine treatment in patients with *GSTM1* deletion: this result is consistent with a recent study by Broekman et al. [29]; this study also supports a lack of effect of *GSTA1* variants on azathioprine TGN and MMPN concentrations. The clinical implications of these observations need to be further explored. Additionally, age may affect the association of GSTs variants with thiopurine effects; indeed, studies in adult patients could not identify a consistent effect of *GSTM1* variants on thiopurines induced adverse effects [30], and therefore, other studies in the pediatric population are needed.

A recent study investigated the association among *GST* polymorphism, enzyme activity and azathioprine-related adverse drug reactions in Chinese Han patients with IBD, finding that the patients who became neutropenic had a significantly higher GSTs activity when compared with patients who did not develop toxicity [31]. The authors found, in the univariate analysis, that *GSTM1* wildtype genotype had a relationship with leukopenia and flue like symptoms, while *GSTP1* variant was strongly associated with leukopenia. Following adjustment for other potential risk factors, it was shown that *GSTP1* variants only were associated with increased risks of leukopenia. In our current study, we did not consider the effect of *GSTP1* polymorphisms on azathioprine effects and metabolism, since in our previous studies no significant association with adverse effect [13] or biotransformation [21] could be detected. The lack of association may be due to the tissue distribution of *GSTP1* and *GSTT1*, which are not highly expressed in the liver, but even to the lack of specific activity of these enzymes toward the catalysis of the reaction of azathioprine with glutathione [9,16]. Since the GSTP isoform does not catalyze the biotransformation of azathioprine to mercaptopurine [16], other mechanisms could be involved in the association observed by Liu and collaborators, such as induction of oxidative stress or modulation of apoptosis [24]. The study by Liu et al. reporting an effect of *GSTP1* variants on azathioprine induced adverse events, with a milder effect of *GSTM1*, seems to underline that the

effects of GSTs on azathioprine pharmacogenetics may be influenced by ethnicity. Indeed, it is already known that variants frequent in Asian patients but uncommon in other ethnic groups, are associated with increased sensitivity to thiopurines, such as *MRP4* and *NUDT15* [32,33].

One limitation of our study is its retrospective design and the consequent difficulty of properly assessing phenomena such as drug interaction, which should rely on data collection from patients' charts. Drug interactions between azathioprine and other agents employed in IBD have been described; in particular, a significant decrease in TGN levels after discontinuation of aminosalicylates has been previously reported [23]. Moreover, for *GSTA1*, a marginal effect on thiopurine efficacy was observed but this was not supported by an effect on thiopurine pharmacokinetics: this may be related to the limited number of patients homozygous for the *GSTA1* variant enrolled. Another limitation is the fact that the current assay for thiopurine metabolites quantifies two main species (thioguanine nucleotides and methyl-mercaptopurine nucleotides), without distinguishing between monophosphate, diphosphate, and triphosphate nucleotides. Innovative mass spectrometry based assays are now available to quantify thiopurine metabolites [34], allowing quantification of phosphorylation of thionucleotides [35] and they could be applied to evaluate differences in thiopurine biotransformation in patients with various GST genotypes. Evaluation of the combined effects of genotypes in this study is limited. Indeed, multivariate analysis indicates independency in the effects of the candidate genotypes considered on the pharmacological variables in the present cohort. A larger cohort is needed to detect significant effects by combined genotypes. For demographic covariates, in particular gender, no significant effect was identified in the univariate analysis; to further evaluate interactive effect of gender and the considered genotypes, a larger cohort is needed.

In conclusion, GSTs variants were associated with azathioprine efficacy and pharmacokinetics; more studies, both clinical and molecular are still needed to apply this evidence to improve outcomes of therapy with azathioprine in young patients with IBD.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/10/4/277/s1, Figure S1: inflammatory bowel disease (IBD type) and azathioprine thioguanine nucleotide (TGN) metabolites, Figure S2: Patient's age and azathioprine daily dose, Figure S3: Patient's age and ratio between concentration of azathioprine thioguanine-nucleotide metabolites and azathioprine daily dose.

**Author Contributions:** Conceptualization, M.L., G.S., E.C., G.T., A.V., and G.D.; Data curation, M.L., G.S., A.L., M.B., S.N., E.C., and G.D.; Formal analysis, M.L., G.S., D.F., N.M., E.C., and G.D.; Funding acquisition, S.M., M.B.; Methodology, M.L., S.M., D.F., R.F., N.M., A.L., M.B., S.N., and E.C.; Supervision, S.M., N.M., E.C., G.T., A.V., and G.D.; Writing—original draft, M.L., G.S. and G.D.; Writing—review & editing, M.L., G.S., S.M., E.C., A.V., and G.D.

**Funding:** This study was supported by the Italian Ministry of Health (Progetti Ricerca Corrente 23/2005, 12/2013 and 21/17).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **MicroRNAs Mediated Regulation of Expression of Nucleoside Analog Pathway Genes in Acute Myeloid Leukemia**

#### **Neha S. Bhise 1,2, Abdelrahman H. Elsayed 1, Xueyuan Cao 3, Stanley Pounds 4 and Jatinder K. Lamba 1,\***

Received: 26 March 2019; Accepted: 20 April 2019; Published: 24 April 2019

**Abstract:** Nucleoside analog, cytarabine (ara-C) is the mainstay of acute myeloid leukemia (AML) chemotherapy. Cytarabine and other nucleoside analogs require activation to the triphosphate form (ara-CTP). Intracellular ara-CTP levels demonstrate significant inter-patient variation and have been related to therapeutic response in AML patients. Inter-patient variation in expression levels of drug transporters or enzymes involved in the activation or inactivation of cytarabine and other analogs is a prime mechanism contributing to development of drug resistance. Since microRNAs (miRNAs) are known to regulate gene-expression, the aim of this study was to identify miRNAs involved in regulation of messenger RNA expression levels of cytarabine pathway genes. We evaluated miRNA and gene-expression levels of cytarabine metabolic pathway genes in 8 AML cell lines and The Cancer Genome Atlas (TCGA) data base. Using correlation analysis and functional validation experiments, our data demonstrates that miR-34a-5p and miR-24-3p regulate DCK, an enzyme involved in activation of cytarabine and DCTD, an enzyme involved in metabolic inactivation of cytarabine expression, respectively. Further our results from gel shift assays confirmed binding of these mRNA-miRNA pairs. Our results show miRNA mediated regulation of gene expression levels of nucleoside metabolic pathway genes can impact interindividual variation in expression levels which in turn may influence treatment outcomes.

**Keywords:** nucleoside analogs; microRNAs; gene expression; drug resistance; AML

#### **1. Introduction**

Nucleoside analogs (NA) are a class of chemotherapeutic agents that structurally resemble the endogenous purine or pyrimidine nucleosides. These therapeutic agents mimic the endogenous nucleosides with respect to their uptake and metabolism and are incorporated into the newly synthesized DNA leading to inhibition of DNA synthesis and chain termination. Some of the nucleoside analogs also inhibit or block the enzymes that are required for the synthesis of purine or pyrimidine nucleotides and RNA synthesis, leading to the activation of the caspase cascade and cell death. The nucleoside analogs are extensively used for the treatment of both hematological malignancies and solid tumors.

The pyrimidine nucleoside analog, cytarabine, is one of the most widely used chemotherapeutic drugs for the treatment of acute myeloid leukemia (AML).

One of the major obstacles in the treatment of AML is development of resistance to nucleoside analogs. There is a growing need to understand the mechanisms that lead to development of resistance to these nucleoside analogs in order to help identify strategies that would effectively treat patients with relapsing or refractory diseases. One of the primary mechanisms of resistance to nucleoside analogs is insufficient intracellular concentration of the active triphosphate metabolite. This insufficient triphosphate levels could be due to inefficient cellular uptake of the drug, reduced levels of the activating enzyme, increased levels of inactivating enzymes and/or due to increased levels of endogenous deoxynucleotide (dNTP) pools [1–4]. Resistance could also develop due to inability to achieve sufficient alterations in the DNA strands or the dNTP pools, either due to altered interaction with DNA polymerases or by a lack of inhibition of ribonucleotide reductases, or due to inadequate p53 exonuclease activity. Since the expression and activity of drug transporters and metabolizing enzymes in the activation pathway of nucleoside analogs plays an important role in development of resistance to the NAs, it is essential to understand the factors influencing the expression and activity of these proteins.

MicroRNAs (miRNAs) are a group of novel gene regulators, which have been recently recognized to play an important role in cancers due to their tumor suppressive and oncogenic functions [5]. MiRNAs are known to regulate the expression of the target genes by binding to the specific sequence mainly on the 3' untranslated region on the genes. Role of miRNAs in regulating the expression of various drug-metabolizing enzymes like cytochrome P450 3A4 (CYP3A4) etc., drug transporters like BCRP and various drug targets [6–8] have been established. However, there have not been any studies that have comprehensively evaluated the effect of miRNAs on the important genes involved in the transport, activation and inactivation of nucleoside analogs. Hence, the aim of this study was to assess the effect of miRNAs on the expression of nucleoside analog pharmacokinetic (PK) and pharmacodynamic (PD) pathway genes and in turn assessing their potential impact on resistance to nucleoside analogs.

#### **2. Materials and Methods**

#### *2.1. Cell culture and Reagents*

The AML cell lines HL-60, MV-4-11, Kasumi-1, THP-1, AML-193 and KG-1 cell lines were obtained from American Type Culture Collection (ATCC) (Manassas, VA, USA), while MOLM-16 and ME-1 cell lines were obtained from DSMZ (Braunschweig, Germany). Kasumi-1, ME-1 and MOLM-16 cell lines were cultured in Roswell Park Memorial Institute (RMPI)-1640 medium supplemented with 20% fetal bovine serum (FBS), THP-1 cell line was cultured in RPMI-1640 medium supplemented with 10% FBS, HL-60 and KG-1 cell lines were cultured in IMDM medium supplemented with 20% FBS, while AML-193 cell lines was cultured in IMDM medium supplemented with 5% FBS, 0.005 mg/mL insulin, 0.005 mg/mL transferrin and 5 ng/mL granulocyte/macrophage colony stimulating factor (GM-CSF). All the cell lines were maintained in a 37 ◦C humidified incubator with 5% CO2. The cells were passaged every two to three days in order to maintain them in logarithmic growth phase.

#### *2.2. RNA Isolation*

Total RNA was isolated from the AML cell pellets using RNeasy Plus Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's protocol and stored in −80 ◦C until further analysis. The RNA quality and concentration were measured using NanoDrop 2000 UV-Vis spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The ratio of absorbance at 260 nm and 280 nm was used to assess RNA sample purity and A260/A280 ratio of 1.8–2.1 was considered to be indicative of highly purified RNA. RNA was normalized to 0.2 μg/μL with nuclease-free water before being used for performing reverse transcription reactions, as recommended by the manufacturer. The total RNA was

reverse transcribed to complementary DNA (cDNA) using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA) according to manufacturer's protocol.

#### *2.3. Gene Expression Analysis*

The expression of nucleoside analog genes was determined using the TaqMan® Low Density Array (TLDA) cards (Applied Biosystems). Each TLDA card was custom designed with pre-loaded gene expression assays for measuring the messenger RNA (mRNA) expression of selected nucleoside analog metabolic pathway genes- (*n* = 14) deoxycytidine kinase (*DCK*), cytidine deaminase (*CDA*), solute carrier family 29, member 1 (*SLC29A1*), solute carrier family 28, member 1 (SLC28A1), solute carrier family 28, member 3 (SLC28A3), deoxycytidylate deaminase (*DCTD*), 5'-nucleotidase, cytosolic II (*NT5C2*), 5'-nucleotidase cytosolic III (*NT5C3*), cytidine 5'-triphosphate synthase (CTPS), cytidine monophosphate kinase (*CMPK*), nucleoside diphosphate kinase 1 (*NME1*), ribonucleotide reductase M1 (*RRM1*), ribonucleotide reductase M2 (*RRM2*), ribonucleotide reductase M2B (*RRM2B*). Each TLDA card consists of eight separate loading ports that fill into 48 separate wells, for a total of 384 wells per card. Thus, each card could analyze the expression of 24 different genes for eight different samples in duplicates. Each cDNA sample was added to equal volume of 2X TaqMan Universal PCR Master Mix (Thermo Scientific) and 100 μL of the sample-specific PCR mix was added to the fill reservoir on the TLDA card. The card was centrifuged twice for one minute at 1200 rpm to distribute the sample-specific PCR reaction mix to the reaction wells. The card was sealed using the TaqMan Array Micro Fluidic Card Sealer (Thermo Scientific) and placed on microfluidic card thermal cycling block of Applied Biosystems 7900HT Fast Real-time PCR System (Applied Biosystems). Thermal cycling conditions were as follows: 2 min at 50 ◦C, 10 min at 94.5 ◦C, 30 s at 97 ◦C, 1 min at 59.7 ◦C for 40 cycles. The target mRNA expression levels were normalized to GAPDH and the expression values of nucleoside analogs pathway genes were calculated using ΔΔCT method [9].

#### *2.4. MicroRNA Expression Analysis*

For determination of miRNA expression, total RNA was isolated using mirVana™ miRNA Isolation kit (Life Technologies, Carlsbad, CA, USA) as per the manufacturer's protocol. The RNA quality and concentration were measured using NanoDrop 2000 UV-Vis spectrophotometer (Thermo Scientific). A total of 100 ng of purified total RNA was used for nCounter miRNA sample preparation reactions according to manufacturer's instructions and was assayed for determination of 800 human miRNA expression using the nCounter Human v2 miRNA Expression Assay kit (Nanostring Technologies, Seattle, WA, USA). Preparation of small RNA samples involved multiplexed ligation of specific tags (miRtags) to the target miRNAs that provide unique identification for each miRNA species. After ligation, the detection was done by hybridization to microRNA: tag specific nCounter capture and barcoded reporter probes. Data collection was carried out using the nCounter Digital Analyzer (Nanostring Technologies) at The University of Minnesota Genomics Center, following manufacturer's instructions to count individual fluorescent barcodes and quantify the target miRNA molecules present in each sample. MiRNA expression data normalization was performed using the nSolver™ Analysis Software (Nanostring Technologies) according to the manufacturer's instructions. In particular, initially the data was normalized using the expression of the top 100 code sets. Further, to account for the background correction, mean of negative controls plus two-standard deviation (SD) method was used. In order to avoid using the miRNAs with a very low expression, we further filtered out the miRNAs with expression counts < 30 (2 times the mean ± 2 SD of negative control value), in order to account for the background noise. Total 412 miRNAs with expression counts > 30 were included for further analysis.

#### *2.5. Acute Myeloid Leukemia Patient Sample Data from The Cancer Genome Atlas*

The miRNA expression and mRNA expression of the nucleoside analog pathway genes in AML patients was extracted from The Cancer Genome Atlas (TCGA) Data Portal (cancergenome.nih.gov) [10]. Out of the 200 AML patients in TCGA database, 197 patients had gene expression profiling data available and 187 patients had miRNA expression data available. 186 patients had both gene expression and miRNA expression data available.

#### *2.6. Electrophoretic Mobility Shift Assays*

The functional validation for binding efficiencies between miRNAs and mRNAs was performed using the electrophoretic mobility shift assays (EMSAs). The binding free energy between the respective mRNA and miRNA pair was predicted using the RNAhybrid software. The miRNA oligonucleotides were labeled with cy5™ dye on their 5' ends. The 2' O-methyl-modified mRNA oligonucleotides were labeled with IRDye®800 (LI-COR Biosciences, Lincoln, NE, USA) dye on their 5' ends. The labeled oligonucleotides were synthesized by Integrated DNA Technologies (Coralville, IA, USA). RNA EMSA experiment was performed using the LightShift Chemiluminescent RNA EMSA Kit (Thermo Scientific) according to the manufacturer's protocol. The mRNA oligonucleotide was heated for 10 min at 80 ◦C and then placed on ice in order to relax the secondary structures. In each 20 μL binding reaction, 200 nM miRNA oligonucleotide and/or mRNA oligonucleotide were mixed with RNA EMSA binding buffer and incubated at 25 ◦C for 25 min. The reaction mixtures were separated on a 12% polyacrylamide gel by electrophoresis at 4 ◦C. The binding reactions were transferred onto nylon membrane and the resulting mobility shifts were imaged using and Odyssey CLx Infrared System (LI-COR Biosciences).

#### *2.7. Bioinformatic Analysis*

Prediction of miRNA binding sites was performed using multiple prediction programs, which use different criteria for prediction of binding sites: TargetScan (www.targetscan.org), miRanda (www.microRNA.org), PICTAR (pictar.mdc-berlin.de), miRWalk (www.umm.uni-heidelberg.de/apps/ zmf/mirwalk). Binding free energy calculations were performed using RNAhybrid software [11]. The 3'UTR (3' untranslated region) sequence of mRNA was obtained from the UCSC Genome browser (https://genome.ucsc.edu/) and miRNA sequence was obtained from miRBase software (http://www.mirbase.org/).

#### *2.8. Statistical Analysis*

The nonparametric Spearman correlation was used to measure the correlation of mRNA expression with miRNA expression. Statistical significance was determined when *p*-value was < 0.01.

#### **3. Results**

#### *3.1. E*ff*ect of Micro RNA on the Expression of Nucleoside Analog Pathway Genes in Acute Myeloid Leukemia Cell Lines*

We determined expression of 800 miRNAs and 13 genes involved in PK/PD pathway of nucleoside analogs (Figure 1) in cytogenetically different AML cell lines (*n* = 8). In order to identify the miRNAs associated with the expression of nucleoside analog pathway genes, we correlated the miRNA expression and mRNA expression using the spearman correlation analysis. Table 1 lists the negative correlations of nucleoside analog pathway genes and respective miRNAs at *p* <0.01. We further used CyTargetLinker software [12] to establish the network of miRNAs correlated with the respective PK/PD pathway genes of the nucleoside analogs. Figure 2 shows the miRNA-mRNA pairs identified by CyTargetLinker.

The expression of *DCK* (the rate-limiting enzyme in the nucleoside analog pathway) correlated with the expression of miR-34a-5p expression (spearman *r* = −0.88; *p*-value < 0.01) and miR-96-5p expression (spearman *r* = −0.91; *p*-value < 0.01). The expression of deactivating enzyme *DCTD* was found to be correlated with miR-24-3p expression (spearman *r* = −0.93; p-value < 0.01). Interestingly, in our previous study, we identified that expression of miRNA miR-24-3p was correlated with cytarabine-induced cell cytotoxicity (spearman *r* = −0.81, *p*-value < 0.05) [13] which is in agreement

with the current observation. The expression of *CMPK*, a kinase responsible for phosphorylation of the monophosphate form of nucleoside analog was negatively correlated with the expression of miR-1301, miR-1323, miR-320e, miR-381, miR-507, miR-584-5p, miR-605, miR-762, miR-769-3p, miR-891a (all *p*-values < 0.01). RRM2 expression was found to be negatively associated with the expression of miR-151a-3p (*p*-value < 0.01). Figure 3 shows correlation plots between *DCTD*-miR-24, *DCK*-miR34 and *NT5C3*-miR149 pairs.

**Figure 1.** Disposition pathway of nucleoside analogs, cytarabine and clofarabine. DCTD: Deoxycytidylate deaminase, DCK: Deoxycytidine kinase, CDA: Cytidine deaminase, NT5C2/3: 5'-Nucleotidase, cytoplasmic, CTPS: CTP synthase, RRM1/2: Ribonucleotide reductase, SLC29A: Solute carrier family 29, SLC28A: Solute carrier family 28, MPK: Monophosphate kinase, NDK: Nucleoside diphosphate kinase

**Figure 2.** MicroRNA-mRNA network constructed using CyTargetLinker [12]. MiRNAs associated with the nucleoside analog metabolic pathway genes are shown in yellow and the mRNAs are depicted in red. Blue lines indicated positive and red indicated negative associations.

**Figure 3.** Correlation of miRNAs with nucleoside analog pathway gene expression in acute myeloid leukemia (AML) cell lines and patient samples from The Cancer Genome Atlas (TCGA) database. (**A**) Correlation between DCTD and hsa-miR-24a-3p. (**B**) Correlation between DCK and hsa-miR-34a-5p. (**C**) Correlation between NT5C3 and hsa-miR-149.

#### *3.2. Bioinformatic Prediction of Binding of Micro RNAs and Messenger RNAs*

MiRNAs are known to guide the RNA-induced silencing complex (RISC) to the specific sequence (usually in located in the 3'UTR) on the target mRNA. Using various bioinformatic prediction programs (TargetScan, miRanda, PICTAR, miRWalk) we determined if the miRNAs that were correlated with gene expression had binding sites on 3'-UTR of target genes. MiRNAs miR-1323, miR-30d-5p, miR-381, and miR-605 were predicted to have binding sites on *CMPK* gene, while miRNA miR-24-3p was found to have binding site on *DCTD* by multiple prediction programs. Supplementary Figure S1 shows comparisons of different miRNA prediction programs for genes involved in nucleoside analog pathways.

#### *3.3. E*ff*ect of Micro RNAs on the Expression of Nucleoside Analog Pathway Genes in Acute Myeloid Leukemia Patients*

In order to validate the significant correlations between miRNAs and mRNAs identified in AML cell lines, we evaluated the correlation between miRNA expression and nucleoside analog pathway gene expression in AML patient samples from TCGA database (*n* = 186). We extracted the miRNA expression data and nucleoside analog pathway gene expression data for AML patients from TCGA database and performed spearman correlation to identify the significant mRNA-miRNA pairs. Consistent with results from AML cell lines miR-24-3p was inversely correlated with the expression of both the probes for *DCTD* in AML patient samples (*r* = −0.22; *p*-value < 0.01 and *r* = −0.21; *p*-value < 0.01, Figure 3A). *DCK*-mir-34a pair unfortunately did not show significant association within the TCGA data-set. miR-149 was significantly correlated with expression of *NT5C3* in AML patients (*r* = −0.25; *p*-value < 0.01, Figure 3C).

Validation of Binding Efficiencies between Messenger RNAs and Micro RNAs

MiRNAs are known to bind to the specific seed sequence on the 3'UTR of the mRNAs, thereby regulating the expression of their target genes. In order to validate the binding between the mRNA-miRNAs identified from the in vitro studies and in AML patient samples, we performed RNA EMSA assays. We validated the interaction between *DCTD* and miR-24, since we identified this mRNA-miRNA pair to be significantly inversely correlated in both AML cell line and AML patient samples. In addition, various prediction databases predicted miRNA miR-24 to have binding site on the 3'UTR of *DCTD* mRNA. In silico analysis predicted miR-24-3p and miR-34a-5p might form complexes with target sequences in the 3'UTR of *DCTD* and *DCK* respectively with minimum free energies of binding of −27.2 kcal/mol for *DCTD* and miR-24-3p (Figure 4A and Table 1) and -25.6 kcal/mol for *DCK* and miR-34a-5p (Figure 5A and Table 1). The RNA EMSA results for IRD-800®-labeled *DCTD* and Cy-5-labeled miR-24-3p show miR-24a-3p was able to bind to its target sequence on *DCTD* 3'UTR (Figure 4B, lane 3) as seen by the band shift. The thermodynamic stability of this complex correlated with binding observed in the RNA EMSA assays. In addition, we found that mRNA-miRNA complex formed by *DCTD* and miR-24a-3p could be eliminated by adding excess unlabeled hsa-miR-24a-3p probe (Figure 4B, lane 4), but not by adding excess unlabeled non-specific probe (Figure 4B, lane 6). Adding excess unlabeled mRNA probe resulted in binding of all the labeled miRNA giving a greater intensity signal (Figure 4B, lane 5).

Similar to the interaction observed between *DCTD* and miR-24a-3p, we observed that miR-34a-5p binds to *DCK* 3'UTR as seen by shift in the band (Figure 5B, lane 3) and this interaction was eliminated by addition of unlabeled probe (Figure 5B, lane 4). Since the minimum free energy of binding for *NT5C3*- miR-149 pair was not strong (Table 1), we did not pursue EMSA assays for this pair.

**Figure 4.** Validation of binding interaction between DCTD mRNA and has-miR-24-3p by RNA electrophoretic mobility shift assays (EMSAs). RNA EMSA with cy5-labeled has-miR-24-3p oligonucleotide and 2'-O-methyl modified and IRD-800 labeled DCTC mRNA oligonucleotide. Lanes 1 and 2 show the mobility of the labeled mRNA or miRNA oligonucleotide. Lane 3 shows the mobility of the labeled has-miR-24-3p oligonucleotide with *DCTD* mRNA oligonucleotide. Lanes 4 and 6 show the mobility of labeled *DCTD* mRNA oligonucleotide in presence of unlabeled excess specific competitor (has-miR-24-3p) ad excess unlabeled non-specific competitor (NC).

**Figure 5.** Validation of binding interactions between *DCK* and mRNA and has-miR-34a-5p by RNA EMSAs. RNA EMSA with cys5-labeled has-miR-34a-5p oligonucleotide and 2'-O-methyl modified and IRD-800 labeled DCK mRNA oligonucleotide. Lanes 1 and 2 show the mobility of the labeled mRNA or miRNA oligonucleotide. Lane 3 shows the mobility of the labeled has-miR-34a-5p oligonucleotide with DCK mRNA oligonucleotide. Lanes 4 and 6 show the mobility of labeled DCK mRNA oligonucleotide in presence of unlabeled excess specific competitor (has-miR-34a-5p) and excess unlabeled non-specific competitor (NC).

**Table 1.** MiRNAs with significant negative association (as reflected by *r* value) with nucleoside analog pathway genes in AML Cell lines.

We further evaluated the association of the in vitro chemosensitivity of ara-C (defined previously Bhise et al, 2015 [13] with the three top miRNAs in the AML cell lines. As shown in Supplementary Figure S2. consistent with our results of miR-24-*DCTD* pair, we observed cell lines that were sensitive to ara-C has significantly higher levels of miR-24 as compared to cell lines that are resistant to ara-C (*p* = 0.03). These results suggest that high miR-24 in ara-C sensitive cell lines might be resulting in lowering *DCTD* levels and given that *DCTD* is involved in inactivation of ara-C, its low levels will result in better response.

#### **4. Discussion**

Nucleoside analogs are synthetic analogs of endogenous nucleosides that largely used for the treatment of hematological malignancies and solid tumors. Cytarabine, a pyrimidine nucleoside analog is the backbone of AML chemotherapy, while clofarabine is a second-generation purine nucleoside analog that is currently being investigated for treatment of AML in various clinical trials. Both cytarabine and clofarabine require active transport into the cell by nucleoside transporters, followed by activation by various kinases to form active di- and tri-phosphate metabolites that are incorporated in growing DNA strand and/or inhibit various enzymes involved in synthesis of endogenous nucleotides (Figure 1). However, despite being the backbone of treatment regimen used in AML patients, there is variability in response to cytarabine and other nucleoside analogs. In our previous study, we have demonstrated that miRNA expression is predictive of response to cytarabine therapy in AML patients and is also significantly associated with in vitro chemosensitivity of cytarabine in AML cell lines [13]. Recent studies are expanding on therapeutic relevance of miRNAs in AML. Recent review article by Wallace et al, 2017, have summarized the miRNAs that are deregulated in AML and thus hold potential as AML biomarkers. Further, mimics for miR-22, miR29b and miR-181 and antagomiRs for miR-21/miR-196b and miR-126 are currently under investigation for their therapeutic potential in AML [14]. However, given that cytarabine a nucleoside analog is the mainstay of AML chemotherapy for decades, in the current study, we wanted to determine if miRNA mediated regulation of expression of the transporters, activating and inactivating genes involved in the metabolic pathway of nucleoside analogs, could contribute to development of resistance to nucleoside analogs in AML patients. We hypothesized that miRNAs could bind to the 3'UTR of the mRNAs of nucleoside analog metabolic pathway genes, thereby altering their expression, which in turn would result in lower intracellular levels of active nucleoside analog triphosphate, resulting in chemo-resistance.

In our current study, DCK expression was negatively correlated with the expression of hsa-miR-34a-5p and hsa-miR-96-5p in AML cell lines (*p*-value < 0.01). DCK is a rate-limiting enzyme that is involved in activation of cytarabine and clofarabine. Studies have reported that decreased or complete loss of DCK activity results in cellular resistance to cytarabine [15–17]. Also, *DCK* mRNA expression has been shown to be positively associated with AML patient outcome, AML patients with higher *DCK* mRNA expression demonstrated longer event-free survival than those with lower *DCK* mRNA expression [18]. Using EMSA, we were able to show that hsa-miR-34a-5p by binding to the 3'-UTR regulates expression of *DCK*. MiRNA hsa-miR-34a has been extensively studied in various cancers [19–23] and it has been shown to play an important role as a tumor suppressor by targeting various genes. However, the effect of miR-34a on *DCK* expression has not yet been studied. Identification of this additional regulatory mechanism for an important enzyme in the activation of nucleoside analogs could help in better prediction of chemosensitivity of these drugs. Unfortunately, we did not observe any significant correlation between *DCK* mRNA and miR-34a levels in TCGA dataset.

We also identified hsa-miR-24-3p to be negatively correlated with the expression of *DCTD* in both AML cell lines and in AML patient samples (*p*-value < 0.01). In addition, multiple bioinformatic prediction programs identified a binding site for hsa-miR-24-3p on the 3'UTR of DCTD. Our RNA EMSA results confirmed the binding interaction between DCTD and miR-24-3p. DCTD is an enzyme involved in deamination of the monophosphate form of the nucleoside analog, thus inactivating the drug. The levels of DCTD could thus affect the levels of the intracellular active triphosphate metabolites of nucleoside analogs. However, the role of *DCTD* in chemosensitivity of nucleoside analogs is poorly defined. Various studies have demonstrated a significant role of *DCTD* in metabolism of the monophosphate metabolite of the nucleoside analogs in human leukemia cells [24–26]. Sequencing of this gene identified a nonsynonymous SNP affecting the activity of *DCTD* in vitro [27]. Hence, the limited data on the regulation of *DCTD* gene warrants the need to evaluate additional mechanism regulating gene-expression. MiRNA miR-24 is has also been extensively studied in various cancers [28–34] and has been shown to enhance metastasis and invasion. Increased expression of miR-24 has been associated with increased risk of relapse and poor survival in acute lymphoblastic leukemia (ALL) [32]. In addition to miR-24 and hsa-miR-34a-5p, we also identified multiple other miRNAs that correlated with multiple genes in the PK/PD pathway (Table 1). We acknowledge the limitation of use of different platforms of miRNA and mRNA quantification between cell lines and the TCGA data-set which warrant the need for future prospective study to validate our results. Additionally, given that gel-shift assays do not consider interactions with argonaute which can impact thermodynamics of target recognition future in-depth mechanistic validation considering these factors are needed to establish miRNA-mRNA regulatory pairs of significant therapeutic implications in AML.

In summary, we identified several miRNAs, which were significantly associated with the expression of nucleoside analog pathway genes. Identification of these additional mechanisms of regulation would help provide a better understanding of the variability in the expression of these enzymes and transporters and in turn, help in better prediction of therapeutic response in AML patients. While additional functional studies are required to gain mechanistic understanding of these miRNA-mRNA interactions and its effect on the protein levels and activity, this study helps identify candidate miRNAs for further studies.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/10/4/319/s1, Figure S1: Comparison of microRNA prediction programs for predicting binding sites on nucleoside analog pathway genes. Figure S2: A) Eight AML cell lines were treated with varying concentration of ara-C for 48 hrs followed by measuring cell viability using MTT assays using as described previously (Bhise et al, 2015). Area under the survival was calculated using Graphpad Prism and cell lines were classified as sensitive or resistant to ara-C. B) miR24-3p levels were observed to be higher in cell lines sensitive to ara-C as compared to ara-C resistant cell lines.

**Author Contributions:** J.K.L. and N.S.B. were involved in designing the study; N.S.B., J.K.L., A.H.E., X.C. and S.P. performed statistical analysis. All authors contributed to manuscript writing. All authors read and approved the final manuscript.

**Funding:** This research was funded by National Institutes of Health-National Cancer Institute R01-CA132946 (Lamba and Pounds).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **Diagnostic and Therapeutic Strategies for Fluoropyrimidine Treatment of Patients Carrying Multiple** *DPYD* **Variants**

**Carin A. T. C. Lunenburg 1, Linda M. Henricks 2,3, André B. P. van Kuilenburg 4, Ron H. J. Mathijssen 5, Jan H. M. Schellens 2,3, Hans Gelderblom 1, Henk-Jan Guchelaar 6 and Jesse J. Swen 6,\***

Received: 9 November 2018; Accepted: 26 November 2018; Published: 28 November 2018 -

**Abstract:** *DPYD* genotyping prior to fluoropyrimidine treatment is increasingly implemented in clinical care. Without phasing information (i.e., allelic location of variants), current genotype-based dosing guidelines cannot be applied to patients carrying multiple *DPYD* variants. The primary aim of this study is to examine diagnostic and therapeutic strategies for fluoropyrimidine treatment of patients carrying multiple *DPYD* variants. A case series of patients carrying multiple *DPYD* variants is presented. Different genotyping techniques were used to determine phasing information. Phenotyping was performed by dihydropyrimidine dehydrogenase (DPD) enzyme activity measurements. Publicly available databases were queried to explore the frequency and phasing of variants of patients carrying multiple *DPYD* variants. Four out of seven patients carrying multiple *DPYD* variants received a full dose of fluoropyrimidines and experienced severe toxicity. Phasing information could be retrieved for four patients. In three patients, variants were located on two different alleles, i.e., in *trans*. Recommended dose reductions based on the phased genotype differed from the phenotype-derived dose reductions in three out of four cases. Data from publicly available databases show that the frequency of patients carrying multiple *DPYD* variants is low (< 0.2%), but higher than the frequency of the commonly tested *DPYD*\*13 variant (0.1%). Patients carrying multiple *DPYD* variants are at high risk of developing severe toxicity. Additional analyses are required to determine the correct dose of fluoropyrimidine treatment. In patients carrying multiple *DPYD* variants, we recommend that a DPD phenotyping assay be carried out to determine a safe starting dose.

**Keywords:** pharmacogenomics; pharmacogenetics; genotype; phenotype; alleles; precision medicine

#### **1. Introduction**

Fluoropyrimidines (including 5-fluorouracil (5-FU) and capecitabine) are the cornerstone of treatment for various types of cancer and are used by millions of patients worldwide each year [1–3]. However, up to one-third of treated patients experience severe toxicity (Common Terminology Criteria for Adverse Events (CTC-AE) grade ≥ 3), such as diarrhea, hand–foot syndrome, or mucositis upon treatment with fluoropyrimidines [4,5]. These adverse events can lead to mortality in approximately 1% of patients who experience severe toxicity [4,6]. Dihydropyrimidine dehydrogenase (DPD) is the key enzyme in the metabolism of 5-FU and its decreased activity is strongly associated with toxicity [7,8]. Variants in *DPYD*, the gene encoding DPD, can lead to decreased DPD enzyme activity [9–12]. Prospective *DPYD* genotyping of four main *DPYD* variants followed by dose reductions in patients carrying any of these four *DPYD* variants is safe, cost-effective, and feasible in clinical practice [13–15]. These *DPYD* variants are *DPYD*\*2A (rs3918290, c.1905+1G>A, IVS14+1G>A); *DPYD*\*13 (rs55886062, c.1679T>G, I560S); c.1236G>A/HapB3 (rs56038477, E412E); and c.2846A>T (rs67376798, D949V). For these four variants, convincing evidence has been provided warranting implementation in clinical practice [4,5,12,15–17].

An increasing number of hospitals apply prospective *DPYD* genotyping when treating patients with fluoropyrimidines [18]. Individual dosing guidelines for the abovementioned four *DPYD* variants are provided by the Dutch Pharmacogenetics Working Group (DPWG) and the Clinical Pharmacogenetics Implementation Consortium (CPIC) [19,20]. Dosing guidelines are based on the expected remaining DPD enzyme activity and can be applied to patients who are heterozygous carriers of a single *DPYD* variant. For homozygous *DPYD* variant allele carriers (two identical variants) and compound heterozygous *DPYD* variant allele carriers (two or more different variants), dosing guidelines are not yet available (or treatment with an alternative drug is advised), although safe treatment with low-dose fluoropyrimidines in these homozygous *DPYD* patients was demonstrated by a recent case series [21].

Patients who carry multiple variants (compound heterozygous) can carry the variants on a single allele (in *cis*) or on different alleles (in *trans*). In the first case, one functionally active allele remains, whereas in the latter case, both alleles are affected, which may result in a proportionally decreased enzyme activity (Figure 1). With currently used genotyping techniques, the allelic location of variants (phasing) cannot be determined. This uncertainty hampers adequate interpretation of the pharmacogenetic test result in compound heterozygous patients and makes it impossible to give an appropriate dose recommendation based on the genotype alone. Since it is likely that in the future, even more *DPYD* variants will be tested, the probability of finding compound heterozygous *DPYD* variant allele carriers will increase. The aims of this study are to examine diagnostic and therapeutic strategies for fluoropyrimidine treatment of patients carrying multiple *DPYD* variants and the frequency and phasing of variants of compound heterozygous *DPYD* patients in publicly available databases.

**Figure 1.** Illustration of zygosity and clinical interpretations. Black stars represent variants; boxes represent alleles. A wild-type patient carries no variants, resulting in normal-activity alleles (green). A heterozygous patient carries one variant, resulting in one reduced or inactive allele (red) and one active allele (green). A partly reduced enzyme activity is expected, since there is still one active allele left. For homozygous patients, both variants result in a reduced or inactive allele (red). Depending on the effect of the variants on the protein, a reduced or absent enzyme activity is expected. Compound heterozygous patients can carry variants on different alleles (in *trans*) or on one allele (in *cis*), resulting in differences in enzyme function, either like that of a heterozygous patient or a homozygous patient.

#### **2. Materials and Methods**

In this study, we present seven compounds heterozygous *DPYD* variant allele carriers as clinical cases. In addition, we have performed in silico research in publicly available databases.

#### *2.1. Patients*

Data and DNA from patient cases carrying multiple *DPYD* variants were collected. Patients were identified either after development of severe toxicity from fluoropyrimidine-containing therapy, by additional retrospective genotyping in a clinical trial (clinicaltrials.gov identifier NCT00838370, [13]), or prior to treatment in routine clinical care. The study was reviewed and approved by the institutional review board of the Leiden University Medical Centre (LUMC, G18.015). Patient data from the electronic medical records was handled following the codes of proper use and proper conduct in the self-regulatory codes of conduct [22]. Toxicity to fluoropyrimidine-containing therapy was graded by the treating physicians using the National Cancer Institute CTC-AE version 4.03 [23], and severe fluoropyrimidine-induced toxicity was defined as CTC-AE grade ≥ 3. In some cases, additional patient material to determine the phasing of the *DPYD* variants was collected. In these cases, additional patient consent was obtained.

#### *2.2. Dihydropyrimidine Dehydrogenase Enzyme Activity Measurements*

For all patients, DPD enzyme activity was determined. This could be either prior to treatment or retrospectively after the occurrence of severe toxicity. DPD enzyme activity measurement in peripheral blood mononuclear cells (PBMCs) [24,25] was used as a reference to assess DPD activity, and has been used previously to determine dosages in *DPYD* variant-carrying patients [21,26]. A validated method [27] was used, containing radiolabeled thymine as a substrate and consisting of high-performance liquid chromatography (HPLC) with online radioisotope detection using liquid scintillation counting. Normal values for healthy volunteers are 9.9 ± 2.8 nmol/(mg×h),

for DPD-deficient patients are 4.8 ± 1.7 nmol/(mg×h), and reference values range from 5.9 to 14 nmol/(mg×h) [28]. Dose reductions based on DPD enzyme activity were performed in a one-to-one ratio, as was previously described by Henricks et al. [21]. Thereafter, toxicity-guided dosing was used.

#### *2.3. Molecular Methods for Estimation of Phasing*

In regard to the size of the *DPYD* gene, the location of the variants, and the material available (DNA, RNA) from the patients, three molecular methods to determine the phasing of the variants could be used in this study. In four patients, we could execute one or more of these methods. These methods are explained and illustrated in the supplementary material (Figure S1). Details on these techniques have been published elsewhere [29–31].

#### *2.4. Frequencies of Compound Heterozygous DPYD Carriers*

To investigate the incidence of compound heterozygous *DPYD* variant allele carriers (of the four genotyped *DPYD* variants), large databases were queried [32,33]. The incidence was calculated using minor allele frequencies (MAFs) of each variant identified in the databases separately. Since the determined variants are not in the same haplotype, it was assumed that the inheritance of these individual *DPYD* variants is independent. All genotypes from the databases were calculated to be in Hardy–Weinberg equilibrium, except for *DPYD*\*2A and c.1236G>A for the Exome Aggregation Consortium (ExAC, http://exac.broadinstitute.org/) [32] and Genome Aggregation Database (gnomAD, http://gnomad.broadinstitute.org/) [33] due to a slight overrepresentation of homozygous cases. The calculated frequencies were compared to frequencies from databases in which phasing could be determined.

#### Exome Aggregation Consortium and Genome Aggregation Database

Both the ExAC [32] and gnomAD [33] databases collect exome sequencing data and aggregate the data for public use. The ExAC dataset (v0.3.1) contains sequenced data of 60,706 unrelated individuals. The gnomAD dataset (v2.0) contains sequenced data of 123,136 exomes and 15,496 genomes from unrelated individuals. In ExAC, 2791 *DPYD* variants, and in gnomAD, 2190 *DPYD* variants were found. MAFs of *DPYD* variants from these databases reflect those of the population due to the large group size in the databases. Since both ExAC and gnomAD do not contain individual matched or phased data, it is not possible to search for compound heterozygous patients in these databases.

#### *2.5. Phasing in Compound Heterozygous DPYD Carriers*

Three databases were used to identify compound heterozygous *DPYD* variant allele carriers and determine the phasing, i.e., allelic location, of variants.

#### 2.5.1. Genome of the Netherlands Datasets

The Genome of the Netherlands (GoNL, http://www.nlgenome.nl/) trio datasets contain information of related fathers, mothers, and children, and phasing information is therefore available. Datasets were previously processed and phased using trio-aware variant calling [34]. After the exclusion of children, phased variant call format (VCF) files for 496 subjects (fathers and mothers) were obtained from the GoNL repository. The toolset Bedtools (https://bedtools.readthedocs.io/en/latest/, v2) was used to extract all variants found in the *DPYD* locus (chr1:97,543,300–98,386,615). Next, for all individuals, the carrier status of *DPYD*\*2A, *DPYD*\*13, c.1236G>A, and c.2846A>T was examined. Individuals who carry at least one of the four actionable *DPYD* variants were identified, and using a custom Python [35] script, the phasing of variants was assessed for individuals with multiple variants.

#### 2.5.2. 1000 Genomes Database

The 1000 Genomes Project (http://www.internationalgenome.org/) is the largest publicly available catalogue of human variation and genotyped phased data. It originally ran from 2008 until 2015, and thereafter it was maintained and expanded by the International Genome Sample Resource (IGSR) [36]. On 27 October 2016, phased data of the *DPYD* gene (chr1: 97,543,300–98,386,615) was downloaded from the 1000 Genomes ftp server (phase 3; GRCh37; chr1: 97,543,300–98,386,615) using Tabix (v1.1) [37]. The statistical program R (v3.2.5) [38] was used to select the genotypes at four *DPYD* risk alleles in unrelated individuals of Caucasian descent.

#### 2.5.3. Exome Trios Leiden University Medical Centre Database

This diagnostic database of the clinical genetics department of LUMC contains 433 complete exome trios (father, mother, and child). The exome was enriched by the Agilent sureselect v5 kit and sequenced using various Illumina (San Diego, CA, USA) sequencers (Hiseq 2000, 2500, 4000, Nextseq). Carrier status of the abovementioned *DPYD* variants was established by querying the trio VCF files. We also investigated all samples with sufficient coverage of this region to obtain a reliable frequency estimate. In the case of trios, only parents were taken into account.

#### **3. Results**

#### *3.1. Patient Cases and Clinical Implications*

Details of the demographics and clinical characteristics of the seven cases are described in the supplementary material (patient cases). All patients received treatment with fluoropyrimidines and were identified as compound heterozygous *DPYD* variant allele carriers either prior to the start of treatment or retrospectively. Table 1 shows an overview of the cases. Table 2 shows all genotype and phenotype results. With additional genetic testing, phasing could be determined in four out of seven patients. In three patients, the variants were located in *trans*, and one patient carried the variants in *cis*. With the phasing information available, it is possible to calculate a dose recommendation using publicly available pharmacogenetic dosing guidelines [19,20]. For example, patient 1 carried *DPYD*\*2A and c.1236G>A in *trans*. The gene activity values range from inactive (0) to fully active (1). *DPYD*\*2A and c.1236G>A/HapB3 have values of 0 and 0.5, respectively. As this patient carries the variants in *trans*, each allele contains one variant and no fully functional allele remains. Therefore, the cumulated gene activity score (GAS) is 0.5. The GAS can be used to determine dose recommendations according to the genotype, as was previously described [12]. The GAS ranges from 0 to 2, and a score of 0.5 corresponds to a dose recommendation of 25%. The DPD enzyme activity of patient 1 was 0.9 nmol/(mg×h). This was divided by the mean of the reference value (9.9), which results in a theoretical DPD activity of 9%. For each patient for whom phasing details were known, the GAS was determined and compared to the theoretical DPD activity. Dose recommendations according to the GAS (genotype) and theoretical DPD activity (phenotype) were divergent in almost all cases, as shown in Table 2.

#### *3.2. Preventing Toxicity*

Three of the seven case patients were identified as carriers of one or more *DPYD* variants prior to the start of therapy. For one patient, the DPD enzyme activity was determined prior to the start of therapy. Based on their genotype or phenotype, these three patients received initially reduced fluoropyrimidine dosages of 50%. They experienced limited and reversible toxicity (CTC-AE grades 0–2). The dose of one patient was increased to 70% in the second treatment cycle, after which CTC-AE grade 3 toxicity occurred.

Four of the seven case patients received a full dose, since their genotype was unknown prior to the start of therapy. These patients all experienced severe toxicity (CTC-AE grades 3–5), and three of them were admitted to the hospital for 7–14 days. An overview of cases, including the toxicity, is shown in Table 3.

#### *3.3. Frequencies of Compound Heterozygous DPYD Carriers without Phasing Information*

The ExAC and gnomAD databases revealed an average MAF for *DPYD*\*2A, *DPYD*\*13, c.1236G>A, and c.2846A>T of 0.55%, 0.03%, 1.43%, and 0.27%, respectively. MAFs for ExAC and gnomAD separately are summarized in Table 4. The probability of identifying a compound heterozygous *DPYD* patient for two variants according to these databases was ≤ 0.008%, as was calculated using frequencies of combinations of *DPYD* variants. Results for each combination of *DPYD* variants are shown in Table 5. With several million fluoropyrimidine users each year, thousands of patients worldwide are compound heterozygous for a subset of these four *DPYD* variants.

#### *3.4. Frequencies of Compound Heterozygous DPYD Carriers with Phasing Information*

In the GoNL database, genetic data from 496 subjects (fathers and mothers only) was reviewed. One subject was found who carried two *DPYD* variants. This subject was a carrier of the *DPYD* c.1236G>A and *DPYD* c.2846A>T variants, both of which were located on a single allele (in *cis*). Based upon the data in GoNL, the probability of having compound heterozygosity of the four *DPYD* variants is <0.2%.

In the 1000 Genomes database, data of 2513 individuals were available. After the selection of unique, unrelated individuals, 407 individuals remained. One subject was found who carried two *DPYD* variants. This subject was a carrier of *DPYD* c.1236G>A and *DPYD* c.2846A>T, both of which were located on different alleles (in *trans*). Based upon the data in 1000 Genomes, the probability of having compound heterozygosity of the four *DPYD* variants is <0.3%.

In the LUMC clinical genetics database (exome trios LUMC), the analysis was restricted to the children, since this would allow phasing. None of the 433 children carried more than one *DPYD* variant, thus compound heterozygosity in this database is <0.2%.

Despite the low frequency, compound heterozygous patients were identified in all databases except the LUMC clinical genetics database. However, the low frequency did not allow to determine the probability of in *cis* or in *trans* phasing of variants in a patient.

**Table1.**Characteristicsofpatientcases.Shownperpatientareprimarytumor,treatment,capecitabinedose,executedassays(genotype,dihydropyrimidine

 The reference DPD activity ranges from 5.9–14 nmol/(mg×h) [28], and therefore the percentage of DPD activity can be calculated using the average of the reference (9.9 Thispercentagecouldbeusedasapercentageoftheregulardose.

*Genes* **2018** , *9*, 585

GoNL 489/7/0 0.7% 494/2/0 0.2% 475/21/0 2.1% 490/6/0

**Table 5.** Calculated frequencies for compound heterozygous *DPYD* patients. Using the average MAFs of the ExAC and gnomAD databases (for *DPYD*\*2A, *DPYD*\*13, c.1236G>A, and c.2846A>T, these are 0.55%, 0.03%, 1.43%, and 0.27%, respectively), possible combinations for two out of four currently genotyped *DPYD* variants are shown. *Abbreviations:* MAF: minor allele frequency; ExAC: Exome Aggregation Consortium; gnomAD: Genome Aggregation Database.

#### **4. Discussion**

Prospective genotyping of *DPYD* variants followed by individual dose adjustments is increasingly applied as the standard of care for patients starting fluoropyrimidine therapy. Standard dose reductions from CPIC and DPWG guidelines cannot be applied in patients who carry more than one *DPYD* variant, as the phasing of the variants is unknown. Despite the low population frequency of < 0.2%, the absolute number of identified compound heterozygous patients will increase as the number of genotyped patients increases and the panel of tested variants is expanded. To the best of our knowledge, this is the first study that describes a case series of compound heterozygous *DPYD* variant allele carriers and investigates diagnostic and therapeutic strategies for these patients.

Our study shows the clinical need for further information on the genotype, as four patients were identified as compound heterozygous carriers retrospectively and all of them experienced severe toxicity. These compound heterozygous *DPYD* variant allele carriers have an increased risk of developing severe fluoropyrimidine-induced toxicity if dosages are not adequately adjusted. Previously, compound heterozygous patients have been described with severe or even lethal side effects after fluoropyrimidine treatment [39,40]. Three patients in this study were prospectively identified as compound heterozygous carriers, received initial dose reductions, and developed only mild toxicities.

Out of the four patients for whom we were able to retrieve phasing information, three were in *trans* and one was in *cis* orientation. Data from publicly available databases also showed that both in *cis* and in *trans* orientations exist. However, the recently updated CPIC guidelines on *DPYD* assumes in *trans* phasing for compound heterozygous patients [20]. The DPWG guidelines do not mention phasing; however, the dosing recommendations of the DPWG use the GAS, a score based on the activity of individual alleles [19]. This implies the need for phasing information. The assumption of in *trans* phasing could result in the underdosing of patients with variants phased in *cis*, and thus exemplifies the need for the determination of the phasing of variants.

In this study, we looked at different diagnostic strategies to determine the phasing of *DPYD* variants in compound heterozygous patients. In four patients, the phasing of *DPYD* variants could be determined using one of three different molecular methods. These methods are in the early phases of development, not routinely available, quite expensive, and not always conclusive. For two of these techniques, patient RNA is used, which degrades quickly after the blood draw unless specifically designed blood tubes are used. Compound heterozygous patients are rare, yet here we describe seven patients heterozygous for multiple *DPYD* variants. A limitation of our study is that most patients were identified retrospectively and in different institutions. Because of this, not enough of or not the right material was available for analysis, thus not all genotyping techniques could be executed in each patient. For two samples, tests failed or produced inconclusive results. For this reason, a formal comparison of their suitability to identify phasing was not possible. However, of the three explored molecular methods, PacBio sequencing seems most promising. While phasing improved the prediction of DPD enzyme activity, patients with identical combinations of *DPYD* variants and identical phasing

showed considerable differences in DPD enzyme activity, which could potentially limit the added value of the determination of the phasing of *DPYD* variants. However, larger numbers of compound heterozygous *DPYD* variant allele carriers would be necessary to draw a firm conclusion.

The measurement of DPD enzyme activity in PBMCs was used as a reference to assess DPD activity. The method is well-established, commonly available, and shows limited intra- and interpatient variability [27]. However, recently, differences in intrapatient variability in DPD enzyme activity related to circadian rhythm were shown [41], which can result in the under- or overestimation of DPD enzyme activity. In this study, we present one patient with extremely low DPD enzyme activity, which could possibly be influenced by the presence of severe neutropenia, as DPD activity is normally measured in mononuclear cells. Therefore, DPD enzyme activity can differ depending on the clinical condition of the patient and should thus be measured prior to treatment.

A major question is whether genotyping or phenotyping is the best method to determine DPD activity to guide fluoropyrimidine dosing in patients carrying multiple *DPYD* variants. Despite the low population frequency, we present seven patients carrying multiple *DPYD* variants, of which three received initially reduced fluoropyrimidine dosages. However, based on these data, it is not possible to determine if a dose recommendation based on phased genetic information or DPD enzyme activity measured in PBMCs is safer. In three out of four cases, differences were observed between the theoretically calculated DPD activity using genotyping or phenotyping. These differences would result in different dosing recommendations. For example, there is a considerable interpatient variability in DPD enzyme activity in carriers of the *DPYD* variant c.1236G>A/HapB3 [12]. Due to this variability, genetic dose recommendations are categorized (e.g., 25 or 50%) on the average of the phenotypes. This categorization could explain the observed dosing differences derived from genotyping and phenotyping. Other variants of *DPYD* currently not routinely tested for or variants in other genes, e.g., *MIR27A* [42], might also be involved in reducing DPD activity or explaining fluoropyrimidine-induced toxicity. DPD enzyme activity measurements are well-established, and additional molecular methods to resolve phasing are still in early phases of development. Therefore, in our opinion, the current therapeutic strategy for compound heterozygous *DPYD* variant allele carriers should be to determine initial dose reductions based on a DPD phenotyping test, for example, by measuring enzyme activity in PBMCs. Dosing could be adjusted by the treating physician in subsequent cycles based on observed severe toxicity (or lack thereof).

#### **5. Conclusions**

In conclusion, patients carrying multiple *DPYD* variants are at high risk of developing severe toxicity. Additional analyses are required to determine the correct dose of fluoropyrimidine treatment. In patients carrying multiple *DPYD* variants, we recommend that a DPD phenotyping assay be carried out to determine a safe starting dose. The dose could be titrated in subsequent cycles based on observed toxicity.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/9/12/585/s1; Description of patient cases (1 to 7) and Figure S1 illustration of molecular methods.

**Author Contributions:** Conceptualization, C.A.T.C.L., H.G., H.-J.G. and J.J.S.; Data curation, C.A.T.C.L., L.M.H. and A.B.P.v.K.; Formal analysis, C.A.T.C.L., H.-J.G. and J.J.S.; Methodology, C.A.T.C.L., H.-J.G. and J.J.S.; Supervision, H.G., H.-J.G. and J.J.S.; Writing—original draft, C.A.T.C.L.; Writing—review & editing, L.M.H., A.B.P.v.K., R.H.J.M., J.H.M.S., H.G., H.-J.G. and J.J.S.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank S.Y. Anvar, G.W.E. Santen, E.B. van den Akker and H. Mei for assistance with the databases; R.H.A.M. Vossen, E.J. de Meijer and R.J.H.M van der Straaten for execution of the molecular methods for estimation of phasing; and F.M. de Man for the collection of patient material and clinical data.

**Conflicts of Interest:** The authors declare no conflict of interest. C. Lunenburg was previously supported by an unrestricted grant from Roche Pharmaceuticals. There was no involvement in the study design, data collection, analysis, or interpretation of the data.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **FARMAPRICE: A Pharmacogenetic Clinical Decision Support System for Precise and Cost-E**ff**ective Therapy**

**Rossana Roncato 1, Lisa Dal Cin 1, Silvia Mezzalira 1, Francesco Comello 1, Elena De Mattia 1, Alessia Bignucolo 1, Lorenzo Giollo 2, Simone D'Errico 2, Antonio Gulotta 2, Luca Emili 3, Vincenzo Carbone 3, Michela Guardascione 1, Luisa Foltran 4, Giuseppe To**ff**oli 1,\* and Erika Cecchin 1,\***

Received: 28 February 2019; Accepted: 1 April 2019; Published: 4 April 2019