Widespread aberrant alternative splicing despite molecular remission in chronic myeloid leukemia patients

Patient samples

We retrieved nine diagnostic specimens (total leukocytes from peripheral blood collected in EDTA) from patients enrolled in the RESIST (1x imatinib; ACTRN12610000055000), ENESTxtnd (4x nilotinib, NCT01254188) and PINNACLE (4x nilotinib + pegylated interferon, ACTRN12612000851864) trials, with matched remission specimens (major or complete molecular response at 12 or 24 months). In addition, eight CML patient samples (7 diagnostic and 1 matched remission) from the French Persistem study were obtained. For the validation of IR events and fusion transcripts we retrieved further matched diagnosis/remission samples from the TIDELII trial (20x imatinib, ACTRN12607000325404). RNA was isolated from mononuclear cells (from peripheral blood collected in LiHep) which had been cryopreserved and subsequently thawed into TRIzol. All studies were approved by the institutional review boards (HREC protocol No: 131015, 070718c, 081211, and 101010 - Royal Adelaide Hospital) and are in accordance with the Declaration of Helsinki. Patients provided written informed consent.

RNA isolation and mRNA-seq

Total RNA was isolated from diagnosis, remission, and control samples (peripheral blood mononuclear cells - PBMCs) using Trizol according to the manufacturer’s protocol. For mRNA sequencing, poly-A enriched mRNA libraries were prepared from 1 μg of total RNA using the TruSeq Stranded RNA sample prep kit (Illumina), prior to paired-end sequencing using the HiSeq 2500 platform. RNA sequencing was performed in triplicates for each sample. The raw sequencing data (fastq files) have been deposited at GEO under accession GSE144119.

mRNA sequencing data analysis

Paired-end RNA-sequencing reads (125 nt) were trimmed and mapped to the human reference genome hg38 using STAR¹⁴. Quality control of raw and mapped sequencing reads was performed using FASTQC (github.com/s-andrews/FastQC), RSeQC¹⁵, and multiqc¹⁶. Putative confounders and batch effects were excluded using principal components analyses (Supplementary Figure 1). STAR-FUSION¹⁷ was used for the identification of fusion genes and Fusion Inspector (FusionInspector.github.io) for in silico validation of the predicted gene fusions. Gene expression levels specified as transcripts per million (TPM) were determined using Salmon 0.14.1¹⁸.

General statistics on alternative splicing events were determined using the rMATS software¹⁹. We used our IRFinder algorithm for the detection of IR events in introns extracted from an Ensembl gtf file (regions between two adjacent exons)²⁰. IRFinder estimates the abundance of IR by computing the ratio between gene transcripts retaining an intron and the sum of all transcripts of the respective gene (more information is provided in the Supplementary Materials).

Differentially used exons were determined using the R Bioconductor package DEXSeq²¹. Differential exon usage (DEU) = changes in the relative usage of exons:

IR validation

Extracted RNA (2 µg) was treated with DNAse Turbo (Invitrogen) and reverse transcribed into cDNA using oligo(dT) priming and Superscript III (ThermoFisher) according to the protocol supplied by the manufacturer. For each sample, a corresponding control without reverse transcriptase was used to assess DNA contamination. For qPCR, cDNA templates were amplified, and the Ct values were quantified with SYBR Green Master Mix (ThermoFisher). Beta-2 microglobulin (B2M) was used as the normalization control for cDNA input. Experiments were performed using the CFX96 Real-Time PCR System (BIORAD). The list of primers is provided in Supplementary Table 1.

Whole genome bisulfite sequencing (WGBS)

WGBS libraries were prepared following Illumina’s “Whole-Genome Bisulfite sequencing for Methylation Analysis” protocol. Briefly, 1 μg of genomic DNA was spiked with 0.5% unmethylated lambda DNA and sonicated to generate fragments of size between 150 to 300 bp. Library preparation was performed using Illumina’s Paired-end DNA Sample Prep Kit (discontinued, Illumina, CA, USA) according to the manufacturer’s protocol. The size-selected libraries were then subjected to bisulfite conversion as previously described²². Adaptor-ligated bisulfite treated DNA was enriched by 10 cycles of PCR amplification using the PfuTurbo Cx Hotstart DNA Polymerase (Stratagene). Qualitative and quantitative checks of the libraries were performed using Agilent’s High sensitivity DNA kit (Agilent) and KAPA Library quantification kit (KAPA Biosystems). Three lanes of paired-end 100 bp sequencing was performed for each of the library on the Illumina HiSeq2500 platform using the TruSeq v3 cluster kits and SBS kits to achieve coverage ranging between 25-30x.

WGBS data analysis

Reads were processed and aligned to the human (hg38) reference genome using Meth10X (github.com/luuloi/Meth10X)²³. In short, the Meth10X pipeline takes raw reads in fastq format and trims the adaptors, which are then aligned to the human reference genome using bwa-meth (github.com/brentp/bwa-meth). The generated bam files were marked with duplication and merged if necessary. Estimation of the duplication rate, coverage bias (genomic features) and methylation bias in reads was carried out to provide quality control. All the metrics such as percentage of unmapped/mapped read metrics, mapping quality distribution, GC content distribution, insert size distribution and coverage distribution were generated by Qualimap 2 for further evaluation of the alignment. Finally, a count table and bigwig files of methylated and coverage at each CpG site in the genome was constructed for all samples. Differentially methylated regions (DMRs) were identified using MethPipe²⁴ with count table of all samples as input. The analysis of differential methylation and IR is described in the Supplementary Materials.

ChIP-seq data analysis

ChIP-seq data of histone modifications in K562 cells was retrieved from ENCODE (encodeproject.org) and further processed as described in the Supplementary Materials.

Statistical analyses

All statistical analyses were performed in R v.3.6.2. The Wald test implemented in the R package DESeq2²⁵ was used for the identification of statistically significant differential alternative splicing events, and differential gene expression analysis. P-values were adjusted for the false discovery rate using the Benjamini-Hochberg procedure. Functional enrichment analysis of differentially expressed genes and genes with differentially retained introns was performed using DAVID 6.7²⁶ with the human genome and the set of expressed genes, respectively, used as a background.

Heterogeneous CML transcriptomes converge at remission

RNA sequencing reads mapping to the BCR-ABL1 locus were found in diagnosis samples only (Supplementary Figure 2A), confirming the results of clinical qRT-PCR-based BCR-ABL1 quantification (Supplementary Data 1). Apart from the characteristic BCR-ABL1 fusion we identified other recurring intra- and inter-chromosomal fusions (Supplementary Figure 2A&B; Supplementary Data 2) that occur in diagnosis, remission or control samples. A novel fusion transcript between the myeloid cell-specific transmembrane glycoprotein CLEC12A (12p13.31) and the microRNA miR-223 host gene MIR223HG (Xq12) was expressed at higher levels (FC = 7.16) at diagnosis compared to remission (Supplementary Figure 2C&D). For most patients we observed a global reduction of gene expression in remission compared to diagnosis (Supplementary Figure 3A). This can be attributed to the reconfigured blood cell composition, which resembles lymphoid cell-enriched normal controls (Figure 1A&B). A significant downregulation was observed in genes encoding kinases other than BCR-ABL1 and kinase-like proteins such as CENPE, CDK1, AURKB, MELK, and BUB1B (Supplementary Figure 4). Reduced expression at remission, similar to levels observed in control samples, also affected cyclins such as CCNA1, CCNB2 and other cell cycle regulators as well as genes related to DNA replication and repair (Figure 1A & Supplementary Figure 5A). In this context, no consistent differences were noted between patients treated with imatinib or nilotinib (Supplementary Figure 1).

• Download figure
• Open in new tab

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. 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)⁴⁷. HSC – hematopoietic stem cell, CLP – common lymphoid progenitor, CMP – common myeloid progenitor, GMP - granulocyte macrophage progenitor, MEP - megakaryocytic erythroid progenitor, NK – natural killer, DC – dendritic cell. (C) Principal component analysis separates diagnosis samples from the common features of remission samples and healthy controls.

A set of genes (n=81), that are functionally associated with oxygen and bicarbonate transport (GO:0015671; GO:0015701) as well as the hemoglobin’s chaperone pathway (Biocarta), was significantly differentially expressed in all three comparisons (diagnosis vs control, remission vs control, and diagnosis vs remission; Supplementary Figure 5B). Overall, we observed that patients diagnosed with CML have highly heterogeneous transcriptomes. In contrast, remission transcriptomes are more alike (Figure 1C) and exhibit fewer differences in comparison to healthy control samples (Supplementary Figure 3B).

DNA methylation profiles return to normal at remission

To characterise putative epigenetic causes for differential gene expression in CML patients prior to and after successful TKI treatment, all matched diagnosis/remission samples as well as all control samples were subjected to WGBS. The analysis revealed almost 25,000 differentially methylated regions (DMRs) in the genomes of diagnosis and remission samples (Figure 2A, left). Vast DNA methylome differences were found in every matched patient sample (Supplementary Figure 6) and across all chromosomes (Figure 2B). In contrast, only 710 DMRs were observed between remission and control samples (Figure 2A, right). Indeed, multidimensional scaling showed that DNA methylation profiles of remission samples are closer to those from healthy donors rather than to matched diagnosis samples (Supplementary Figure 7). DMRs between diagnosis and remission samples contain, on average, a larger number of differentially methylated CpG sites (Figure 2C). This is consistent with previous reports that suggest an increase in CpG site methylation during CML progression that affects tumour-suppressor genes and regulators of cell proliferation¹⁰. DMRs are associated with processes of the innate immune response, apoptosis, cell differentiation, and cell migration (Figure 2D).

• Download figure
• Open in new tab

Figure 2. DNA methylation in 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⁴⁸. (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 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.

Our results show that DNA methylation returns to normal levels after successful TKI treatment. Most DMRs can be found in intronic regions, however, the highest density of DMRs was observed in 5′ untranslated regions (Figure 2E), suggesting a key role for differential DNA methylation in orchestrating CML gene expression regulation.

Aberrant alternative splicing distinguishes remission samples from healthy controls

To assess whether major or complete molecular remission also affects gene isoform expression, we conducted a systematic analysis of the five major forms of alternative splicing (Supplementary Figure 8A). In a pairwise comparison of alternative splicing event frequencies between matched diagnosis/remission samples, no consistent trend towards an increase or decrease in either of the two conditions was found (Supplementary Figure 8B). This suggests, that in contrast to gene expression profiles, splicing patterns remain atypical after major or complete molecular remission.

Alternative splicing analysis of the Cancer Genome Atlas (TCGA) cohort by Dvinge et al. (2015) revealed that the frequency of IR events is consistently higher in acute myeloid leukemia compared to normal controls⁸. Although mean IR frequencies are also higher in CML samples (diagnosis, = 3 382; remission, = 3 484) compared to healthy donors ( = 2 822), high inter-donor variability led to differences in IR frequencies that did not reach significance (Supplementary Figure 8C). However, when we subsampled the RNA sequencing data to facilitate a uniform sequencing depth for comparison, we observed a clear separation with significantly increased IR frequencies at diagnosis and remission compared to controls (Supplementary Figure 9). Moreover, we found that the IR ratios are consistently higher in diagnosis and remission samples compared to healthy controls (Figure 3A). Differences between diagnosis and remission samples were less pronounced (Figure 3A; Supplementary Figure 10). While IR profiles are unsuitable for patient stratification (Supplementary Figure 11), differentially retained introns between diagnosis and remission samples are, similarly to differentially expressed genes, associated with cell cycle regulators as well as genes related to DNA replication and repair (Supplementary Figure 12).

• Download figure
• Open in new tab

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: p_adj <= 0.05; red: p_adj <= 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.⁹ 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. (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).

However, in remission IR primarily affects genes involved in splicing processes (Supplementary Figure 12) suggesting the aberrant splicing persists in remission due to auto-regulatory processes. Overall, IR profiles are highly dynamic in CML, which becomes apparent when considering the large overlap of differentially retained introns (n = 419) in all three comparisons (Figure 3B). We confirmed these results with an independent dataset of 59 CML patients from Branford et al.⁹, which includes samples from various stages of disease progression, and from 4 healthy controls (Figure 3C). The number of IR events in this cohort was highly variable. While there are on average more IR events in chronic phase compared to healthy controls, the number drops significantly at blast crisis (Figure 3C, right).

Examining differential exon usage in our samples we observed a similar trend although not as marked as observed in differential IR, i.e. both diagnosis and remission samples exhibit a larger number of differentially used exons (DUEs) compared to normal controls. Most of these DUEs are upregulated (diagnosis, 78%; remission, 81%; Figure 3D). A large fraction of the DUEs at diagnosis (72%) remain differentially expressed at remission as well. Some genes of these mutual DUEs are associated with DNA repair mechanisms (top-enriched GO term; Figure 3E). Only 26 exons in total are differentially used in diagnosis vs remission samples (Figure 3D).

Modulation of IR levels supports lineage-specific gene expression in CML

IR is accepted as a mechanism of post-transcriptional gene regulation triggering nonsense-mediated decay due to the presence of premature termination codons within the retained introns^27,28. We analysed the relationship between IR and gene expression changes but overall the two variables did not significantly correlate (Figure 4A).

• Download figure
• Open in new tab

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 v control; blue - myeloid, magenta – lymphoid cells). The data was retrieved from the Human Protein Atlas (v19.3; http://www.proteinatlas.org)⁴⁹.

However, many of the genes, whose expression does anti-correlate with the retention of their introns, are involved in cell cycle regulation (e.g. CCL3, CKS2, and SPAG5), and DNA replication/repair (e.g. FANCI/FANCD2, and PCNA) (Figure 4A).

We selected three differentially retained introns of cell cycle regulator genes for experimental validation (CKS2 intron 1, CCL3 intron 1, and SERPINB1 intron 2). Based on the RNA sequencing data analysis, all three introns have an increased expression and higher IR ratio at diagnosis compared to healthy controls (Figure 4A). We confirmed these observations via qRT-PCR (Figure 4B). The IR ratios of CKS2 intron 1, and SERPINB1 intron 2 increase further at remission resulting in a reduction of host gene expression, while the IR ratio of CCL3 intron 1 is reduced at remission and associated with an increase in gene expression (Figure 4A & B).

Interestingly, some instances of inverse relationships between IR and gene expression affect lineage-specific genes (Figure 4C). For example, the chemokine Interleukin-8 (CXCL8), which is predominantly expressed in neutrophils, is upregulated at diagnosis. This could be attributed to the increased myeloid progenitor cell abundance (Figure 1B), but also to the reduced IR levels in CXCL8 mRNA (Figure 4A). In control samples, the IR levels increase from 12.9% to 26.7% in CXCL8 (intron 2), while gene expression is reduced 4-fold. CXCL8 promotes CML cell proliferation²⁹ and previous reports suggest that CXCL8 expression is modulated both by BCR-ABL1 expression (causing increase) and TKI treatment (causing inhibition of CXCL8)³⁰. In contrast, PRKCH is downregulated at CML diagnosis along with increased IR levels (Figure 4A). PRKCH is a serine-threonine kinase that regulates hematopoietic stem cell function and is specifically expressed in lymphoid cells. Notably, high PRKCH expression has been associated with poor prognosis in acute myeloid leukemia³¹.

These observations suggest that lineage-specific expression of drivers of cell proliferation and hematopoietic cell differentiation can be mediated through IR level changes in CML.

Epigenetic features modulate intron retention in CML

Characteristics of the retained introns in our CML cohort are in accordance with previously described attributes, such as their relatively short length (Supplementary Figure 13), and higher GC content (Supplementary Figure 14). While these are intrinsic features of IR, we sought potential trans-regulators and found that the gene that most strongly correlates (ρ = 0.64) with IR frequencies is MED21 (Mediator Complex Subunit 21; Figure 5A). MED21 regulates gene transcription by interacting with RNA polymerase II³². RNA polymerase II elongation rates influence splice site recognition and IR¹³. Interestingly, the expression of ZNF160, encoding a Zinc Finger Protein and repressor of transcription, strongly anti-correlates (ρ = −0.62) with IR frequencies (Figure 5A).

• Download figure
• Open in new tab

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-value < 2.2e-16 (adj.-p < 1e-10). (C) The effect of RNA binding proteins on IR frequency. RNA-seq data of 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²⁰. TPM - transcripts per million mapped reads; C - control; D - diagnosis; R – remission; n.s. – non-significant, * p < 0.05, ** 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.

We have previously shown that DNA methylation regulates IR in myeloid cells¹³. To confirm this mode of regulation in our CML cohort, we compared CpG methylation around the splice sites (±200 bp) and in the centre (200 bp) of retained and non-retained introns. We confirmed decreased DNA methylation levels near the 3□ splice junction and within the body of retained introns (Figure 5B). Similar patterns were observed in remission and control samples suggesting that this mode of IR regulation is also present in lymphoid cells (Supplementary Figure 15).

To further investigate the role of epigenetics in the regulation of IR in CML, we retrieved RNA-seq and ChIP-seq data of the K562 CML cell line from ENCODE. We analysed transcriptomics data of RNA binding protein knockdown experiments and observed significantly increased IR occurrences with the most drastic IR increase (∼5 fold) following Poly(RC) Binding Protein 2 (PCBP2) knockdown (Figure 5C). All genes with enriched binding motifs near frequently retained introns²⁰ were downregulated at CML diagnosis, except for SRSF5 (Supplementary Figure 16). Three of these RNA binding proteins, including PCBP2, remained downregulated in remission (Figure 5D).

To identify additional potential explanations for the dysregulation of IR in CML, we analysed histone modifications in K562 cells. While higher levels of H3K36 trimethylation have previously been associated with IR^33,34, we found only marginal differences in H3K36me3 levels in retained versus non-retained introns. Increasing numbers of H3K36me3 peaks around introns were observed in more highly expressed genes (Supplementary Figure 17). However, we found that monomethylation and acetylation of histone H3 lysines 4 and 9, respectively (H3K4me1 and H3K9ac) were significantly enriched near splice sites of retained introns and within their intron bodies in contrast to non-retained introns (Figure 5E). While a significant enrichment of H3K9ac has been found to co-occur at many gene regulatory elements in mouse embryonic stem cells³⁵, we here established a strong association of H3K9ac with IR. Elevated levels of H3K4 methylation in association with IR have been previously reported by Zhou et al.³³.

In summary, our analysis of CML patient transcriptomes revealed aberrant alternative splicing both at diagnosis and remission. Elevated IR ratios in CML can be explained by multiple regulatory mechanisms including the regulation of RNA Pol II elongation, demethylation of CpG sites, and histone modifications (Figure 5F). The TKI treatment itself seems to have no effect on IR (Supplementary Figure 18).