Exome, transcriptome and miRNA analysis don’t reveal any molecular markers of TKI efficacy in primary CML patients

Background Approximately 5–20% of chronic myeloid leukemia (CML) patients demonstrate primary resistance or intolerance to imatinib. None of the existing predictive scores gives a good prognosis of TKI efficacy. Gene polymorphisms, expression and microRNAs are known to be involved in the pathogenesis of TKI resistance in CML. The aim of our study is to find new molecular markers of TKI therapy efficacy in CML patients. Methods Newly diagnosed patients with Ph+ CML in chronic phase were included in this study. Optimal and non-optimal responses to TKI were estimated according to ELN 2013 recommendation. We performed genotyping of selected polymorphisms in 62 blood samples of CML patients, expression profiling of 33 RNA samples extracted from blood and miRNA profiling of 800 miRNA in 12 blood samples of CML patients. Results The frequencies of genotypes at the studied loci did not differ between groups of patients with an optimal and non-optimal response to TKI therapy. Analysis of the expression of 34,681 genes revealed 26 differently expressed genes (p < 0.05) in groups of patients with different TKI responses, but differences were very small and were not confirmed by qPCR. Finally, we did not find difference in miRNA expression between the groups. Conclusions Using modern high-throughput methods such as whole-exome sequencing, transcriptome and miRNA analysis, we could not find reliable molecular markers for early prediction of TKI efficiency in Ph+ CML patients. Electronic supplementary material The online version of this article (10.1186/s12920-019-0481-z) contains supplementary material, which is available to authorized users.


Introduction
Chronic myeloid leukemia (CML) is a clonal myeloproliferative disorder characterized by the presence of [9;22] translocation and BCR/ABL fusion gene with high tyrosine kinase activity which activates MAPK pathway, cell proliferation, blocks apoptosis and leads to genome instability resulting in further development of the disease. Imatinib, BCR/ABL tyrosine kinase inhibitor (TKI) is the standard therapy in CML-Ph+ patients since its FDA approval in 2001. Patient's survival improved from 7.5 years after diagnosis before imatinib era to 17.5 years now-a-days [1].
Despite high efficacy of imatinib the problem of primary resistance persists. Based on the recent report about 21% of CML patients are switched to another TKI because of resistance or intolerance [2]. According to other authors approximately 20 to 30% of patients develop resistance to imatinib [3].
At the same time hematologists have limited instruments to determine which patients will have primary resistance to imatinib and may benefit from other treatment regimens or use of newly developed TKIs as the 1st line therapy. Sokal and Hasford scores were developed in pre-imatinib era and now poorly predict the outcomes of the TKI therapy while EUTOS score gives more reliable prediction [4]. However in some studies it was estimated that all three scoring systems didn't work correctly to predict complete cytogenetic response and survival with imatinib treatment, especially in non-European populations [5,6]. Currently, the main trend to predict the better outcome is to use the rule "deeper and earlier response" [7], but this approach allows only late prediction after starting the therapy.
Among other factors gene expression differences are of a particular interest since expression profiling was demonstrated to be a valuable prognostic tool in many types of cancer. It was previously shown that expression levels of several genes differ between groups with different responses to TKI therapy [14,15]. Frank et al. demonstrated that COCH, ANAPC5 and TPSAB1/B2 were overexpressed in imatinib non-responders while VNN1 and RPH3A were down-regulated. The magnitude of differences was dramatically lowfold change varied from 0.547 to 1.487 and was confirmed by qRT-PCR for only two genes. However 128-gene expression signature was successfully used to correctly classify the subset of test samples [14].
MircoRNAs (miRNAs) are small (18-25 nucleotides in length) non-coding RNAs which regulate gene expression by translational repression or mRNA cleavage [16]. Previously obtained data shows that miRNAs are involved in CML pathogenesis: some miRNAs are up-regulated and some are down-regulated in the peripheral blood of CML patients [17][18][19][20]. Moreover, there is data supporting the idea of different expression levels of miRNAs in CML patients with good and poor response to TKI therapy. San José-Enériz et al. [21] performed analysis of expression profiles of 250 miRNAs in bone marrow mononuclear cells from patients with Ph+ CML at diagnoses and showed that 19 miRNAs were differentially expressed in resistant and responder samples. Similar study was performed in peripheral blood samples by microarray analysis in two groups of patients with response and resistance to TKI. Authors identified 70 differently expressed miRNAs between these groups [22]. In both studies cluster unsupervised analysis of obtained expression levels of miRNAs was able to distinguish clearly both groups. It was also shown that miR-30 reduces BCR/ABL mRNA and protein levels by binding directly to the ABL 3′UTR and increases sensitivity of BCR/ABL-positive cells to imatinib. CML patients expressing low levels of miR-30 were less sensitive to imatinib [23]. High expression of miR-424 suppressed proliferation and induced apoptosis of K562 cells thereby increased sensitivity to imatinib treatment [24]. In another work it was shown that miR-26a, miR-29c, miR-130b and miR-146a were down-regulated in imatinib resistant patients in comparison to responders [25].
Despite the variety of the approaches and findings at DNA and RNA levels none of the potential markers were validated or implemented in clinical practice. The aim of our study is to find new molecular markers of TKI therapy efficacy in CML patients.

Patients
Newly diagnosed patients with Ph+ CML in chronic phase were included in this study. Patients were recruited in the National Research Center for Hematology (Moscow), Regional Children Hospital #1 (Ekaterinburg), Research Institute of Medical Cell Technologies (Ekaterinburg) and Rostov State Medical University (Rostov-on-Don). The study was performed according to the provisions of the Declaration of Helsinki. Informed consent was obtained from all patients included in this study. European Leuke-miaNet (ELN) Recommendations for Management of CML [26] were followed to form criteria of optimal response at 6 month: BCR-ABL IS < 1% and/ or 0% Ph+ metaphases out of at least 20 or more bone marrow cells; non-optimal response: BCR-ABL IS 1-10% and/or > 1% Ph+ metaphases out of at least 20 or more bone marrow cells at 6 months.

Genotyping of CML patients
Direct DNA Sanger sequencing of PCR products was performed as described previously [13]. Briefly, DNA from peripheral blood samples was obtained from CML patients before TKI therapy. Sequencing was performed using ABI Prism 3130XL genetic analyzer (Applied Biosystems, Foster City, CA) with either forward or reverse primer.

Transcriptome and Bioinformatic analysis
Blood samples were collected in Paxgene tubes (Becton, Dickinson and Company, USA) at diagnosis before initiating TKI therapy. RNA was extracted using PAXgene Blood RNA Kit (Qiagen, Germany) according to the manufacture's protocol. Expression profiling was performed using HumanHT Expression BeadChip (Illumina). Data analysis was performed in GenomeStudio_GX_Module (Illumina) and Partek Genomics Suite (Partek).

miRNA and Bioinformatic analysis
Total RNA was extracted from blood samples and lymphocytes cell cultures obtained before beginning the therapy at the time of diagnosis CML. RNA from Paxgene tubes (Becton, Dickinson and Company, USA) was extracted using PAXgene Blood RNA Kit (Qiagen, Germany) according to the manufacture's protocol; RNA from cell cultures was extracted by standard phenolchloroform method. miRNA profiling of 800 miRNA was done with nCounter miRNA Expression Assay (Nanostring Technologies) according to the manufacturer's protocol. Data analysis was performed in nSolver (Nanostring Technologies).

Statistical analysis
Statistical calculations were performed using STATIS-TICA 64 ver.13 (Dell Inc., USA), MS Office Excel, and SPSS Statistics 21.0. Chi-square tests were used to process the obtained genotyping data. A significance level of 0.05 was set and Bonferroni correction for multiple comparisons was applied. Mann-Whitney U test for the two independent groups was used for expression data analysis, patients' characteristics and comparison of numbers of obtained genetic variants. Haplotypes analysis was performed using R Haplo.stats [27]. Analysis of variance (ANOVA) and pathway enrichment analysis were performed in GenomeStudio_GX_Module (Illumina) and Partek Genomics Suite (Partek).

Results
Whole exome sequencing (WES) was used to find new potential markers of TKI efficacy [Lavrov AV, 2016].
With the aim of finding new molecular markers of failure of the TKI therapy we performed validation of the WES findings, transcriptome and miRNA analyses.
To validate these possible prognostic markers 62 patients with Ph+ CML in chronic phase were enrolled in the study. Median age was 46 years (range 19-71). Female:male ratio was 1:1.3. Fifty-nine patients received first-line therapy and threesecond-line due to intolerance to imatinib or sub-optimal response at 3 months. According to ELN 2013 recommendations [26] at 6 months of TKI treatment all patients were divided into two groupswith optimal (n = 32) and non-optimal (n = 30) response to targeted therapy. Description of genotyped patients is available in Table 1.
The results of the genotyping of all the patients are summarized in Table 2. The frequencies of genotypes at the studied loci did not differ between groups of patients with an optimal and non-optimal response to TKI therapy.
Since genotype rates for individual genes did not demonstrate significant differences between groups of patients with optimal and non-optimal response to TKI, we analyzed whether there is an association of response to therapy with haplotypes (different combinations of alleles in selected genes). There was no reliable association with the therapy response of any of the haplotypes presented in patients (p = 0.94). Correction by gender and age also did not change the result.

Transcriptome analysis
Thirty-three Ph+ CML patients with the follow-up period more than 6 months were treated in National Research Center for Hematology (Moscow). Most of the patients demonstrated good response to the TKI therapy. However, five didn't have optimal response. All of them had > 10% Ph+ cells at 6 months (according to ELN2013 recommendations non-optimal response). We also chose 7 patients with optimal response for the comparison. Patient ages ranged from 20 to 77 years, with a median of 37 years. All 12 patients were treated by TKIs: 10 patients treated by imatinib, 2 patients by nilotinib. Detailed characteristics of these patients are available in Table 3.
ANOVA analysis of the expression of 34,681 genes revealed 26 differently expressed genes (p < 0.05) in patients with optimal and non-optimal responses ( Table 4). Ten of these genes (DAZAP2, UBA52, PRR13, ATG7, LOC387820, PAK1, RAB11A, EMP3, GSTM1, GSTM2) are associated with other cancers and another four (RSAD2, MAP3K11, COMMD1, TNFRSF1A) are the key members of regulatory pathways including NF-kappa-B, MAPK, JNK and other thought to be disbalanced in many cancers (see Table 4). Imatinib 600 mg, n (%) 2 (6.7) Nilotinib 400-600 mg, n (%) *Sokal score data for 23 patients were not available **patients were switched to the second line therapy due to imatinib intolerance or sub-optimal response at 3 months TKI tyrosine kinase inhibitor Table 2 Rates of genotypes of analyzed genes in groups of patients with optimal and non-optimal response to TKI therapy Gene and reference SNP Reference allele Variant allele Patients with optimal response at 6 months (n = 32) Patients with non-optimal response at 6 months (n = 30) Pathway enrichment analysis of the differently expressed genes showed the following molecular networks involved (p < 0.05): glutathione metabolism, drug and xenobiotic metabolism, chemical mutagenesis, glycosaminoglycan degradation and MAPK signaling pathway. According to the expression levels we may suggest that these pathways are activated in patients who failed to achieve optimal response to TKI therapy.
From the list of the genes associated with cancers we selected first four with the largest difference in expression levels and analyzed their expression with RT-PCR: DAZAP2, UBA52, PRR13 and ATG7. Expression of these genes was measured in the same samples as whole transcriptome analysis and additional 10 samples from optimal responders were added. All the tested genes demonstrated similar differences in expression levels in optimal vice non-optimal responders as in whole transcriptome analysis (Table 4), however these differences were not significant (ANOVA, p = 0.3).

miRNA analysis
After transcriptome analysis we focused on miRNA expression in patients with different TKI efficacy. We had already 140 patients with the follow-up period more than 6 months treated in three Hospitals in Russia. According to ELN2013 recommendations we selected 6 patients with optimal response and 6 patients with failure of the TKI therapy. Patient ages ranged from 20 to 89 years, with a median of 52.5 years. All 12 patients were treated with 400 mg imatinib. Detailed characteristics of  Coefficients below zero correspond to elevated gene expression in optimal responders and coefficients above zero correspond to elevated gene expression in non-optimal responders these patients are available in Table 5. All patients with therapy failure had > 10% Ph+ cells at 6 months in bone marrow. Four patients were included in both transcriptome and miRNA analysis (#3, 5, 7, and 12 in Tables 3 and 5). ANOVA did not reveal significant difference between responders and failures. However since the number of subjects (samples) was very limited and the number of parameters (genes) was high it was reasonable to pay particular attention to the most significant expression changes of individual miRNAs. We selected miRNAs with individual p-value < 0.05-18 genes (Table 6). We analyzed these miRNAs using DIANA tools [28]. This set of programs allows pathway enrichment analysis with genes, interacting with selected miRNAs.
We found 40 pathways enriched with genes regulated by the selected miRNAs (Additional file 2: Table S2). These differently expressed miRNAs are involved in vital processes in cells including p53 signaling pathway, toll-like receptor signaling pathway, adherence junction, notch signaling pathway, PI3K-Akt signaling pathway, apoptosis, Wnt signaling pathway, MAPK signaling pathway, transcriptional misregulation in cancer and cell cycle. Among these pathways "Pathways in cancer" (hsa05200) has the highest p-value (9.5*10 − 14 ) with 27 genes regulated by 7 miRNAs. Among other cancer specific pathways one can particularly mention "Chronic myeloid leukemia" (hsa05220) (p = 0.00018) with 8 genes regulated by 4 miRNAs and "Acute myeloid leukemia" (hsa05221) (p = 0.04) with 4 genes regulated by 3 miRNAs.

Discussion
Here we describe our attempts to find molecular predictors of TKI therapy efficacy in primary CML patients based on genetic variants and expression.
We applied Sokal score to the patients included in transcriptome and miRNA analyses and were surprised by the poor prediction of the TKI efficacy. Seventy-three percent (8 from 11 patients) were correctly predicted as optimal responders and only 56% (5 from 9 patients) were correctly predicted as non-responders. Taking together only 65% of patients in our study were correctly classified based on Sokal score. Sure we clearly understand that our groups are very small for any conclusions about application of Sokal score. But we also understand Table 5 Characteristics of patients for miRNA study   Patient  Gender  Age, years  Sokal score  TKI first line, dosage  Efficient TKI   3  Female  57  Low risk  Imatinib, 400 mg  Optimal response   5  Female  68  Intermediate risk  Imatinib, 400 mg  Non-optimal response   7  Male  34  Intermediate risk  Imatinib, 400   Coefficients below zero correspond to elevation of miRNA expression in optimal responders and coefficients above zero correspond to elevated miRNA expression in non-optimal responders that it is very difficult to rely on the predictive value of this score for the prediction of early treatment efficacy. At the same time there were no patients with borderline values of BCR-ABL IS or Ph+ metaphases in our transcriptome and miRNA study. All responders really responded well to treatment, and non-responders were > 10% Ph+ cells at 6 months in bone marrow. Therefore, we conducted this study in order to find new molecular predictors of treatment efficacy . To date, researchers performed several attempts to detect genetic markers associated with resistance to TKI or non-optimal response to targeted therapy. It was shown that the resistance of leukemic cells is associated with the enhanced PI3K/AKT, RAF/MEK/ERK and STAT3 signaling pathways [29]. At the same time, studies involving patients with different responses to TKI, show influence of genetic polymorphisms in TKI metabolizers' genes [9], BIM gene [30], FAS gene [10], DNMT3A and ASXL1 genes [11].
We did not confirm different allelic frequencies in patients with different response to the therapy. This can be explained by the fact that the efficiency of TKI therapy can be affected by a variety of both genetic and environmental factors, making the therapy failure a multifactorial feature.
Previously performed studies showed that the expression profiles in patients with optimal and non-optimal responses to TKI therapy varied widely. For example, as published by Frank et al. in 2006, differential expression of some genes was found in two groups of patients and based on this findings authors made conclusions about the discovery of candidate predictors of imatinib resistance [14]. In another study, published in the same year in the same journal, significant difference in the expression of other genes was showed [15]. Finally, one more work in the same journal demonstrated a third list of genes of potential predictors [31]. All these works were done before the next-generation sequencing era, so we conducted our study to find the expression differences by unbiased analysis of the large number of genes. We found some changes in expression levels, but they were very small and it is inappropriate to further look for changes in some specific genes and to develop prognostic tests on their basis. Additionally, we checked all genes that have been previously mentioned in similar studies to be associated with different responses to TKI therapy. We found no difference in allelic frequencies of these genes in our groups. It is possible that revealing specific activation patterns of the whole pathogenic mechanisms is more promising for diagnostic purposes. Another possible solution is to analyze expression of bone marrow samples, containing more stem cells or direct selection and analysis of cancer stem cells which is however difficult for routine implementation.
We revealed 17 differently expressed miRNAs in patients with optimal and non-optimal responses to TKI therapy. For nowadays miRNA profiling is frequently used to find molecular predictors, treatment targets or biomarkers in cancer. Indeed miRNAs are involved in CML pathogenesis [17][18][19]. Recently it was shown that miR-92a-3p functions as oncogene blocking apoptosis in AML cells. Inhibition of this miRNA leads to increase of apoptosis in M-07e cell line [32]. Nevertheless in our study patients with optimal response had 8-fold higher expression of miR-92a-3p than non-optimal responders. We revealed 3-fold decrease of miR-331-5p in treatment failures. This finding is in the agreement with previously reported association of low levels of this miRNA in K562 cells with the resistance to doxorubicin [33]. We also found that expression of miR-324-5p was 3.8-fold higher in patients with non-optimal response to TKI which correlates with its elevated expression in patients with poor outcome (progression-free and overall survival) of the diffuse large B-cell lymphoma [34]. Despite such large differences in miRNA expression levels between groups, we did not continue this work due to the absence of significant difference using ANOVA with strict corrections for multiple comparisons.
We didn't sequence BCR/ABL kinase domain in all patients, because these mutations usually develop during treatment and result in the secondary resistance to TKIs while we examined patients with primary TKI resistance. Nevertheless we checked 7 out of 11 patients included in transcriptome and miRNA expression analysis and didn't find mutations in BCR/ABL kinase domain (data not shown).

Conclusions
In conclusion, using modern high-throughput methods such as whole-exome, transcriptome and miRNA analysis, we could not find reliable molecular markers for early prediction of TKI efficiency in Ph+ CML patients. Apparently, efficacy of targeted therapy is determined not only by genetic factors, the environmental component of therapy response also plays a great role. It is necessary to continue search for mechanisms of tumor cells sensitivity to therapy for their complete elimination.