Molecular Changes in Chronic Myeloid Leukemia During Tyrosine Kinase Inhibitors Treatment. Focus on Immunological Pathways.

Introduction: The aim of our research was to investigate changes in the molecular background of the immune response in the chronic phase (CP) of chronic myeloid leukaemia (CML) during treatment with tyrosine kinase inhibitors (TKIs).
Methods: Global gene and miRNA expression profiles were assessed using genome-wide RNA and miRNA microarray technology in bone marrow mononuclear cells. Fifty-one patients were recruited, and bone marrow samples were taken at diagnosis before treatment with TKIs and after 3, 6, and 12 months of treatment with TKIs. The largest number of upregulated genes was observed when the 0-month group (time of diagnosis) was compared to the 3-month group; 1774 genes were significantly upregulated, and 390 genes were significantly downregulated. Discussion: Upregulated biological processes according to gene ontology (GO) classification involved basic cellular processes such as cell division, cell cycle, cell-cell adhesion, protein transport, mitotic nuclear division, apoptosis, and DNA replication. Differentially expressed miRNAs were annotated using GO classification to several immunity-related processes, including the T cell receptor signalling pathway, T cell costimulation, immune response, and inflammatory response. TKI therapy exerts a significant impact on cellular cycle processes and T-cell activation, which was proven at the molecular level.

Introduction

Chronic myeloid leukaemia (CML) is a clonal myeloproliferative disorder that arises after a balanced reciprocal translocation involving chromosomes 9 and 22, t(9; 22) (q34; q11).1 This rearrangement produces the BCR-ABL1 transcript and consequently oncoprotein, which acts as a constitutively expressed defective tyrosine kinase. The downstream pathways affected include JAK/STAT, PI3K/AKT, and RAS/MEK, which involve cell growth, cell survival, inhibition of apoptosis, and activation of transcription factors.2 The BCR-ABL1 oncoprotein is a target for tyrosine kinase inhibitors (TKIs). The BCR-ABL1 transcript is a precise parameter for evaluating treatment response. The most common BCR-ABL1 transcripts are the e13a2 and e14a2 transcripts.3,4 In the vast majority of patients, treatment with TKIs results in suppression of leukaemic clone proliferation and achievement of complete haematologic remission with restoration of the proportion of the cell population present in the blood. This revolutionary targeted therapy has significantly improved the survival rates of patients affected with CML, which are now comparable to those of age-matched healthy individuals.5 However, the risk of a refractory chronic phase or even progression to an accelerated phase and blast crisis of CML still occurs and remains a serious clinical challenge. Despite attempts at safe therapy cessation, called treatment-free remission (TFR), which is now the therapeutic goal, the disease remains incurable.6 Considering the possibility of discontinuing treatment with TKIs without relapse and the phenomenon of a so-called “functional cure”, much attention has been given to the role of the immune system in CML pathophysiology. The immune system seems to be responsible for the control of residual leukaemic cells and the lack of relapse of CML in patients without treatment. They exhibit low levels of BCR-ABL1 transcript as measured by highly sensitive BCR-ABL DNA detection methods.7 Since leukaemic stem cell LSCs are TKI resistant and genetically unstable, they are a potential target for new immunotherapeutic approaches.8 The important role of immune surveillance over this population of TKI-resistant cells in maintaining a profound molecular response after cessation of therapy has been postulated. To date, cell-mediated immunity in the control of leukaemic cells patients has been most extensively studied.9 Considerably less is known about molecular changes in response to TKI treatment, especially the expression of genes and microRNAs. MicroRNAs (miRNAs, miRs), 20–22 nucleotide long noncoding RNA molecules, have recently emerged as key regulators of various genes, and their role in immune regulation has also been suggested.10 MiRNAs regulate gene expression at the posttranscriptional level and are involved in almost every biological process as well as neoplasm-associated processes. Recently, the role of miRNAs in the progression of CML and the development of TKI resistance, based on their important regulatory function in cell homeostasis, has been implicated.11 Differentially expressed miRNAs have been reported in several diseases and as a response to treatment.12 To better understand the underlying molecular mechanisms of immune restoration during TKI treatment, in the present study, we investigated global mRNA and miRNA expression in bone marrow cells in CML patients at different time points during the first 12 months of therapy with imatinib (IM). To date, a wide panel of mRNAs/miRNAs has not been tested at different time points of CML treatment to provide insight into the immune changes and treatment response. The aim of this study was to elucidate possible correlations between miRNAs and their target gene expression profiles in CML patients treated with TKIs at the time of diagnosis compared with several time points during the first year of therapy.

Materials and Methods

The Study Cohort

The study was conducted on 51 patients with CML enrolled from the Department of Haematology and Bone Marrow Transplantation, Pomeranian Medical University, Szczecin, Poland. Prior to enrolment, all patients signed informed consent in accordance with the Declaration of Helsinki. Before patient enrolment, we obtained appropriate approval from the local ethics committee of the Pomeranian Medical University (approval code: KB-0012/11/2021). We enrolled adult patients in the chronic phase of CML. Consecutive patients in whom a diagnosis of chronic phase CML was established were recruited for the study. The inclusion criteria for the study were as follows: age from 18 to 80 years, diagnosed chronic myeloid leukemia in the chronic phase. The exclusion criteria were as follows: age under 18 or over 80 years, diagnosis chronic myeloid leukaemia in accelerated or blastic phase, coexistence of another myeloid or lymphoid neoplasm or suspected secondary chronic leukaemia to prior conducted therapy (eg radiotherapy). The median age of the patients was 54.5. Patients received IM (76%) or second-generation TKIs as a second-line treatment (24%, nilotinib=3, dasatinib=9). All patients started IM therapy, and some patients changed treatment. The study cohort is summarized in Table 1.

Table 1

Characteristics of Study Cohort

Patients’ Cohort		N=51
Demographic features - Sex
Female		24 (47)
Male		27 (53)
Clinical features – age of diagnosis (years)
Median		54.5
Range		26–79
Sokal Score at diagnosis
Low	Intermediate		High
30	17		4
Treatment
Imatinib (%)		39 (76)
2nd TKI (%)		12(24)
Response
After 3 months
Optimal response		Without optimal response
28		30
After 6 months
Optimal response		Without optimal response
30		21
After 12 months
Optimal response		Without optimal response
20		31

To assess the effect of TKI treatment on the molecular mechanisms regulating the immune response, we collected bone marrow (BM) samples at four consecutive time points. BM was aspirated at 0, 3, 6 and 12 months from the posterior iliac crest from diagnosis. Next, the BM sample was diluted 1:1 with phosphate-buffered saline (PBS) and layered on Lymphocyte Separation Medium (MP Biomedicals, Santa Ana, CA, USA). The obtained mononuclear cell suspension was then subjected to RNA and miRNA isolation.

RNA and miRNA Isolation

Total RNA enriched in miRNAs was isolated from bone marrow mononuclear cells (1 × 106) using the mirVana™ miRNA Isolation Kit (Thermo Fisher, Waltham, MA, USA) following the manufacturer’s instructions. An Epoch spectrophotometer (Biotek, Winooski, VT, USA) was used to assess the concentration and quality of the obtained RNA. For subsequent miRNA and whole transcriptome microarray analysis, total RNA enriched in miRNAs isolated from bone marrow samples from three patients was pooled to generate one sample per group.

Affymetrix GeneChip Microarrays

RNA and miRNAs isolated from BM mononuclear cells of patients at different time points were pooled to generate one sample per group for microarray experiments. Affymetrix miRNA 4.1 Array Strip (Affymetrix, Santa Clara, CA) and Affymetrix Human Gene 2.1 ST Array Strip (Affymetrix) was used for the study. Generation of sense strand cDNA from the total RNA and subsequent fragmentation and labelling steps were performed using the GeneChip™ WT PLUS Reagent Kit (Thermo Fisher Scientific, Waltham, MA, USA). Finally, the sample was hybridized onto an Affymetrix Human Gene 2.1 ST Array Strip. The Affymetrix GeneAtlas System was used for hybridization, fluidics and scanning steps. Subsequent analyses were performed using BioConductor software. We used the Robust Multiarray Average (RMA) normalization algorithm (from the “Affy” library) for normalization, background correction, and calculation of the expression levels of the examined genes. Microarrays were made and analysed according to the methods described in previous studies.37,38

DAVID

For functional annotation and enrichment analysis, we used DAVID Bioinformatics Resources (Database for Annotation, Visualization, and Integrated Discovery) at , as previously described 25, 26. Functional annotation charts generated by DAVID with overrepresented gene annotations are shown as bubble plots from the BACA BioConductor package (). The following criteria were applied to generate bubble plots: p value < 0.5, adjusted method = Benjamini 27, and minimal number of genes per group = 5. Groups of genes that met these criteria are shown in a graph, where the bubble size is indicative of the number of genes represented in the corresponding annotation and their down- or upregulation.

Statistical Analysis

The nonparametric Mann–Whitney test was used to compare values between study groups, because in the majority of cases the quantitative variable distribution notably differed from the normal distribution. Spearman’s rank correlation coefficient (Rs) was used to measure the strength of associations between gestational age and concentrations of each tested complement component. The sign of the Rs value indicates the direction of the association (positive or negative), while a higher absolute value (closer to -1 or +1) indicates a stronger association. We considered p<0.05 statistically significant. Statistica 13 software (Dell Inc., OK, USA) was used for statistical analysis.

Results

Gene Expression Profile in Mononuclear Cells Derived from Bone Marrow Aspirates

We distinguished three separate groups based on the time from CML diagnosis (3, 6, 12 months) and compared their gene expression profiles with the group at diagnosis (“Group 0”). Microarray analysis revealed that the highest number of differentially expressed genes in this comparison was observable for the 0-month group when compared to the 3-month group. In this case, 1774 genes were at least 2-fold upregulated and 390 genes were at least 2-fold downregulated in the 0-month group compared to the 3-month group (Figure 1). When we analysed gene expression profiles in the 6- and 12-month groups, the number of genes that were differentially expressed was smaller. Figure 2 presents 211 upregulated and 152 downregulated genes in the 0-month group compared to the 6-month group. Figure 3 presents 304 upregulated and 138 downregulated genes in the 0-month group compared to the 12-month group.

Figure 1

The scatterplot of global gene expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 3 months after diagnosis (HP3). Red points show downregulated genes (at least 2-fold change, p < 0.05), green points show upregulated genes (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of genes with the largest change in expression.

Figure 2

The scatterplot of global gene expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 6 months after diagnosis (HP6). Red points show downregulated genes (at least 2-fold change, p < 0.05), green points show upregulated genes (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of genes with the largest change in expression.

Figure 3

The scatterplot of global gene expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 12 months after diagnosis (HP12). Red points show downregulated genes (at least 2-fold change, p < 0.05), green points show upregulated genes (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of genes with the largest change in expression.

Next, we looked at the most significantly up- and downregulated genes in the 0-month group compared to the 3-, 6- and 12-month groups (Tables 2–7). Interestingly, several genes showed similar expression profiles. For example, HIST2H3D was significantly upregulated in all analyses; however, the highest fold change (FC) equal to 14.30 was observed when the 0-month group was compared to the 3-month group (FC 4, 37 for 0 vs 6 months and 4, 82 for 0 vs 12 months). On the other hand, the expression of FAM95B1, CXCL8 and SLC28A3 was significantly upregulated only when the 0-month group was compared to both the 6- and 12-month groups. The lists of the ten most downregulated genes when the 0-month group was compared to the 6- and 12-month groups were distinctly similar, as 9 out of 10 genes were the same and differed only slightly in their FCs: IGHV3-23, IGKJ3, TRBC2, TRAJ35, IGLV3-1, IGLV1-40, HLA-DRB1, HLA-DRB5, and IGHJ6 (Tables 5 and 7). Four of these genes were also significantly downregulated in the 0-month group compared to the 3-month group: IGLV3-1, IGLV1-40, TRBC2, and TRAJ35 (Table 3). In Table 8, we summarized genes whose expression changed similarly at diagnosis compared to each subsequent time point (3, 6 and 12 months).

Table 2

The List of 10 Most Upregulated Genes in 0-Month Group in Comparison to 3-Month Group

Gene Symbol	Gene Name	Fold Change
XIST	X inactive specific transcript (non-protein coding)	25.11
HIST1H3A	Histone cluster 1, H3a	16.69
SERPINB10	Serpin peptidase inhibitor, clade B (ovalbumin), member 10	14.63
HIST2H3D	Histone cluster 2, H3d	14.30
CYP4F3	Cytochrome P450, family 4, subfamily F, polypeptide 3	11.46
ANKRD22	Ankyrin repeat domain 22	11.00
CEACAM4	Carcinoembryonic antigen-related cell adhesion molecule 4	10.13
HIST2H2BE	Histone cluster 2, H2be	9.39
HIST1H2BF	Histone cluster 1, H2bf	9.21
TRGJP1	T cell receptor gamma joining P1	9.14

Table 3

The List of 10 Most Downregulated Genes in 0-Month Group in Comparison to 3-Month Group

Gene Symbol	Gene Name	Fold Change
GIMAP7	GTPase, IMAP family member 7	−7.90
CD79A	CD79a molecule, immunoglobulin-associated alpha	−8.33
FAIM3	Fas apoptotic inhibitory molecule 3	−8.82
IGHM	Immunoglobulin heavy constant MU	−11.18
IL7R	Interleukin 7 receptor	−12.88
IGLV3-1	Immunoglobulin lambda variable 3–1	−13.27
IGHG1	Immunoglobulin heavy constant gamma 1 (G1m marker)	−15.53
IGLV1-40	Immunoglobulin lambda variable 1–40	−16.27
TRBC2	T cell receptor beta constant 2	−18.46
TRAJ35	T cell receptor alpha joining 35 (non-functional)	−25.63

Table 4

The List of 10 Most Upregulated Genes in 0-Month Group in Comparison to 6-Month Group

Gene Symbol	Gene Name	Fold Change
RP11-146D12.2	Novel protein	6.12
FAM95B1	Family with sequence similarity 95, member B1	5.44
APOLD1	Apolipoprotein L domain containing 1	5.06
CXCL8	Chemokine (C-X-C motif) ligand 8	4.85
CCDC144A	Coiled-coil domain containing 144A	4.81
HIST1H2BH	Histone cluster 1, H2bh	4.67
HIST2H3D	Histone cluster 2, H3d	4.37
KIAA1524	KIAA1524	4.29
SLC28A3	Solute carrier family 28 (concentrative nucleoside transporter), member 3	4.24
ANKRD36B	Ankyrin repeat domain 36B	4.17

Table 5

The List of 10 Most Downregulated Genes in 0-Month Group in Comparison to 6-Month Group

Gene Symbol	Gene Name	Fold Change
IGLV2-23	Immunoglobulin lambda variable 2–23	−8.26
IGHV3-23	Immunoglobulin heavy variable 3–23	−9.41
IGKJ3	Immunoglobulin kappa joining 3	−10.12
TRBC2	T cell receptor beta constant 2	−10.19
TRAJ35	T cell receptor alpha joining 35 (non-functional)	−11.04
IGLV3-1	Immunoglobulin lambda variable 3–1	−11.27
IGLV1-40	Immunoglobulin lambda variable 1–40	−24.00
HLA-DRB1	Major histocompatibility complex, class II, DR beta 1	−30.04
HLA-DRB5	Major histocompatibility complex, class II, DR beta 5	−40.02
IGHJ6	Immunoglobulin heavy joining 6	−89.13

Table 6

The List of 10 Most Upregulated Genes in 0-Month Group in Comparison to 12-Month Group

Gene Symbol	Gene Name	Fold Change
JUN	Jun proto-oncogene	9.99
CXCL8	Chemokine (C-X-C motif) ligand 8	7.05
FAM95B1	Family with sequence similarity 95, member B1	5.87
EGR1	Early growth response 1	5.58
HIST2H3D	Histone cluster 2, H3d	4.82
SNORA52	Small nucleolar RNA, H/ACA box 52	4.65
SLC28A3	Solute carrier family 28 (concentrative nucleoside transporter), member 3	4.64
PMAIP1	Phorbol-12-myristate-13-acetate-induced protein 1	4.59
IL23A	Interleukin 23, alpha subunit p19	4.28
FOSB	FBJ murine osteosarcoma viral oncogene homolog B	3.95

Table 7

The List of 10 Most Downregulated Genes in 0-Month Group in Comparison to 12-Month Group

Gene Symbol	Gene Name	Fold Change
IGHJ3	Immunoglobulin heavy joining 3	−8.39
IGLV3-1	Immunoglobulin lambda variable 3–1	−9.04
IGKJ3	Immunoglobulin kappa joining 3	−9.12
TRBC2	T cell receptor beta constant 2	−9.14
IGHV3-23	Immunoglobulin heavy variable 3–23	−9.21
TRAJ35	T cell receptor alpha joining 35 (non-functional)	−13.75
HLA-DRB5	Major histocompatibility complex, class II, DR beta 5	−13.82
IGLV1-40	Immunoglobulin lambda variable 1–40	−20.60
HLA-DRB1	Major histocompatibility complex, class II, DR beta 1	−27.44
IGHJ6	Immunoglobulin heavy joining 6	−196.35

Table 8

Summary of Genes Whose Expression Changed Significantly and Similarly at Diagnosis Compared to Each Subsequent Time Point (3, 6 and 12 Months)

Gene Symbol	Upregulated Genes
	Fold Change: 3 Month	Fold Change: 6 Month	Fold Change: 12 Month
HIST2H3D	14.30	4.37	4.82
Gene Symbol	Downregulated Genes
	Fold Change: 3 Month	Fold Change: 6 Month	Fold Change: 12 Month
IGLV3-1	−13.27	−11.27	−9.04
IGLV1-40	−16.27	−24.00	−20.60
TRBC2	−18.46	−10.19	−9.14
TRAJ35	−25.63	−11.04	−13.75
IGHV3-23	—	−9.41	−9.21
IGKJ3	—	−10.12	−9.12
HLA-DRB1	—	−30.04	−27.44
HLA-DRB5	—	−40.02	−13.82
IGHJ6	—	−89.13	−196.35

Since these significantly up- and downregulated genes encode proteins important to the immune system (immunoglobulins, MHC, T cell receptors), we aimed to take a broader look at the gene expression profiles in our patients and investigate whether immunity-related processes were dysregulated. Therefore, in the next part of our analysis, the differentially expressed genes were classified according to the gene ontology (GO) classification of biological processes (Figure 4). The largest number of upregulated processes was observed when the 0-month group was compared to the 3-month group, which corresponds to the high number of differentially expressed genes between these groups, as discussed earlier (Figure 1). These upregulated processes mainly included basic cellular processes: cell division, cell cycle, cell-cell adhesion, protein transport, mitotic nuclear division, autophagy, apoptosis and DNA replication. Next, we analysed the dynamics of changes in different groups of genes associated with biological processes. The changes in gene expression and their biological processes at the 6th month of treatment revealed a similar trend as at the 3rd month. There was a noticeable upregulation of genes involved in regulation of the cell cycle, DNA repair processes. The t-cell signalling pathway and regulation of the adaptive immune response were the biological processes most altered at the 12-month time point. The majority of downregulated processes in all of the comparisons were related to immunity: T cell receptor signalling pathway, regulation of immune response, T cell costimulation, immune response and adaptive immune response. Interestingly, the number of these downregulated processes and genes involved was the highest when the 0-month and 3-month groups were compared and decreased in subsequent time points.

Figure 4

The bubble plot of overrepresented biological processes according to gene ontology (GO) classification in 0-moths group (HP0) compared to 3- (HP3), 6- (HP6) and 12-month (HP12) groups. Genes assigned to individual processes fulfilling the criteria: adjusted p < 0.05, method = Benjamini, and minimum number of genes per group = 5, are presented. The bubble size indicates the number of genes represented in the corresponding annotation. The up- and downregulation is indicated by colors, with red corresponding to upregulation and blue corresponding to downregulation.

Next, we implemented the obtained results to assess global shifts in gene expression at the subsequent time points analysed. The K-clustering approach distinguished five clusters of significantly changed genes in four analysed groups (0, 3, 6, 12 months). This approach enabled us to analyse groups of genes similarly expressed at subsequent time points. As shown in Figure 5, genes grouped in Cluster 1 and 4 showed the opposite trend in expression, whereas Cluster 2 and 3 were similar and opposite to Cluster 5.

Figure 5

Changes in the expression of genes assigned to each of the 5 clusters at each time point: 0 month (HP0), 3 month (HP03), 6 month (HP06), 12 month (HP12). On the right, heatmaps representing the expression of genes assigned to each gene cluster in all analyzed time points.

The assignment of differentially expressed genes to relevant Gene Ontology terms revealed that Cluster 2 and 3, which shared a similar pattern of expression at subsequent time points, also contained genes involved in related processes: cell division, regulation of the cell cycle, cell proliferation and protein phosphorylation (Figure 6). The cluster that grouped genes with increasing expression in subsequent points, Cluster 1, contained genes involved in processes related to immunity: T cell receptor signaling, T cell activation, T cell costimulation, antigen processing and presentation.

Figure 6

The bubble plot depicts the assignment of genes in each cluster to biological processes based on gene ontology (GO) classification. Groups of genes fulfilling criteria: adjusted p < 0.05, method = Benjamini, and minimum number of genes per group = 5, are presented. Bubble size indicates the number of genes represented in corresponding annotation.

MiRNA Expression Profile in Mononuclear Cells Derived from Bone Marrow Aspirates

The results of the miRNA expression analysis showed a similar number of differentially expressed genes (at least 2-fold) in the group at diagnosis compared to the 3-, 6-, 12-month group: 146 significantly upregulated and 32 downregulated miRNAs compared to the 3-month group (Figure 7), 116 significantly upregulated and 31 downregulated miRNAs compared to the 6-month group (Figure 8) and 91 significantly upregulated and 39 downregulated miRNAs compared to the 12-month group (Figure 9).

Figure 7

The scatterplot of miRNAs expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 3 months after diagnosis (HP3). Red points show downregulated miRNAs (at least 2-fold change, p < 0.05), green points show upregulated miRNAs (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of miRNAs with the largest change in expression.

Figure 8

The scatterplot of miRNAs expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 6 months after diagnosis (HP6). Red points show downregulated miRNAs (at least 2-fold change, p < 0.05), green points show upregulated miRNAs (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of miRNAs with the largest change in expression.

Figure 9

The scatterplot of miRNAs expression in mononuclear cells derived from bone marrow aspirates at the diagnosis (HP0) compared to 12 months after diagnosis (HP12). Red points show downregulated miRNAs (at least 2-fold change, p < 0.05), green points show upregulated miRNAs (at least 2-fold change, p < 0.05). The grey dotted lines mark the border of significantly higher or significantly lower expression (log2(fold)>2), whereas the grey dashed lines correspond to log2(fold)>4. The graph also comprises the names of miRNAs with the largest change in expression.

Since miRNAs regulate gene expression, we decided to conduct an analysis comparing the outcomes of miRNA microarray expression data to the obtained mRNA data, with a particular focus on immunity-related processes. Figure 10 presents a network of connections between miRNAs and their target mRNAs, all involved in T cell activation.

Figure 10

The diagram shows significantly up- or downregulated miRNAs and their target genes involved in T cell activation in 0-month group when compared to 3-month group. Target genes are assigned to each miRNA with a marked change in expression (at least fold 2 or −2) (red indicates upregulation and green indicates downregulation). The solid lines correspond to the validated genes, and the dashed lines corresponds to the predicted genes.

Discussion

The latest achievements in the treatment of CML allow for the discontinuation of pharmacotherapy and the maintenance of a deep molecular response. It is thought that recovery from CML is the consequence of two parallel processes: halting the transduction of proliferative signals and renewing the immune response. We hypothesized that TKIs also play an important role also in this second process. The first step of our study was to analyse the global gene expression in BM cells from CML patients at different time points during TKI treatment and to compare it with the time of diagnosis. The most pronounced gene expression upregulation was observed in comparison to 0-month and 3 months after introducing treatment with TKI. The gene with the highest fold change, equal to 25.11 in this comparison was XIST, a long noncoding RNA (lncRNA) essential for the initiation of X chromosome inactivation (XCI), which in turn plays an important role in silencing one X chromosome in female mammals to balance gene expression between the sexes.13 Selective excision of XIST in murine haematopoietic stem cells (HSCs) leads to progressive myelofibrosis, proliferation of granulocytic lineage, myelodysplastic changes and appearance of extramedullary haematopoiesis revealed as splenomegaly as presented by Yildirim et al.14 The changes in haematopoiesis accompanying the knockdown of XIST are consistent with the clinical presentation of myeloproliferative neoplasms, including CML. Our study provides quite different data and we demonstrated XIST downregulation three months after treatment with TKIs. This is the moment of the most significant recovery of normal haematopoiesis, which is particularly evident in the first months of treatment and is related to the achievement of a complete haematological response. Wang et al15 proved that XIST was highly expressed in acute myeloid leukaemia (AML) bone marrow cells and that silencing XIST could reduce the drug resistance of AML bone marrow cells via downregulation of MYC. However, in light of the aforementioned data, the role of XIST expression in bone marrow niche is not clear; it is noteworthy, that XIST is overexpressed in gliomas compared to normal brain tissues and is related to miR-204-5p downregulation.16 In other human cancers, such as non-small cell lung, cancer breast cancer, pancreatic and gastric cancer, XIST also acts as an oncogene.17 Conversely, in breast cancer, XIST has been found to act as a tumour-suppressor. In this context, our study provided evidence against the hypothesis developed by Yildirim et al regarding the suppressive role of XIST in haematologic malignancies. Here, we demonstrated for the first time that X–inactive specific transcripts were significantly downregulated after introduction of TKI treatment in CML patients. In further analyses compared to diagnosis, the fold change in XIST expression was as follows - at 6-months, XIST expression increased with a fold 0.73 and at 12-month XIST expression slightly decreased with a fold change of -1, 33. We suggest that the potential clinical implications of the results we presented are related to the high level of XIST gene expression at diagnosis and its reduction following treatment. Given the observations of other research teams that have demonstrated that XIST inhibition can reduce the viability and drug resistance of acute myeloid leukemia cells by reducing MYC expression and increasing miR-29a expression.15 Strategies leading to decreased expression of this gene may improve the efficacy of new therapies.

Another significantly upregulated gene at the moment of diagnosis was HIST1H3A (histone cluster 1, H3a gene). Its role in the pathomechanism of CML is unknown, but it seems to play a role in the progression of Fanconi anaemia to acute myeloid leukaemia.18 Of the ten most downregulated genes in the 0-month group compared to the 3-month group, a significant proportion of genes were responsible for forming immunoglobulins. This downregulation trend deepened 6 and 12 months after the initiation of TKI treatment. Our results obtained at the molecular level contradicted clinical observations in which imatinib-treated patients showed a significant decrease in IgA and IgG levels, and dasatinib-treated patients showed decreased IgM levels.19 Then, using a clustering approach, we were able to detect groups of genes that were significantly changed at the assessed time points. The analysis revealed five clusters of differentially expressed genes (DEGs) in all four analysed groups included in the analysis (0, 3, 6 and 12-months). The two clusters presented a similar pattern of dynamics of change over time in patient treatment (Clusters 2 and 3). These clusters comprised genes assigned to corresponding processes involved in cell division, regulation of the cell cycle, cell proliferation and protein phosphorylation. Patients treated with TKIs had significantly lower expression of these genes after 3 months of treatment, resulting in apparent suppression of these processes. Considering the clinical features, this is the period where the most significant reduction of leukocytosis occurs, to achieve a complete haematological response in the vast majority of patients. It is also the time where the most pronounced hematologic toxicity of TKIs is observed, expressed by cytopenias visible in the peripheral blood. At months 6 and 12 of therapy, the expression of these genes increased and reached a level of approximately half the activity compared to that seen at the time of diagnosis. Cluster 1 contained genes involved in processes related to cell immunity: T cell receptor signalling ~GO:0050852, T cell activation ~GO:0042110, T cell costimulation ~GO:0031295, antigen processing and presentation and signalling pathway mediated by interferon gamma ~GO:006033. Interestingly, we observed an increase in the expression of genes grouped in Cluster 1 after 3 months of TKI introduction. The activation of signalling pathways responsible for immunity related to T cells is not reflected in changes in the populations of these lymphocytes. Lu et al demonstrated that during treatment with TKIs, the absolute number of total T cells and CD4+ T and CD8+ T cells was decreased compared with the time of diagnosis.20 Analysis of genes classified in Cluster 5 in the context of biological processes in GO revealed changes in cell immunity mediated mainly by CD8+ T cells and NK cells. In Cluster 5, we observed a renewal of the pathways responsible for T cell receptor signalling ~GO:0050852, positive regulation of natural killer cell mediated toxicity ~GO:0045086, and T cell costimulation ~GO:0031295. The gene expression landscape is consistent with the observations of other authors, as briefly discussed below. It has been well proven that programmed cell death protein 1/programmed death ligand-1 (PD-1/PD-L1), an immune checkpoint regulator, plays a significant negative role, inhibiting the activation, anergy, and eventually causing the apoptosis of T cells. TKIs decrease PD-1 expression levels on CD8+ T cells. PD-1 expression in CD8+ T cells is lower at the time of achieving a complete haematological response than in patients in the chronic phase of CML at diagnosis.21 Importantly, in the aforementioned study, a significant change in PD-1 expression occurred when patients achieved complete haematologic response (CHR). CML-specific cytotoxic T lymphocytes are another population of lymphocytes that increase in parallel to the depth of the molecular response.22 Likewise, the analysis of the expression of analogous signaling pathways in our patient group indicated a similar interaction. Our observations indicate that in Clusters 1 and 5, which contain genes responsible for the immune system, the greatest dynamics of changes occurred during the first three months of treatment. In subsequent observation, gene expression either reached a plateau phase as in Cluster 1 or decreased as in Cluster 5. On this basis, we conclude that the first period of TKI use is the most important for the achievement of immunological effects, and the continued use of TKIs only enhances this effect. Perhaps this observation will be helpful in further research related to TFR. NK cells constitute a significant component of innate immunity, acting against malignant cells by recognizing the downregulation of HLA-class I molecules.23 In vitro studies with immune cell exposure to different TKIs have suggested their immunosuppressive effects on T cell and NK-cell activation.24 In another study, it was demonstrated that all TKIs exert an inhibitory effect on CD56+ cell restoration, and dasatinib reduces the cytolytic activity of NK cells.25 Hughes et al demonstrated an increased number of NK cells with a more mature phenotype, an increased number of CML-specific cytotoxic T lymphocytes, and reduced immune suppressor myeloid-derived suppressor cells (MDSCs) and regulatory T cells (Tregs) in TKI-treated patients compared to time at diagnosis.22 We suggest that the increased expression of genes involved in the positive regulation of NK cell-mediated immunity is involved in the maturation of the phenotype of NK cells during TKI treatment, as previously demonstrated by Hughes.22 An increasing NK cell population is observed in the bone marrow niche during treatment and impacts the clinical response.26 On the other hand, a higher number of NK cells at the moment of TKI discontinuation is a predictive factor for maintaining a response despite treatment withdrawal.24 Similarly, positive changes in humoral immunity were observed in the form of an upregulation of the interferon-gamma production pathway in Cluster 5. As in other clusters, the most important changes in gene expression occurred at 3 months of TKI treatment, and the dynamics of changes in the expression of these genes in the following months were substantially lower. Emerging reports about the combination of TKIs with interferon-alpha suggest not only different sites of action of these drugs but also additional immunomodulatory effects of interferon.27 It has been well documented that interferon pretreatment is significantly associated with the success of imatinib cessation.28 This observation corresponds to the hypothesis of regaining control by the immune system. In the future, it may also lead to the control of a residual clone of leukaemic cells. In summary, in the analysis of global gene expression, we reviewed studies in this field, and we did not find studies with a similar design. We found that Radich et al29 identified genes that can influence CML progression and annotated HIST1H2AE as a “progression” gene. This observation is consistent with ours. We found deregulation of several histone-coding genes, such as HIST1H3A, HIST2H3D, HIST2H2BE, and HIST1H2BF, during treatment with TKIs.

In contrast, the current knowledge of humoral immune response-related molecular events in CML is very limited. MiRNAs, which are 20–22 nucleotide long noncoding RNA fragments, have an impressive regulatory impact on gene expression. The search for markers of various malignancies among miRNAs has continued since their discovery in the 1990s, and it has brought encouraging results.10 MiRNAs that have fundamental regulatory impacts on both innate and adaptive immune cells in health and disease are called “immune-miRs”.10 Immune miRs regulate immune processes as well as central inflammatory signalling pathways. Thus, we analysed the global expression of miRNAs in BM cells from CML patients at different time points during TKI treatment. The miRNA with the highest fold change at diagnosis compared to 3 months of treatment was miR-361-5p. This molecule suppresses angiogenesis by reducing the expression of VEGF-A, and it modulates the inflammatory response by downregulating IL-6 and IL-8 expression.30 Both mechanisms are involved in the pathogenesis of CML; thus, downregulation of miR-361-5p positively correlated with CML treatment. The expression fold change of miR-145-5p was 17 times greater in diagnosis than after 3 months of treatment. This miRNA silences Semaphorin 3A, which acts as a terminator of T and B cells. MiR-145-5p decreases the concentration of proinflammatory cytokines, and this cascade is observed in autoimmune diseases.31 It has also been shown that radiation injury of BM upregulates the expression of the gene encoding Semaphorin 3A, which can negatively regulate vascular recovery following radiation injury.32 In our study, the miRNA expression profile of BM cells revealed significant upregulation of miR-652-3p in naïve CML patients. A similar tendency was observed in bladder cancer tissue, where miR-652-3p promotes the proliferation, migration, and invasion of cancer cells.33 The opposite trend was observed in a cohort of paediatric patients treated for acute lymphoblastic leukaemia (ALL).34 In that study, it was shown that miR-652-3p was downregulated at diagnosis compared with healthy controls, and the level of miR652-3p was restored in complete remission of ALL, negatively correlating with disease severity. However, the BCR/ABL1 rearrangement was rare in the study group and accounted for 2.3%. It was suggested that miR-652-3p may be a potentially useful biomarker for ALL in children, but the clinical usefulness of this miRNA in CML remains to be elucidated. Flamant et al observed increased expression of miR-150 after two weeks of imatinib therapy; the results of our study provided evidence for decreased expression of miR-150 at the time of diagnosis compared to 3, 6 and 12 months of treatment with IM, with the highest fold change of -10.64 at 6 months.35 Several miRNAs may have predictive value. We found that miR-221-3p was significantly upregulated at diagnosis compared to 3-month treatment. It has been also demonstrated that miR-221 is upregulated in CML blast crisis. These observation may suggest its potential role as therapeutic target. When considering different miRNA, our results are in line with previous reports showing that miR-451 can be useful as a diagnostic factor, since we observed downregulation of this molecule at diagnosis compared to 3-month.36

Analysis of functional annotations of differentially expressed miRNAs was assigned to specific biological processes according to the gene ontology (GO) classification (Figure 6). Interestingly, similar to our mRNA analysis, differentially expressed miRNAs were annotated to several immunity-related processes, including the T cell receptor signalling pathway, T cell costimulation, immune response, and inflammatory response. The list of the most up- or downregulated miRNAs detected in our study are specified in Tables 9–11. The list of the most dysregulated miRNAs with their target genes, which are involved in the immune response at diagnosis compared to 3, 6 and 12 months of treatment, is presented in Table 12.

Table 9

The List of the Most Upregulated/Downregulated miRNAs in Mononuclear Cells of Bone Marrow in the 0-Month Group Compared to the 3-Month-Group

miRNA Name	Role	Fold Change	Reference
mir-361-5p	Suppress angiogenesis via reducing the expression of VEGF-A and downregulation of expression of IL-6 and IL-8	17.73	[30]
mir-145-5p	Silencing Semaphorin 3A which acts as a terminator of T and B cells and lowering concentrations of proinflammatory cytokines	16.96	[31]
miR-3135b	Downregulated in colorectal cancer cells, play protective role for Golgi apparatus fragmentation induced by chemotherapy	12.51	[36]

Notes: The table comprises miRNA name, fold change in our study and potential role described in the literature.

Table 10

The List of the Most Upregulated/Downregulated miRNAs in Mononuclear Cells of Bone Marrow in the 0-Month Group Compared to the 6-Month-Group

miRNA Name	Role	Fold Change	Reference
mir-1246	MiR-1246 has been identifed as a tumor promoter in gall blader cancer	22.09	[37]
hsa-let-7e-5p	Potential biomarker for rectal carcinoma-associated heterochrony hepatic metastasis	11.86	[38]
miR-4505	miR‐4505 is associated with the transcription of heat shock protein family members	11.51	[39]

Notes: The table comprises miRNA name, fold change in our study and potential role described in the literature.

Table 11

The List of the Most Upregulated/Downregulated miRNAs in Mononuclear Cells of Bone Marrow in the 0-Month Group Compared to the 12-Month-Group

miRNA Name	Fold Change	Role	Reference
mir-150-5p	In 3 month: −5.88	mir-150-5p inhibits the CD28/B7 co-stimulatory pathway which result in decreasion of production of inflammatory cytokines, such as IL-2 and TNF	[40]
	In 6 month: −10.64
	In 12 month: −8.93	mir-150-5p can be considered as a biomarker in the diagnosis of AML	[41]
mir-146a-5p	In 12 month: 8.41	Target the TNF receptor-associated factor 6 (TRAF6) and IL-1 receptor-associated kinase 1 (IRAK1) genes	[42]

Notes: The table comprises miRNA name, fold change in our study and potential role described in the literature.

Table 12

The List of Most Dysregulated miRNAs with Their Target Genes, Which are Involved in Immune Response at Diagnosis Compared to 3, 6 and 12 Month of Treatment

miRNA Name	Fold Change	Direction of Target Gene Regulation	Target Gene	Fold Change Target Gene mRNA
hsa-miR-145-5p	16.96 vs 0 month	↓	CD28	−2.01
hsa-miR-199a-3p	8.17 vs 3 month	↓	PYHIN1	−3.65
hsa-miR-4505	11.51 vs 6 month	↓	CD28	−2.24
hsa-miR-145-5p	16.95 vs 3 month	↓	BACH2	−2.45
hsa-let-7e-5p	12.79 vs 3 month	↓	CYP46A1	−2.07
hsa-let-7f-5p	7.11 vs 3 month	↑	DUSP1	1.72
hsa-miR-6763-5p	5.15 vs 12 month	↓	POU2F2	−1.36

Notes: If presented miRNAs downregulate expression of their target genes, the direction of target genes regulation is marked as “↓” and the fold change of target gene is negative, indicating decreased expression. If presented miRNAs upregulate expression of their target genes, the direction of target genes regulation is marked as “↑” and the fold change of target gene is negative, indicating increased expression.

Conclusions

Our research provides new data on significant molecular changes in BM cells during TKI treatment in CML patients, using analysis of global expression of RNA as well as miRNAs and their target genes. We observed the greatest changes in gene expression at month 3 compared to diagnosis (0-month), with the XIST gene being the most affected. The extremely high expression of this gene at the time of diagnosis suggests possible implications for future therapies affecting the inhibition of expression of this gene. Cluster analysis revealed downregulation of genes involved in cell proliferation after 3 months of TKI treatment. However, it is more interesting that using this approach, we have demonstrated upregulation of genes responsible for T cell functioning eg, T cell activation, T cell costimulation and antigen processing and presentation. This phenomenon can be interpreted as the gradual renewal of immune system function. miRNAs are being considered as candidates for biomarkers of early response, which is significant for achieving good treatment outcomes. Our results indicate that major changes in the immune system in CML patient occurred in the first three months of treatment with TKIs. We found that miRNAs expression such as miR-361-5p, miR-145-5p, miR-652-3p, are markedly changed during successful TKI treatment. Other findings suggest miR-221 as prognostic factor of the success of therapy. Our results also implied diagnostic relevance of miR-451a in CML. It is possible that the application of our results in the future will be helpful in selecting patients with resistant disease or in predicting the success of the TKI treatment or even treatment discontinuation. However, further research in this interesting field is needed to provide a basis for novel approaches in treatment of this group of patients.

Study Limitations

Although our study provided valuable results, it did have some drawbacks. The first is the small number of samples in the array experiment. Moreover, the experiment was carried out on mononuclear bone marrow cells, thus it remains unknown whether the upregulation of T-cell related gene expression is due to T cell percentage increase in the sample or actual cellular response to TKI treatment. However, the presented study is basic and preliminary in nature, therefore we must be rather cautious when interpreting the obtained data. Further studies are required to fully elucidate the role of immunological pathways in the course of TKI treatment for CML.