Widespread Aberrant Alternative Splicing despite Molecular Remission in Chronic Myeloid Leukaemia Patients

Abstract

Vast transcriptomics and epigenomics changes are characteristic of human cancers, including leukaemia. At remission, we assume that these changes normalise so that omics-profiles resemble those of healthy individuals. However, an in-depth transcriptomic and epigenomic analysis of cancer remission has not been undertaken. A striking exemplar of targeted remission induction occurs in chronic myeloid leukaemia (CML) following tyrosine kinase inhibitor (TKI) therapy. Using RNA sequencing and whole-genome bisulfite sequencing, we profiled samples from chronic-phase CML patients at diagnosis and remission and compared these to healthy donors. Remarkably, our analyses revealed that abnormal splicing distinguishes remission samples from normal controls. This phenomenon is independent of the TKI drug used and in striking contrast to the normalisation of gene expression and DNA methylation patterns. Most remarkable are the high intron retention (IR) levels that even exceed those observed in the diagnosis samples. Increased IR affects cell cycle regulators at diagnosis and splicing regulators at remission. We show that aberrant splicing in CML is associated with reduced expression of specific splicing factors, histone modifications and reduced DNA methylation. Our results provide novel insights into the changing transcriptomic and epigenomic landscapes of CML patients during remission. The conceptually unanticipated observation of widespread aberrant alternative splicing after remission induction warrants further exploration. These results have broad implications for studying CML relapse and treating minimal residual disease.

Keywords: BCR-ABL1; CML; DNA methylation; alternative splicing; cancer; epigenetics; histone modifications; intron retention; transcriptomic complexity.

Figure 1

Molecular remission is reflected in transcriptomic profiles. (A) The Top 200 differentially expressed genes with lowest adjusted p-value in likelihood ratio test. (B) Predicted cell type enrichment in peripheral blood samples (diagnosis/red: D1–27; remission/green: R1–10; control/blue: C1–6). The stacked bar graph shows normalized cell type enrichment scores for patient and control transcriptomes based on bulk RNA sequencing data. Predictions were made with the xCell R package (github.com/dviraran/xCell) [16]. HSC—hematopoietic stem cell, CLP—common lymphoid progenitor, CMP—common myeloid progenitor, GMP—granulocyte/macrophage progenitor, MEP—megakaryocyte/erythroid progenitor, MK—megakaryocyte, NK—natural killer, DC—dendritic cell. (C) Principal component analysis separates diagnosis samples from the common features of remission samples and healthy controls.

Figure 2

DNA methylation in chronic myeloid leukaemia (CML). (A) Frequency of differentially methylated regions (DMRs) based on delta values between pooled diagnosis and remission samples (left) as well as remission and control samples (right). (B) Circos plot illustrating genome-wide differences in the frequency of DMRs (central ring with blue scatter plot: diagnosis vs. remission, inner ring with red scatterplot: remission vs. control). This plot was generated using the Circos software [17]. (C) Scatterplot illustrating the number of significantly differentially methylated CpG sites per DMR and associated delta values (red points, diagnosis vs. remission; blue points, remission vs. control). (D) Enriched Gene Ontology (GO) terms associated with DMRs (remission vs. diagnosis) in gene promoter regions (red dashed line marks p = 0.05). (E) Distribution and density of DMRs across gene features. RC—remission vs. control, DR—diagnosis vs. remission, DC—diagnosis vs. control.

Figure 3

Global alternative splicing patterns are similar despite major molecular remission. (A) Scatterplots illustrating differential IR events between pooled diagnosis/remission/control samples binned into hexagons to avoid overplotting (yellow: adj.-p ≤ 0.05; red: adj.-p ≤ 0.01). (B) Venn diagram showing intersections of mutually differentially retained introns. (C) Dynamics of IR across CML phases. CML patient data from Branford et al [11]. was analysed for the occurrence of IR events. Left: Sample sizes for Control, CP: Chronic Phase, AP: Accelerated Phase, MBC: Myeloid Blast Crisis and LBC: Lymphoid Blast Crisis samples. Middle: Sample numbers matched from different phases. Right: Boxplot of IR frequencies in samples from different CML phases. Mann–Whitney U test; ** (p < 0.01), ns (p > 0.05). (D) Differential exon usage in pooled diagnosis vs. control (left), diagnosis vs. remission (middle) and remission vs. control samples. Numbers following the red arrows indicate how many upregulated DUEs in the first of the two conditions compared. Numbers following the green arrows indicate how many downregulated DUEs in the first of the two conditions compared. (E) Intersection of differentially used exons in diagnosis and remission samples compared to controls (the top-enriched GO term and associated genes are listed in the middle).

Figure 4

Differentially retained introns inhibit the expression of cell cycle regulators following major molecular remission. (A) Scatterplots illustrating the relationship between differentially retained introns (∆IR) and differential gene expression (log2 expression fold change). Introns (|∆IR| ≥ 0.1) with inverse relationship to host gene expression (|log(FC)| ≥ 1) are highlighted in red. (B) Coverage plots of selected differential IR events (left) and qRT-PCR validation of differential IR between diagnosis and control (right). qRT-PCR validation was performed in triplicates on independent diagnosis samples and results were normalised to B2M expression. D—diagnosis, C—control. (C) Blood cell type expression of genes with inverse IR/gene expression relationships (diagnosis vs. control; blue—myeloid, magenta—lymphoid cells). The data were retrieved from the Human Protein Atlas (v19.3; http://www.proteinatlas.org) [24].

Figure 5

Regulation of IR in CML. (A) Two genes that most strongly correlate/anti-correlate with IR frequencies (Spearman correlation). (B) CpG methylation around intron splice sites (± 200 bp) and in the centre (200 bp) of retained- (red) and non-retained introns (blue) in CML diagnosis samples. General linear hypothesis test; p < 2.2 × 10⁻¹⁶ (p-adj. < × 1 × 10⁻¹⁰). (C) The effect of RNA binding proteins on IR frequency. RNA-seq data of small hairpin RNA (shRNA)-mediated knockdown of indicated RNA binding proteins in K562 cells was retrieved from the ENCODE project (encodeproject.org). WT—wild type. (D) Expression of RNA binding proteins with enriched binding motifs near frequently retained introns [27]. TPM - transcripts per million mapped reads; C—control; D—diagnosis; R—remission; NS.—non-significant, ** p < 0.01, *** p < 0.001. (E) Histone modifications with increased numbers of peaks near splice sites of retained introns (red: retained; blue: non-retained introns). The analysis was repeated for each expression quartile to exclude the possibility that increased histone modifications are associated with gene expression. (F) Possible mechanisms of IR regulation in CML discussed in the text.