Aberrant splicing in B-cell acute lymphoblastic leukemia.

Aberrant splicing is a hallmark of leukemias with mutations in splicing factor (SF)-encoding genes. Here we investigated its prevalence in pediatric B-cell acute lymphoblastic leukemias (B-ALL), where SFs are not mutated. By comparing these samples to normal pro-B cells, we found thousands of aberrant local splice variations (LSVs) per sample, with 279 LSVs in 241 genes present in every comparison. These genes were enriched in RNA processing pathways and encoded ∼100 SFs, e.g. hnRNPA1. HNRNPA1 3'UTR was most pervasively mis-spliced, yielding the transcript subject to nonsense-mediated decay. To mimic this event, we knocked it down in B-lymphoblastoid cells and identified 213 hnRNPA1-regulated exon usage events comprising the hnRNPA1 splicing signature in pediatric leukemia. Some of its elements were LSVs in DICER1 and NT5C2, known cancer drivers. We searched for LSVs in other leukemia and lymphoma drivers and discovered 81 LSVs in 41 additional genes. Seventy-seven LSVs out of 81 were confirmed using two large independent B-ALL RNA-seq datasets, and the twenty most common B-ALL drivers, including NT5C2, showed higher prevalence of aberrant splicing than of somatic mutations. Thus, post-transcriptional deregulation of SF can drive widespread changes in B-ALL splicing and likely contributes to disease pathogenesis.

INTRODUCTION

Despite advances in the treatment of pediatric B-ALL, children with relapsed or refractory disease account for a substantial number of childhood cancer-related deaths. Adults with B-ALL experience even higher relapse rates and long-term event-free survival of <50% (1). Recently, significant gains in the treatment of B-ALL have been achieved through the use of immunotherapies directed against CD19, a protein expressed on the surface of most B-cell neoplasms (2,3). These gains culminated in the recent FDA approval of tisagenlecleucel and axicabtagene ciloleucel, CD19-redirected chimeric antigen receptor (CAR) T-cell immunotherapies, for patients with refractory/relapsed B-cell malignancies. However, relapses occur in 10–20% of patients with B-ALL treated with CD19-directed immunotherapies, often due to epitope loss and/or B-cell de-differentiation into other lineages (4–7). Other targets for immunotherapy include CD20 and CD22 (8–11). However, neither antigen is uniformly expressed in B-ALL, and factors accounting for this mosaicism are poorly understood (3).

We previously reported a new mechanism of pediatric B-ALL resistance to CD19-directed immunotherapy. We discovered that in some cases, resistance to CD19 CAR T cells was generated through alternative splicing of CD19 transcripts. This post-transcriptional event was mediated by a specific splicing factor SRSF3 and generated a CD19 protein isoform invisible to the immunotherapeutic agent via skipping of exon 2 [(12,13), reviewed in (14)].

Our discovery of a resistance mechanism based on alternative splicing prompted us to investigate the extent of this phenomenon in additional B-ALL cases. While driver mutations in splicing factors such as SRSF2, SF3B1 and U2AF1 have recently been discovered in myelodysplastic syndrome/acute myelogenous leukemia (15–17) and chronic lymphocytic leukemia (18,19), SF mutations have not been reported in B-ALL. Nevertheless, our prior work suggested the possibility that SRSF3 (and by inference other SFs) could be deregulated in B-ALL (12), bringing about wide-spread splicing aberrations.

This model would be particularly attractive because B-ALL is a chromosome translocation-driven disease where the prevalence of somatic mutations and copy number variations is relatively low. For example, the commonly mutated IKZF1 gene (which encodes the Ikaros transcription factor) is affected by missense mutations in just ∼20% of B-ALL cases. Similarly, mutations in the key tumor suppressor gene (TSG) TP53 are found in only ∼7% of B-ALLs (per COSMIC database) (20,21). In addition, both genes are robustly transcribed across individual B-ALLs and thus are not epigenetically silenced. This raises the possibility that they and other TSGs are dysregulated by post-transcriptional events, such as alternative splicing.

MATERIALS AND METHODS

Bone marrow fractionation

Isolated mononuclear cells and whole bone marrow aspirates were obtained, respectively, from the University of Pennsylvania Stem Cell and Xenograft Core facility and CHOP Hematopathology Laboratory. For pediatric bone marrow samples, mononuclear cells were isolated by spinning over Ficoll gradient, as described earlier (22). Residual red blood cells were lysed with ammonium chloride lysis buffer with gentle rocking at room temperature for 10 min. Cells were pelleted by spinning at 250 × g for 10 min at 4°C and washed once with cold PBS/2% FBS. Cells were resuspended in 1 ml PBS/2%FBS and incubated with 500 μl FC Block on ice for 10 min. Cells were stained with 1 ml CD34-PE, 500 μl CD19-APC, and 500 μl IgM-FITC for 30 min on ice. Cells were pelleted at 1300 RPM for 6 min at 4°C and washed twice in cold PBS. Cells were resuspended in 3 ml PBS/2%FBS containing 1 µl/ml of 0.1 mg/ml DAPI. Cells were sorted 4-ways using MoFlo ASTRIOS directly into RLT Lysis buffer (Qiagen) at a ratio of 1:3.

Primary B-ALL samples acquisition

Twenty four primary pediatric B-ALL samples were obtained from the CHOP Center for Childhood Cancer Research leukemia biorepository. Mononuclear cells from fresh bone marrow or peripheral blood specimens were purified via Ficoll gradient separation (22) and cryopreserved for downstream experimental use.

RNA-seq of bone marrow fractions, primary samples and cell lines

RNA was isolated using Qiagen RNeasy Mini Kit. RNA integrity and concentration were determined using Eukaryote Total RNA Nano assay on BioAnalyzer (CHOP NAPCore). RNA-seq was performed on 10 ng–1 μg of total RNA according to the GeneWiz Illumina Hi-seq protocol for poly(A)-selected samples (2 × 150 bp pair-end sequencing, 350M raw reads per lane).

RNA-Seq alignment, quantification and differential expression

Fastq files of RNA-seq obtained from GeneWiz were mapped using STAR aligner (23). STAR was run with the option ‘alignSJoverhangMin 8′. We generated STAR genome reference based on the hg19 build. Alignments were then quantified for each mRNA transcript using HTSeq with the Ensembl-based GFF file and with the option ‘-m intersection-strict’. Normalization of the raw reads was performed using the trimmed mean of M-values (TMM). Principal component analysis (PCA) was done on normalized count values of the samples using a correlation matrix and a calculated score for each principal component. Differential expression of wild-type and knock-down or pro-B and B-ALL RNA-Seq datasets were assessed based on a model using the negative binomial distribution, a method employed by the R package DESeq2 (24). Subsequent bar charts were generated using the R package ggplot2. Those differential genes that had a P-value of <0.05 were deemed as significantly up or down-regulated.

Splicing analysis

In order to detect LSVs, the MAJIQ tool (version 1.03) was used (25). MAJIQ was run on the Ensembl-based GFF annotations, disallowing de novo calls (26). LSVs that had at least a 20% change at a 95% confidence interval between two given conditions were chosen for further analysis. The 95% confidence interval in this context is defined as follows: given observed reads in both experiments there is a 95% posterior probability, according to the MAJIQ model, that there is a change of 20% or more in PSI. In-house customized Perl and Ruby scripts were used to filter for events that corresponded to exon inclusion events only, forcing one event per LSV. For the genes containing differential LSVs in all 18 comparisons, enriched gene ontologies (GO) were identified using the gene functional classification tool DAVID (v. 6.7), reporting the top most significant hits based on P-value and false discovery rate. Heatmaps were generated for each ΔPSI value of each LSV that passed the 20% change threshold between Pro-B and B-ALL and additionally filtered for a frequency of n ≥ 2 in B-ALL samples. The R package gplots was used for this purpose. LSVs detected by MAJIQ were also validated using Leafcutter version 0.2.7 (27). In addition, splicing of HNRNPA1 3′ UTR was validated by RT-qPCR using a forward primer spanning exons 9–10 and a reverse primer in exon 11. Expression levels were normalized to actin and the PHL156 leukemia.

Spearman's correlations

Correlations and their significance were computed using the nonparametric Spearman's rank-order correlation implemented in R function cor.test(). Correlation analysis was performed between a) RQ values of HNRNPA1 exon 11 inclusion as determined by RT-qPCR and PSI values of HNRNPA1 exon 11 inclusion as determined by MAJIQ; b) RQ values of HNRNPA1 constitutive exons 2–3 and RQ values of HNRNPA1 exon 11 inclusion as determined by RT-qPCR (normalized to actin and PHL156), and c) normalized counts of HNRNPA1 constitutive exon 9 and cassette exon 11, as determined by RNA-seq.

Motif analysis

The motif detection pipeline utilized RBPMap (v. 1.1) with a medium stringency level (P-value ≤ 0.005) and conservation filter (28). The analysis was performed on affected exons and ±200nts into adjacent intronic regions. hnRNPA1 motifs (29) were identified and average z-scores and P-values for all significant hits were calculated.

Actinomycin D mRNA stability assay

2 × 106 Nalm6 B-ALL cells were treated with 5 μg/ml of actinomycin D for various time intervals between 0 and 24 h followed by RNA extraction. RT-qPCR was performed to assess the abundance of HNRNPA1 3′UTR isoforms using forward primer spanning exons 9–10 and reverse primers in exon 11 or spanning exons 12–13. Data was normalized to actin and 0 hr time point. Linear regression was used to plot the decay curves.

Nonsense-mediated decay inhibition via cyclohexamide treatment and UPF1/2 knockdown

2 × 106 Nalm6 B-ALL cells were treated in duplicates with either DMSO or 30 μg/ml cyclohexamide for 6 h followed by RNA extraction. Alternatively, they were electroporated with 133 nM non-targeting control pool of four siRNAs or 750nM ON-TARGET Plus Human UPF1 SMARTpools of four siRNAs + 300–750 nM ON-TARGET Plus Human UPF2 SMARTpool of four siRNAs (all from Dharmacon). To confirm inhibition of NMD, gel-based PCR with primers spanning SRSF3 poison exon 4 to detect accumulation of a canonical NMD target was performed. To determine if either HNRNPA1 3′UTR was an NMD substrate, RT-qPCR with forward primer spanning exons 9–10 and reverse primers in exon 11 or spanning exons 12–13 were performed. Data were normalized to HNRNPA1 constitutive exons 2–3 and DMSO.

siRNA knock-down of hnRNPA1 and SRSF3

Biological replicate experiments were performed on 2 × 106 P493-6 B-lymphoblastoid cells (30) electroporated using Amaxa Program O-006 with either 133 nM non-targeting control pool of four siRNAs, 300nM ON-TARGET Plus Human HNRNPA1 SMARTpool of four siRNAs or 50–300 nM ON-TARGET Plus Human SRSF3 SMARTpool of four siRNAs (all from Dharmacon). Cells were re-plated in 2 ml warm tetracycline-free RPMI for 24 h. Knockdown was validated through RT-qPCR and western blotting with anti-hnRNPA1 antibody (Abcam#ab5832) or as described previously for SRSF3 (12).

LNA-GapmeR-mediated knock-down of HNRNPA1

To knock-down hnRNPA1 by a non-RNAi mechanism, 2 × 106 P493-6 B-lymphoblastoid cells were electroporated with either 100 nM custom non-targeting control pool of four GapmeRs or 100 nM of custom HNRNPA1-ATG-targeting GapmeR using Neon Transfection System (ThermoFisher). All GapmeRs were hybrids of phosphorothioate-bonded DNA and 2′-O-methyl modified RNA. Cells were plated in 2 ml warm tetracycline-free RPMI for 48 h. RNA isolation was performed as described above. Knockdown was validated through RT-qPCR and western blotting with anti-hnRNPA1 antibody (Abcam #ab5832). Alternative splicing of DICER1 (LSV ID ENSG00000100697:95623567-95623607:target) was validated using RT-qPCR with primers amplifying the alternative exon 3–4 junction normalized to primers amplifying the constitutive exon 25. Alternative splicing of NT5C2 (LSV ID ENSG00000076685:104865463-104865558:target) was validated using RT-qPCR with primers amplifying the exon4-exon4a junction and normalized to primers amplifying constitutive exons 5–7.

Radioactive gel-based RT-PCR

To validate the regulation of alternative splicing of TP53 (LSV ID: ENSG00000141510:7573927-7574033:source) 26 cycles of radioactive RT-PCR were performed as described previously (12) using a forward primer in TP53 exon 9 and reverse primer in exon 10.

Oligonucleotide sequences

Sequences of all oligonucleotides are shown in Supplementary Table S4.

Datasets

Cancer gene symbols and annotations were downloaded from the COSMIC database (20,21). Known splice factors were annotated and obtained from published studies or Ensembl-annotated databases. Pediatric B-ALL samples from the TARGET consortium (phs000218.v19.p7) were accessed via NCBI dbGaP (the Database of Genotypes and Phenotypes) Authorized Access system (project #10088). Pediatric B-ALL samples from the St Jude Children's Research Hospital (EGAD00001002704 and EGAD00001002692) were accessed by permission from the Computational Biology Committee though The European Bioinformatics Institute (EMBL-EBI). CHOP B-ALL samples used in this submission were passed to investigators with a coded identifier and no protected health information (PHI). Basic demographic, treatment, relapse, and survival outcome data for delivered specimens were provided through an online honest broker system. Specimens were obtained via informed consent on institutional research protocols in accordance with the Declaration of Helsinki.

RESULTS

RNA-seq analysis of bone marrow-derived human B-cells

To determine if patterns of splicing dysregulation occur in B-ALL, we first generated normal B-lymphocyte datasets corresponding to potential cells of origin. To this end, we obtained from the University of Pennsylvania Stem Cell and Xenograft Core facility two healthy adult bone marrow biopsies and from the Children's Hospital of Philadelphia (CHOP) Hematopathology Laboratory two bone marrow biopsies from children without leukemia. We enriched for mononuclear cells using Ficoll gradient separation, stained cells for combinations of stage-specific surface markers, and sorted B-cell subsets using flow cytometry. Specifically, we fractionated bone marrow progenitors into early progenitors (CD34+/CD19−/IgM−), pro-B (CD34+/CD19+/IgM−), pre-B (CD34−/CD19+/IgM−), and immature B (CD34−/CD19+/IgM+) populations (Figure 1A top and Supplementary Figure S1a). We then extracted RNA, performed RNA-seq, and quantified transcript levels. Concordant with flow cytometric profiles, CD19 mRNA levels increased throughout differentiation stages, CD34 mRNA was confined to early progenitors and pro-B fractions, and IgM transcript was expressed only in the immature fraction (Figure 1A, bottom). Furthermore, we performed principal component analysis (PCA) on all expression datasets and found tight clustering of the four fractions from different donors, suggesting that at the level of mRNA expression B-cell differentiation supersedes individual variations (Figure 1B). To ensure that the clustering of fractions was not driven solely by CD19 and CD34 expression, we repeated the PCA, removing CD19 and CD34 expression contributions and observe very similar clustering of fractions (Supplementary Figure S1B).

Figure 1.

RNA-seq analysis of bone marrow-derived B-cell progenitors and pediatric B-ALL. (A) Top: Lymphocytes were isolated from normal whole bone marrow aspirates and fractionated into CD34+/CD19-/IgM-, CD34+/CD19+/IgM-, CD34-/CD19+/ IgM-, and CD34-/CD19+/IgM+ populations. Bottom: Validation of surface marker expression by RNA-Seq in early progenitors, pro-B, pre-B and immature B cells fractions. (B) Principal component analysis (PCA) of RNA expression data from bone marrow fractions obtained from 2 adult (solid circles) and two pediatric (open circles) donors. Cumulative variances are shown for each PC. (C) Quantification of CD19 and CD34 expression in bone marrow fractions and B-ALL samples. (D) PCA on RNA expression data from 18 B-ALL samples and four bone marrow samples.

RNA-Seq analysis of primary B-ALL samples

We obtained 24 primary pediatric B-ALL samples from the CHOP Center for Childhood Cancer Research leukemia biorepository. Ficoll-purified mononuclear cells from fresh bone marrows or peripheral blood were used for total RNA extraction. Based on the presence of intact RNA (evident by 28S and 18S bands on gels), an RNA integrity score >5.2, and RNA concentration of >40ng/µL, we successfully extracted high quality RNA for sequencing from 18 out of 24 frozen samples (Supplementary Figure S1C, red asterisks indicate samples that did not pass QC and were not sequenced). These samples were comprised of several different phenotypic and genetic B-ALL subtypes (31) (Supplementary Table S1 and Supplementary Figure S1D). We performed RNA-seq of these leukemia samples and first compared them to adult and pediatric normal bone marrow fractions with respect to CD19 and CD34 expression. We found that B-ALL samples closely resembled the pro-B fractions in that CD19 and CD34 transcripts were readily detectable in the low-to-medium range (Figure 1C). To extend this analysis to the entire transcriptome, we performed PCA on B-ALL samples and four sets of normal bone marrow fractions (two pediatric and two adult). Once again, the 18 B-ALL samples clustered most closely with the pro-B fractions (Figure 1D). This similarity was reflected in shortest Euclidian distances between PC1 coordinates for the centroids of the B-ALL and pro-B cell samples (Supplementary Figure S1E). Therefore, we chose the pro-B fractions as cell-of-origin controls for B-ALL splicing analysis.

Global patterns of aberrant splicing in pediatric B-ALL

To detect patterns of alternative splicing in B-ALL, we utilized the MAJIQ (Modeling Alternative Junction Inclusion Quantification) algorithm (25). MAJIQ offers the ability to detect, quantify, and visualize complex splicing variations, including de-novo variations, while comparing favorably in reproducibility and false discovery rate to common alternatives (25,32). Using MAJIQ, we performed 18 independent pairwise comparisons between the average of pediatric pro-B cell fractions and leukemia samples. To assess heterogeneity in splicing across samples, we measured the number of differential LSVs (minimum of 20% change in splicing, 95% CI—see Materials and Methods) in each sample and compared their identities. We observed that each B-ALL specimen had >3000 LSVs but most LSVs were either unique (first bar in Figure 2A) or shared by only a small number of B-ALL patient samples (subsequent bars in Figure 2A). However, we also found a total of 279 aberrant LSVs in 241 genes that were detected in all of our 18 pairwise comparisons (last bar in Figure 2A). Provocatively, when this 241-gene set was analyzed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) (33), the top gene ontology (GO) categories (34) were related to RNA splicing (Figure 2A, inset). When we investigated individual SF transcripts, we found that of the 167 well-characterized SF genes [Supplementary Table S2, from (35)], 101 were alternatively spliced in B-ALL compared to pro-B cells (Figure 2B, bottom left, proB- versus B-ALL bar chart). Moreover, these changes were highly specific to the malignant phenotype and not to states of normal B-cell differentiation (Figure 2B, upper and middle panels). For example, only six SF transcripts were alternatively spliced during the early progenitors-to-proB transition and only nine during the proB-to-preB transition (Figure 2B, middle).

Figure 2.

Global patterns of aberrant splicing in pediatric B-ALL. (A) Enumeration of local splicing variations (LSVs). The vertical bar chart shows the numbers of LSVs detected by MAJIQ per B-ALL sample compared to normal pro-B cells. The horizontal bar chart (insert) showing enrichment for RNA splicing pathways among genes with most consistent LSVs. (B) Overlap of LSVs in SF genes between B-ALL samples and bone marrow fractions. The left panel shows the different sets and number of genes included in them. The middle section visualizes the different combinations of intersection (including each set by itself). The upper panel shows how many genes of the total set are in each intersection. (C) Heat maps showing changes in inclusion of exons corresponding to SRSF and hnRNP families members in 18 B-ALL samples compared to pro-B cell fractions. Green color denotes increased inclusion, red color – increased skipping. (D) Venn diagram showing overlap of LSVs from hnRNP and SRSF families when B-ALL samples are compared to in-house pediatric pro-B, in-house adult pro-B, and publicly available bone marrow data (GSM1695856 and GSM1695857), as indicated. (E) Differential expression analysis of hnRNP and SRSF family members to decouple detection of LSVs with changes in RNA expression. Majority of samples (Y-axis) showed constant expression of indicated genes compared to normal pro-B samples (green). Several samples showed downregulation of a few genes (red), such as HNRNPA1.

For each LSV, MAJIQ provides a ΔPSI (percent-spliced-in) value to indicate changes in isoform abundance. We observed that many members of the hnRNP and SRSF families, which play key roles in exon inclusion or skipping (36), exhibited profound changes in ΔPSI values (Figure 2C). These 25 detectable LSVs included increases in the so-called ‘poison exons’ with in-frame stop codons in several SRSF proteins (37). Of note, some transcripts contained multiple aberrant LSVs. For example, SRSF11 had three distinct LSVs with different start/end UCSC Genome Browser coordinates.

To increase confidence in the detection of splice factor related LSVs in B-ALL compared to normal pediatric pro-B samples we also compared B-ALL to normal adult pro-B cells and to publicly available RNA-seq data corresponding to the CD19-positive normal bone marrow cells (BM-GEO) (38). We overlapped these comparisons in a three-way Venn diagram to show that of the 25 LSVs detected in B-ALL versus pediatric pro-B 21 LSVs were also detected in the two other independent and publicly available datasets [B-ALL vs. Adult pro-B and B-ALL vs. BM-GEO (38)] (Figure 2D). Further, we searched for LSVs in the SRSF and hnRNP families in larger, independent TARGET (39) and St Jude Children's Research Hospital (SJCRH) (40) B-ALL datasets using pediatric pro-B cells for comparison. We detected 75 and 78 LSVs in the TARGET and SJCRH datasets, respectively (Supplementary Figure Ss2, left bars, middle and bottom). Of these LSVs, 43 were exclusive to these datasets (Supplementary Figure S2A, third column labeled 43 and 2 connected dots below for TARGET and SJCRH) and 22 were detected in all three datasets (Supplementary Figure S2a, middle, first column labeled 22 and 3 connected dots below for CHOP, TARGET and SJCRH) (see also Supplementary LSV Lists).

To decouple observed changes in splicing from changes in gene expression we performed differential expression analysis on genes with frequently detected LSVs. We found that most genes containing LSVs did not have associated changes in expression, with the notable exception of HNRNPA1, which was observed to be downregulated at least two-fold in approximately half of the B-ALL samples (Figure 2E). Furthermore, after lowering the threshold to 60% change, we found that all 18 samples showed decreased expression relative to proB cells. This interesting expression pattern and the pervasive nature of HNRNPA1 splicing alterations (Figure 2C, right) prompted us to investigate its regulation further.

Dysregulated splicing of HNRNPA1 in B-ALL

According to the Ensembl database (26), there is evidence for HNRNPA1 transcripts with alternatively spliced 3′ UTRs, notably the canonical long proximal exon 11 and two shorter, distal exons (exons 12/13, ENST00000547566) (Figure 3A, top). We observed that the exon 11-containing transcript predominated in normal pro-B cells (Figure 3B, red), but in a typical B-ALL specimen (representative sample PHL50) the predominant event was the skipping of exon 11 to exons 12 and 13 (Figure 3B, blue). To ensure validity of this event, we also ran the LeafCutter splicing algorithm to detect splice events present in the averaged RNA-seq datasets from 18 B-ALL samples and four pro-B samples. Confirming our MAJIQ analysis we were able to detect alternative splicing of HNRNPA1 3′ UTR using this independent analysis (Supplementary Figure S2B). While exon 11 percent-spliced-in (PSI) values varied across leukemia samples, most leukemias had increased skipping of exon 11 compared to pro-B samples (Figure 3C, blue). We also observed intron 10 inclusion in all pro-B and B-ALL samples, although this event was not significantly different between normal and malignant samples (Figure 3C, gray stacks).

Figure 3.

Alternative splicing of HNRNPA1 in pediatric B-ALL. (A) Exon–intron structure of the HNRNPA1 transcript and location of RT-qPCR primers. (B) Violin plots depicting alternative splicing of the HNRNPA1 transcript. (C) Stack plot of PSI values of exon 11 corresponding to pro-B cell fractions and 18 B-ALL samples. (D) Analysis of exon 11 inclusion in B-ALL samples by RT-qPCR. (E) Spearman's correlation of RT-PCR vs. MAJIQ (rho = 0.06, P = 0.00737). (F) Spearman's correlations of expression of an HNRNPA1 constitutive exon and alternative exon 11 measured by RNA-Seq (rho = 0.82, P = 2E–05). (G) Measurement of mRNA stability of HNRNPA1 3′ UTR isoforms (proximal exon 11 in gray squares and distal exon 12 in black circles) following actinomycin D treatment at indicated timepoints. The dashed line indicates the half-life of exon 12-included isoform of approximately 1.5 h. (H) Relative expression of alternative HNRNPA1 3′UTR exons 11 or 12 in Nalm6 cells after treatment with control siRNA (black) or siRNAs targeting UPF1 and UPF2 (gray) for 48 h.

We then validated and quantitated exon 11 inclusion by RT-qPCR using a forward primer spanning exons 9–10 and a reverse primer in exon 11 (Figure 3A, D). Using Spearman's rho statistic we observed a strong positive association between MAJIQ predictions and RT-qPCR validations (0.61, P-value = 0.00737) (Figure 3E). We next applied the same correlation analysis to HNRNPA1 mRNA expression levels versus exon 11 inclusion and found a positive correlation between the two measurements based on both RNA-seq (0.82, P-value = 0.00002) (Figure 3F) and RT-qPCR (0.6, P-value = 0.00878) (Supplementary Figure S2C). This suggested that preferential splicing from exon 10 to the distal 3′ UTR exons 12/13 (skipping exon 11) may decrease HNRNPA1 RNA steady state levels. To assess the half-lives of HNRNPA1 isoforms we performed actinomycin D time course treatment followed by RT-qPCR. We observed strong differences in mRNA stability between the UTR isoforms. Specifically, the half-life of the exon 11-skipped isoforms was ∼1.5 h (Figure 3G) while the half-life of the exon 11-included isoform was longer than 24 h (Figure 3G and Supplementary Figure S2D).

Based on the exon junction present between exons 12 and 13 downstream from the stop codon we hypothesized that this transcript might be a target of nonsense-mediated decay [NMD; (41,42)]. To test this, we first treated Nalm6 B-ALL cells with the translation inhibitor cyclohexamide to transiently suppress NMD (41,43,44). To show effective suppression of NMD, we performed gel-based PCR for primers spanning SRSF3 exon 4, a canonical target of NMD (45) (Supplementary Figure S2E). We then performed RT-qPCR for HNRNPA1 3′ UTR exons and observed accumulation of the transcript that skips to exon 12/13, but not the transcript including exon 11 (Supplementary Figure S2F).

To further confirm the exons 12/13 transcript as an NMD substrate, we co-knocked down two factors required for NMD, UPF1 and UPF2 (46,47), in Nalm6 B-ALL cells (Supplementary Figure S2g). Despite modest co-knockdown of the factors (∼35% at most), quantification of the canonical NMD substrate SRSF3 exon 4 by gel-based PCR confirmed inhibition of NMD at the highest concentrations of siRNA for UPF1 and UPF2 (750 nM each) (Supplementary Figure S2H). Next, we performed RT-qPCR for HNRNPA1 3′ UTR exons and observed preferential accumulation of the exon 12/13 transcript as compared to the exon 11 transcript (Figure 3H) confirming that only the former is a robust target of degradation by the NMD machinery.

To model downregulation of this splicing factor and to identify affected transcripts, we knocked down HNRNPA1 mRNA with siRNA in the Epstein-Barr Virus (EBV)-transformed human B-lymphoblastoid P493-6 cell line (30). Both the non-targeting scrambled control and the HNRNPA1-targeting reagent were Smart Pools of four siRNAs, designed to target distinct sites within the specific gene of interest. This approach allowed us to reduce the reliance on and the concentration of each individual siRNA and average out their off-target effects, as documented in (48). Knockdown of HNRNPA1 was confirmed by RT-qPCR (Figure 4A) and Western blotting (Figure 4B). Following RNA-Seq and MAJIQ analysis, we identified 213 LSVs in 184 genes associated with knockdown of HNRNPA1 (minimum of 20% ΔPSI, 95% CI). Of these hnRNPA1-dependent LSVs, 74 (or >30%) were present in at least one B-ALL sample (Figure 4C), with more than half of these 74 LSVs altered in 10 or more samples. DAVID analysis of the genes with overlapping LSVs revealed enrichment for macromolecular metabolic pathways (Figure 4D), supporting the well-established role of hnRNP proteins in RNA metabolism (49). We next searched for the presence of the hnRNPA1 binding motif (50) in exons and flanking intronic sequences (±200nt) involved in overlapping LSVs. We were able to identify the UAGGG motif in 68 out of 74 LSVs in question (Figure 4E, top).

Figure 4.

hnRNPA1-regulated splicing of DICER1. (A) Expression of HNRNPA1 RNA in P493-6 cells as determined by RT-qPCR (ex2–3–4 primers). A pool of four non-targeting control (black) and four HNRNPA1-targeting (grey) siRNAs were used in duplicate experiments. (B) hnRNPA1 protein expression determined by immunoblotting with anti-hnRNPA1 and anti-tubulin antibodies in cells from panel A. (C) Overlap in LSVs detected in B-ALL versus pro-B and HNRNPA1 knockdown versus control cells. (D) The horizontal bar chart showing enrichment for macromolecular metabolic pathways among genes with LSVs detected in both HNRNPA1 knockdown and B-ALL. (E) hnRNPA1 motif analysis of exons and adjacent intronic sequences (±200nt) surrounding LSVs detected in HNRNPA1 knockdown and B-ALL samples (top). Significant hnRNPA1 motifs found ∼200 nts upstream of DICER1 exon 1 (bottom). (F) Violin plots depicting alternative splicing of DICER1 in Pro-B vs. PHL149 B-ALL (top) and HNRNPA1 control and knockdown P493-6 cells (bottom). (G) Expression of HNRNPA1 RNA in P493-6 cells as determined by RT-qPCR (ex2-3-4 primers). A pool of four non-targeting control (black) GapmeRs and HNRNPA1-ATG-targeting GapmeR (grey) were used in triplicate experiments. (H) hnRNPA1 protein expression determined by immunoblotting with anti-hnRNPA1 and anti-tubulin antibodies in cells from panel G. (I) DICER1 alternative exon 4 inclusion measured by RT-qPCR following transfection with control or HNRNPA1-targeting GapmeRs.

Among these genes characterized by overlapping LSVs, involvement in macromolecular metabolism, and the presence of hnRNPA1 binding motifs, was DICER1 (Figure 4E, bottom, purple arrows indicate hnRNPA1 motifs). Interestingly, one of the LSVs in the DICER1 gene was altered in B-ALL samples in the same manner as in the HNRNPA1 KD cells (Figure 4F). We also confirmed the presence of this LSV using LeafCutter algorithm (Supplementary Figure S3A). Notably, this LSV maps to the 5′UTR of the DICER1 transcript, potentially diminishing its translation efficiency. We validated this LSV in Nalm6 B-ALL cells by RT-PCR and Sanger sequencing and noted that this type of deregulation could phenocopy loss-of-function mutations and copy number variations in DICER1 in several types of cancer including leukemias (51).

Since Dicer itself is involved in the RNAi pathway (52), we were concerned that this involvement might affect the results of our siRNA experiments. Thus, we knocked-down HNRNPA1 in a non-RNAi way by using the GapmerR technology (53). To ensure that the knock-down was well-controlled, we custom-designed a pool of four GapmeRs with non-targeting scrambled sequences and compared it to the GapmeR targeting the translation start codon of HNRNPA1. The efficiency of knock-down was validated by RT-qPCR (Figure 4G) and Western blotting (Figure 4H). We were then able to confirm the increased inclusion of DICER1 alternative 5′ UTR exon 4 by RT-qPCR following HNRNPA1 knock-down (Figure 4I).

Aberrant splicing of leukemia drivers in pediatric B-ALL

We next aimed to identify alternative splicing in other genes that contribute to leukemogenesis. We retrieved all genes (141) with known somatic mutations in hematologic malignancies from the COSMIC database (v.82) (20,21) (Supplementary Table S3) and searched for aberrant LSVs affecting these genes in B-ALL samples using the very stringent ΔPSI >50% cutoff for exon inclusion/skipping. We identified 81 aberrant LSVs in 41 genes that were present in at least two B-ALL samples. They accurately separated leukemia samples from the non-transformed pro-B cell counterparts following hierarchical clustering using Euclidian distances (Figure 5A). Of note, these LSVs affected roughly a third of B-ALL driver genes. Some (in genes such as FBXW7) were present in almost all B-ALL samples, attesting to their potential significance in leukemia pathogenesis. Again, to decouple changes in splicing from changes in expression we performed differential expression analysis on genes containing LSVs in the COSMIC genes. Consistent with SF data analysis, the majority of genes in the majority of samples did not have changes in expression associated with changes in splicing (Supplementary Figure S3A).

Figure 5.

Alternative splicing of leukemia drivers in pediatric B-ALL. (A) Heat maps showing changes in splicing in B-ALL of genes commonly mutated in hematologic malignancies. The dendrogram on top of the heatmap represents results of hierarchical clustering. The dotted blue line denotes LSVs in the DICER1 transcript and also demarcates the pro-B samples. (B) Overlap between the LSVs from panel (A) and those detected in TARGET and SJCRH B-ALL datasets. (C) Bar graph representing frequencies of mutations (COSMIC dataset) and major LSVs as defined by MAJIQ in the indicated datasets. (D) Violin plots depicting alternative splicing of NT5C2 exon 4a (red) in control siRNA- vs. HNRNPA1 siRNA-electroporated P493-6 cells (top). The same analysis was applied to pro-B cells vs. the representative B-ALL sample SJBALL021849. (E) NT5C2 alternative exon 4a inclusion as measured by RT-qPCR following transfection with control (black) or HNRNPA1-targeting (grey) GapmeRs.

To confirm our findings in independent datasets, we searched for aberrant LSVs in the COSMIC genes in the TARGET (39) and SJCRH (40) B-ALL datasets using pediatric pro-B cells for comparison. We found 229 LSVs in 80 genes in the TARGET data and 362 LSVs in 95 genes in the SJCRH data (Figure 5B). Importantly, 64 of 81 LSVs identified in the CHOP dataset significantly overlapped with TARGET (Fisher's Exact Test, P = 7e–124) and SJCRH datasets (Fisher's Exact Test, P = 1e–139), 13 LSVs were present only in CHOP and SJCRH datasets and another two LSVs were present only in CHOP and TARGET datasets (Figure 5B). We then narrowed down our analysis to LSVs in top 20 B-ALL tumor suppressors and oncogenes as defined in the COSMIC database. In both TARGET and SJCRH datasets, we discovered largely overlapping LSVs corresponding to 15 out of the top 20 mutated genes in B-ALL. Of note, LSV frequencies were much greater than frequencies of somatic mutations (Figure 5C), suggesting that in B-ALL these driver genes are preferentially affected by post-transcriptional mechanisms, such as alternative splicing. For example, we detected increased inclusion of NT5C2 alternative exon 4a in ∼65% of SJCRH samples analyzed (Figure 5D, top, red trace and violin plot). Apparently, this event is also under the regulation of HNRNPA1, as evidenced by increased inclusion of exon 4a in both HNRNPA1 knock-down experiments (siRNA in Figure 5D, GapmeR in Figure 5E). Additionally, we detected increased inclusion of TP53 exon 9i in B-ALL compared to pro-B (Supplementary Figure S3B), which would promote the expression of the SRSF3-regulated p53β isoform (Supplementary Figure S3C) involved in cell cycle and cell death (54–57).

DISCUSSION

Personalized cancer diagnostics employ carefully selected next-generation sequencing panels, which can identify mutations in specific genes known or suspected to be drivers in human malignancies. Hematologic malignancy panels typically include dominant oncogenes (e.g. FLT3 and IL7R), recessive tumor suppressors (e.g. TP53 and FBXW7), and haploinsufficient DNA/RNA caretakers (e.g. splicing factor SRSF2). Ostensibly, genomic profiling of diagnostic cancer specimens should be able to identify all prognostic mutations and inform treatment selection for patients. However, our data demonstrate that in B-ALL deregulation of cancer genes occurs primarily at the level of splicing, rather than via genetic mutations. For example, MAJIQ analyses of B-ALL transcriptomes revealed that several SRSF genes controlling exon inclusion (36), while unmutated, show widespread variations in their own splicing patterns, some of which are bound to decrease protein levels (37). Consequently, we observed aberrant splicing of some of their known target transcripts such as TP53 (57), which encodes the key tumor suppressor gene. Of note, the majority of aberrant LSVs were highly concordant across different datasets. This reproducibility validates our conclusions and alleviates the potential concern that sample preparation conditions (e.g. storage at ambient temperature) could be affecting RNA surveillance and thus impacting analysis of alternative splicing (58).

One of the most consistent changes in exon usage was the non-canonical selection of 3′ UTR exons of HNRNPA1, a splice factor implicated in cancer progression (59,60). While selection of alternative poly(A) sites is emerging as an important mechanism of gene expression in leukemia (61), that mechanism is distinct from alternative splicing involved in HNRNPA1 regulation. Indeed, as transcription and splicing occur on similar timescales (62), with splicing estimated to be half-complete when Pol II is less than 50nt downstream of introns (63), the decision to utilize 3′ splice site in exon 11 or 12 of HNRNPA1 has to be made prior to RNA polymerase II approaching either of the alternative poly(A) sites. Regardless of the underlying mechanism, while HNRNPA1 is overexpressed in some cancers, in the B-ALL model HNRNPA1 LSV correlated with a decrease in HNRNPA1 mRNA abundance. This is in agreement with our data showing that this multi-exon 3′UTR triggers non-sense mediated decay of the transcript. It is also possible that new 3′UTR sequences could create additional sites for targeting by microRNAs, which are known to play a key role in HNRNPA1 downregulation in chemotherapy-resistant ovarian cancer cells (64). Interestingly, knockdown of HNRNPA1 in an Epstein-Barr Virus (EBV)-transformed human B-cell lymphoblastoid cell line resulted in aberrant splicing of Dicer, a key enzyme in microRNA biogenesis [reviewed in (52)]. Moreover, four non-hnRNPA1-dependent LSVs present in B-ALL samples (Figure 5A) affected Drosha, another key enzyme in the microRNA pathway (52), suggesting a strong link between splicing and miRNA machineries.

Of even greater importance is the fact that B-ALL-specific LSVs affect 15 out of 20 top leukemia driver genes (including the aforementioned NT5C2 and TP53), with frequencies far exceeding those of somatic mutations. This discovery could explain why the prevalence of somatic mutations and copy number variations in B-ALL is low compared to other human cancers. It remains to be determined whether splicing alterations in oncogenes and tumor suppressor genes are functionally equivalent to gain-of-function and loss-of-function mutations, respectively. If this proves to be the case, interfering with splicing using RNA-based therapeutics and/or available small molecule inhibitors could be used to inhibit oncogenes such as FLT3 and to activate dormant tumor suppressor gene such as TP53. Such strategies could yield tangible therapeutic benefits across a broad spectrum of childhood B-ALL subtypes.

DATA AVAILABILITY

The original RNA-Seq data sets have been deposited in the NCBI GEO database under accession number GSE115656.

Click here for additional data file.