Dissecting the tumorigenesis trajectory by single-cell RNA-seq reveals critical roles of RNA splicing

doi:10.21203/rs.3.rs-2105579/v1

Download PDF

Research Article

Dissecting the tumorigenesis trajectory by single-cell RNA-seq reveals critical roles of RNA splicing

https://doi.org/10.21203/rs.3.rs-2105579/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Leukemogenesis is proposed to be a multistep process by which normal hematopoietic stem and progenitor cells are transformed into full-blown leukemic cells, the details of which are not fully understood.

Methods

We performed serial single-cell transcriptome analyses(10X genomics) of preleukemic and leukemic cells (PLCs) and constructed the cellular and molecular transformation trajectory in a Myc-driven AML model in mice, which represented the transformation course in patients. Another single-cell approach, Smart Seq2, was used to analyze alternative splicing events in leukemic cell.

Results

We found that the Myc targets were gradually upregulated along the trajectory. Among them were splicing factors, which showed stage-specific prognosis for AML patients. Furthermore, we dissected the detailed gene network of a tipping point for HSPCs to generate initiating PLCs, which was characterized by dramatically increased splicing factors and unusual RNA velocity. In the late stage, PLCs acquired explosive heterogeneity through RNA alternative splicing. Among them, the Hsp90aa^hi subpopulation was conserved in both human and mouse AML and associated with poor prognosis. Exon 6 skipping of Tmem134 was identified in these cells. While the exon skipping product Tmem134β promoted the cell cycle, full-length Tmem134α delayed tumorigenesis.

Conclusion

Our study emphasized the critical roles of RNA splicing in the full process of leukemogenesis and identified the role of a novel tumor specific alternative splicing event, Exon 6 skipping of Tmem134, during the occurrence and maintenance of leukemia.

Tumorigenesis, a process to accumulate enough genetic alterations for normal cells to be transformed into malignant cells, can take decades in patients(1–5). Similarly, leukemogenesis is conceptualized as the multistage process for normal hematopoietic stem and progenitor cells (HSPCs) to become full-blown leukemic cells (6, 7). It has been shown that largely normal clonal hematopoiesis (CH) can progress into low-risk myelodysplastic syndrome (MDS) and then high-risk MDS or myeloproliferative neoplasm (MPN), eventually leading to full-blown acute myeloid leukemia AML(7). However, AML is among the human malignancies with the lowest mutation burdens(8), and there is accumulating evidence suggesting that the majority of genetic alterations exist in preleukemic conditions such as CH (clonal hematopoiesis) and MDS(9–12). Therefore, we wondered whether molecular rewiring other than mutations plays a significant role during leukemogenesis. Recent advances in omics studies, especially single-cell RNA sequencing (scRNA-seq), combined with our murine AML model, have provided unique opportunities to reveal the trajectory of leukemogenesis and further molecular characterization after acquiring initiating mutations.

In this study, we constructed a transformation route for HSPCs into full-blown AML by analyzing the transcriptomes of 18,900 preleukemic and leukemic cells at different stages and revealed progressively abnormal RNA splicing during leukemogenesis.

Mice

C57BL/6 mice (Jackson Lab, Cat# 000664) were used as HSPC donor mice and transplantation recipients. CAG-Cas9-EGFP mice (Cat# 026179) were purchased from the Jackson Laboratory.

Murine AML model and AML cell line

Whole bone marrow (WBM) cells were obtained from 6-week-old C57/BL6 mouse tibiae, femora and ilia bones and stained with mCD117 magnetic microbeads (Miltenyi Biotec, Cat# 130-091-224). Then, cKit + HSPCs were enriched from MidiMACS™ Separator (Miltenyi Biotec, Cat# 130-042-302) and cultured in BCM medium (50% DMEM + 50% IMDM medium) supplemented with 20% FBS, 0.34% beta-mercaptoethanol, 10 ng/mL mIL3/mIL6, and 50 ng/mL mSCF. After 24 hrs of cultivation in the incubator, the retrovirus harboring pMIG-Myc was spin infected and centrifuged at 2000 rpm for 1 hour. After 24 hours, cKit + HSPCs with an ~ 20% infection ratio were injected via the tail vein into sublethally irradiated (5.5 Gy) 7-week-old female wild-type C57/BL6 mice (1*10^6 cells per mouse). Bone marrow-transplanted mice were maintained on medicated water (ciprofloxacin, 100 µg/mL) for 2 weeks and monitored for leukemogenesis. After 7 ~ 8 weeks, AML cells were harvested from the bone marrow or spleen and cultivated into a cell line in BCM medium.

Cas9 AML cell line

pLentiCRISPV2 (Addgene plasmid #52961) was transduced into AML cells by lentivirus generated by the 293T cell line, and then Cas9-expressing AML cells were selected by puromycin administration (1 µg/mL) for 1 week.

Retroviral constructs

In Fig. 6E-F, the cDNA elements overexpressing TMEM134α and TMEM134β were cloned into pMSCV-IRES-GFP-Myc (pMIG-Myc). In Fig. 6G, 6I, S8N and S8O, the cDNA elements overexpressing MYC, TMEM134α, and TMEM134β (synthesized by Sangon Biotec, Shanghai) were cloned into pMSCV-IRES-GFP (pMIG) as previously reported(13) and pMSCV-IRES-mCherry (pMIC, constructed in house) retroviral vectors harboring a GFP/mCherry fluorescent reporter. Retroviruses were generated in the 293T cell line, and cDNA expression constructs were transduced into the cKit + HSPC/AML cell line by a retrovirus.

CRISPRs

CRISPRs were designed at https://www.atum.bio and then cloned into V2TC (modified by pLentiCRISPV2; briefly, the Cas9 component in the vector was replaced with a mCherry fluorescent reporter). The target sgRNA sequences are shown in Supplementary Table 9. Lentiviruses were generated in 293T cell lines, and sgRNAs were transduced into Cas9-expressing cKit + HSPCs/AML cell lines by lentivirus. The PCR primers are shown in Supplementary Table 10.

Cell growth curve

Cells were seeded at a concentration of 1*10^4 or 2*10^4 cells per well of a 96-well plate with three wells of each line. Cells were counted every 2 days, and subsequent cell counts were normalized to the seeding density after each count. All curves were obtained at least 2 independent times.

Cell proliferation and cell differentiation

Cell proliferation assays were performed with a Cell-Light™ EdU Apollo488 In Vitro Kit (RIBOBIO, Cat# C10310-3). Briefly, cells were infected with viruses encoding Tmem134α or Tmem134β cDNA or empty vector. On day 3 after infection, cultured cells were stained with 1:4000 EdU staining solution for 2 hours, and then the cells were washed and fixed in 4% PFA. Cells were then incubated with 0.5% Triton X-100 and stained with Apollo^® reaction buffer.

For the cell differentiation experiment, Cas9-expressing HSPCs were infected with viruses encoding Myc and sgRNAs targeting scramble or Npm1. On days 2, 4 and 6 of infection, cultured cells were analyzed by staining for the differentiation surface marker Mac-1 (Biolegend, clone: M1/70, Cat# 101223)/Gr-1 (Biolegend, clone: RB6-8C5, Cat# 108430) and the stemness surface marker c-Kit (Biolegend, clone: 2B8, Cat# 105814). Double-positive and single-positive populations were gated following standard flow analysis to detect c-Kit and Mac-1/Gr-1 expression levels.

All flow cytometry analyses were performed with an LSRFortessa instrument (BD Biosciences).

RNA-seq analyses

RNA-seq libraries were constructed by an Illumina Stranded mRNA Sample Preparation Kit (NEB, Cat# E7770) according to the manufacturer’s protocol and were sequenced by an Illumina NovaSeq 6000 sequencing machine with 150-bp paired-end reads. The RNA-seq reads were aligned to the mouse reference genome (GRCm38) by STAR_2.6.0a. Transcript abundance was normalized and measured by transcripts per kilobase million (TPM). Gene set enrichment analysis (GSEA) was applied to show statistically significant similarities and differences between two given clusters by identifying a priori defined gene sets.

SMART-seq2 library preparation and sequencing

SMART-seq2 was performed on AML cells at T₃, following a modified SMART-seq2 protocol published previously(14). The libraries were sequenced by an Illumina HiSeq X™ Ten Sequencing machine with 150-bp paired-end reads. The paired-end reads were aligned to the UCSC mm10 references by STAR_2.6.0a.

Single-cell collection, library preparation and sequencing

Whole Bone Marrow (WBM) were isolated from mouse tibiae, femora and ilia leg bones by crushing in 2%FBS(PBS) buffer, add 1:1000 DNase1 and then translate into 15 mL BD tube, centrifuge at 1200rpm, RT for 3 min. Then, the supernatant was discarded, 3 mL of 1x ACK lysis buffer was added to resuspend the bone marrow pellet on ice for 3 ~ 5 min, and PBS (more than lysis buffer) was added to terminate the lysis reaction. The supernatant was discarded, and 0.04% PBS-BSA buffer was added to resuspend the cell pellet, which was then filtered through a 40 µm cell strainer. More than 80% of the living cells in the sample and a density of living cells ranging from 600 to 1200 cells/µl were preferable. Sample quality was measured using a Bio-Rad TC20 automated cell counter. Libraries were prepared using Chromium Single Cell 3' Reagent Kits v2 according to the manufacturer’s protocol (10x Genomics) and were sequenced using Illumina NovaSeq 6000 sequencing.

Single-cell gene expression quantification

The ectopic Myc sequences were added to the reference genome, which could distinguish exogenous and endogenous cells in AML leukemogenesis. Cellranger (v2.1.1) was used to align clean reads with mm10. The Seurat (v.3.2.0) pipeline was applied to analyze scRNA data. Genes expressed in less than 5 cells were not considered, and cells that contained lower than 10% mitochondrial genes and contained at least 200 but not more than 6,000 genes were retained for the following analysis. Normalized data and scale data were generated from the Seurat workflow.

Batch effect reduction, dimension reduction and unsupervised clustering

To integrate all the single cells from different samples for unsupervised clustering, integration anchor algorithms in Seurat were used to reduce the batch effect. The top 4000 variable genes were used as highly variable genes by FindIntegrationAnchors, and an integrated assay matrix was used to build the nearest neighbor graph.

For dimension reduction, UMAP, t-SNE, diffusion-map and force-directed graph were calculated on a batch-reducing matrix. UMAP and t-SNE were calculated by Seurat, and diffusion maps and force-directed graphs were calculated by URD. The cell subtypes were identified by the classic markers in previous studies(15, 16).

To unify the identities of cell subtypes from a different platform, we used TransferData in Seurat. The 10X genomics data were used as the reference files, and smartseq2 data were used as query files to carry out transfer label pipelines. Then, unsupervised clustering algorithms were used to identify each cell identity in smartseq2.

Trajectory construction and pseudotime calculation

To calculate the trajectory and pseudotime on the UMAP plot of tipping point cells, the position matrix of scRNA data was extracted and processed to feed on Slingshot. Three hundred fixed points were performed to fit an approximate trajectory. In addition, the general additive model was used to identify the genes with dynamic expression based on pseudotime. The force-directed layout and pseudotime of PLCs were determined by URD packages.

GSVA, gene regulatory network analysis and molecular interaction network analysis

GSVA was used to calculate the enrichment score per cell or sample. ClusterProfiler was used to enrich the GO and KEGG pathways. Pathway network analysis and visualization were performed on the modified cnetplot function. The molecular interaction network was calculated by RSCORE, with mouse-specific protein–protein interaction networks downloaded from BioGRID. Gene interaction matrix, generated from the gene expression matrix by R. SCORE function. The gene interaction network was visualized by the Steiner tree. The genes with highly significant differential expression and genes with highly significantly different network interactions are highlighted in the network map.

Generation of gene signatures and visualizations

Cell-type-specific gene signatures were generated by the FindMarker/FindAllMarker function implemented in Seurat. Stemness gene signatures and differentiation gene signatures were generated from HSPC-specific expressed genes or neutrophil-specific expressed genes with significant statistical power. Gene sets of Myc targets and splicing factors were downloaded from the Molecular Signatures Database (MSigDB). The mean value of the log2-transformation matrix was used to quantify the global activation of gene signatures or pathways in each cell/sample. The signature score was projected on the UMAP plot by using the modified FeaturePlot function.

Splicing events calculation

rMATS v4.0.2 was used to analyze splicing events on Smart-seq2 single-cell RNA-seq data. bam files were divided into two groups based on the cluster information. We summarized all splicing event types of rMATS results, including SE, A3SS, A5SS, RI and MXE.

Cell cycle score calculation

To calculate the cell cycle score in each single cell, cell cycle phase signatures for humans and mice were obtained from the Seurat pipeline(17). In each cell of single-cell data, CellCycleScoring with nbin = 28 was used to calculate G2 M scores and S scores. Cells with positive scores were divided into the G2 M phase or S phase, and their G1 scores were assigned to 0. For cells with negative scores, G1 scores were calculated by the following formula.

$$G1 Scores=\left(-1\right)*\left(\frac{G2M Scores+S Scores}{2}\right)$$

After the G1 scores were calculated, the R package ggpubr was used for visualization. To calculate G1 scores in each bulk RNA-seq sample, the z score of the mean TPM value of G2 M signature scores and S phase signatures was evaluated.

Bulk data processing and signature identification in MDS and AML patients

Low-risk, high-risk MDS and MDS_AML patient data GSM2430698(18) and AML patient data GSE1159(19) were downloaded from the GEO database. We processed MDS patient data by the RNAs-seq workflow as previously mentioned. DESeq2 was used to identify the differentially expressed genes (DEGs). Three types of signatures were obtained, including low-risk MDS signatures generated by comparing low-risk samples to normal samples (p.adj < 0.05 and log2FoldChange > 3), high-risk MDS signatures (p.adj < 0.05 and log2FoldChange > 4) and MDS-AML signatures (p.adj < 0.05 and log2FoldChange > 4) generated by the same kind of comparison. For AML patient data, DESeq2 was used to identify the DEGs (p.adj < 0.05 and log2FoldChange > 1.5). To project the patients’ signature on mouse data, we used the R package biomaRt to convert the human genes to homologous mouse genes.

TCGA data analysis

Clinical information and transcriptome data of TARGET-AML and TCGA-AML were downloaded by the R package TCGAbiolinks. The percent-splice-in (PSI) value of TCGA-LAML was downloaded from TCGA SpliceSeq. Then, the clinical data were processed by filtering patients who could not match the transcriptome and clinical information. Patients whose follow-up days and days to death were not unknown were filtered.

To test the correlation of selected signature expression and patient survival, we calculated the signature scores in each patient based on the log2 transformation of the FPKM matrix. Then, surv_cutpoint was used to group patients based on the signature scores. Then, the Surv function and survfit function were used to fit reasonable survival curves, and the ggsurvplot function was used for visualization. The statistical power of survival differences was calculated by chi-squared distribution.

To test the correlation between splicing events and patient survival, the PSI value from TCGA SpliceSeq was used to divide patients into two groups. In the same way mentioned above, the survival curve and statistical power were visualized and quantified.

Identification of malignant cells with CNV estimation

To distinguish malignant cells from nonmalignant cells in human AML generated by 10x Genomics, we inferred large-scale chromosomal copy-number variations (CNVs) in each single cell based on moving averaged expression profiles across chromosomal intervals(20).

RNA Velocity Analysis

Velocyto was used to quantify the unspliced and spliced counts from BAM files in scRNA data by Run10x mode. The mm10 genome and mm10_rmsk references were downloaded from UCSC and used for alignment. RNA velocity analysis was quantified by velocyto. R program.

Construction of a single-cell transcriptomic map of stepwise leukemogenesis in a murine AML model

To dissect the whole process of leukemogenesis, we performed a series of scRNA-seq analyses of bone marrow (BM) cells from mice transplanted with leukemogenic HSPCs at different stages (Figure S1A). Since Myc is frequently amplified and/or overexpressed in human AML (12%), we constructed a Myc-induced murine AML model(21). All the recipient mice developed full-blown AML approximately 8 weeks after transplantation with minimal variation (Figure S1B). The leukemic cells could be identified by the expression of GFP, which was linked with Myc (Figure S1A). A majority (80.3%) of the bone marrow in all transplanted mice was leukemia cells that expressed markers of myeloid lineage (B220-; CD3-; Mac-1+; c-Kit+), and over time, the proportion of myeloid cells increased. (Figure S1C-D). Peripheral blood smears showed step-by-step leukocytosis with increased blasts (Figure S1E). The histological analysis showed that the resulting mouse AML recapitulated the features of human disease (Figure S1F). During leukemogenesis, we harvested the initiating Myc-overexpressing HSPCs (pretransplanted, T₀), preleukemic BM cells at two weeks (T₁), four weeks (T₂) after transplantation and full-blown AML BM cells (T₃, eight weeks after transplantation) for 10x Genomics scRNA-seq analyses, and a total of 41,078 cells were retained for analysis. Two independent repeats for each time point were highly similar to each other (Figure S2A).

Accordingly, we calculated the UMAP plot of all these preleukemic and leukemic cells (PLCs) together with normal bone marrow cells (Figure 1A). The normal cell subpopulations were clustered by their distinct gene signatures, and the remaining GFP-positive cells were defined as the PLC subtype (Figure S2B and Supplementary Table 1). While all the normal cell lineages from each time point largely overlapped, PLCs displayed a shift from being close to neutrophils to being close to HSPCs on the UMAP (Figure 1B). The ratios of normal populations gradually decreased, while those of PLCs increased over time (Figure S2C).

Since it has been proposed that uncontrolled proliferation and blocked differentiation are two major mechanisms of leukemogenesis(7), we analyzed the cell cycle and differentiation status of PLCs. The results showed a gradual increase in the cell cycle for PLCs from T₀ to T₃(Figure 1C and Supplementary Table 2). PLCs also acquired progressively increased differentiation blocks, as indicated by reduced differentiation scores and enhanced stemness scores from T₀ to T₃(Figure 1D and Supplementary Table 2). We also found that the number of UMIs and genes progressively increased in GFP-positive cells, while those of GFP-negative cells did not, which was consistent with a previous report that RNA abundance is associated with elevated expression of stem genes(22) (Figure 1E, S2D-E). Overall, we delineated a single-cell landscape of leukemogenesis that recapitulated the stepwise transition of this disease.

Clinically, more than 30% of MDS cases progress to AML. To visualize the molecular switches in PLCs during leukemogenesis, we performed correlation analysis of the transcriptome between PLC populations at different time points and patients with high-risk MDS, MDS-AML and de novo AML. Of note, the gene signature of T₁ PLCs was highly expressed in patients with high-risk MDS, that of T₂ PLCs was highly expressed in AML patients who developed from MDS and that of T₃ PLCs was highly expressed in de novo AML (Figure 1F and Supplementary Table 3). The correlation of the stepwise molecular switches of PLCs with those of patients suggested that the molecular trajectory of PLC transformation in mice represented a potential mechanism of leukemogenesis from MDS to AML patients.

Progressively deteriorated RNA splicing abnormality during leukemogenesis

To dissect the molecular events underlying stepwise leukemogenesis, we performed gene set variation analysis (GSVA) for specific pathways of PLCs of each stage. We found that the Myc target genes gradually increased from T₀ to T₃, along with genes involved in oxidative phosphorylation and DNA repair, and the genes of the p53 pathways gradually decreased over time, all of which were consistent with the gradually increased aggressiveness of PLCs over time (Figure 2A).

Since myc was the sole driver in this AML model, we first checked the expression levels of Myc. We found that the expression of either ectopic or endogenous Myc was not constitutively increased from T₀ to T₃ (Figure S3A). However, we observed that the expression levels of Myc targets constitutively increased from T₀ to T₃(Figure S3B). Since the expression level of Myc itself was not progressively changed in PLCs, it was unlikely that the increase in Myc targets was just a selection of cells with high levels of Myc, which further suggested that stepwise leukemogenesis was not a result of simply selecting Myc expression. Importantly, by analyzing the transcriptome of the TCGA AML cohort, we found that the signature of Myc target, but not the expression of Myc itself, was associated with the poor prognosis of AML patients (Figure S3C).

Furthermore, the KEGG analysis of the Myc target genes revealed that the spliceosome pathway was the most enriched among the Myc targets in PLCs, along with the cell cycle pathway and RNA transport pathway(23) (Figure 2B). The expression level of splicing factors gradually increased from T₀ to T₃ (Figure 2C and Supplementary Table 4), and stage-specific splicing regulatory factors (SFs) were identified by gene expression (Figure 2D). We found that the expression of these stage-specific SFs associated with PLCs, but not those associated with normal cells, had prognostic value in human AML patients, and the prognostic value of PLC SFs progressively increased from T₁ to T₃ (Figure 2E and Supplementary Table 5). Taken together, these data suggested that aggressively increased expression of splicing factors, independent of the Myc expression level, might underlie the progression from preleukemic to fully transformed leukemia in our stepwise model of leukemogenesis.

Abnormal RNA splicing kinetics distinguish leukemogenic initiating cells from their normal counterparts at the tipping point

Diverging from the normal development route is one of the critical steps at the initiation of tumorigenesis(24). In the landscape of leukemogenesis, we observed a diverging point where HSPCs would differentiate into GMP cells or be transformed into GMP-like preleukemic cells (Figure 3A). The differentiation route from HSPCs to GMP cells was characterized by reduced expression of leukemia signature genes and increased expression of myeloid/neutrophil lineage genes, while in sharp contrast, the transformation route from HSPCs to GMP-like cells displayed increased expression of leukemia signature genes and decreased expression of differentiation genes (Figure 3B). The kinetics of the expression of GMP-specific and GMP-like-specific genes over pseudotime indicated that there was a fate-determining tipping point for Myc-expressing HSPCs to either the GMP differentiation route or GMP-like leukemogenesis route (Figure 3C).

Given the importance of RNA splicing during the whole leukemogenesis process, we analyzed RNA splicing kinetics at the tipping point of leukemogenesis by RNA velocity, which has been suggested to indicate developmental potential and cellular dynamics(25). Consistent with the differentiation tendency of HSPCs into GMPs, HSPCs displayed relatively stationary RNA velocity, and GMPs had high velocity. There was a strong directional flow from HSPCs to GMP. Strikingly, the RNA velocity of GMP-like cells was more stationary than that of HSPCs, and more importantly, the directional flow was from GMP-like cells to HSPCs, which indicated that GMP-like cells at the tipping point may have more potential than HSPCs (Figure 3D). Importantly, we found that at the onset of leukemogenesis, GMP-like cells did not gain a significantly increased AML or leukemic stem cell (LSC) signature (Figure 3E) compared to HSPCs. However, we observed a significant decrease in the expression of RNA splicing factors from HSPCs to GMPs but a significant increase from HSPCs to GMP-like cells (Figure 3F). These active splicing factors might give preleukemic cells a selection pool for leukemia-promoting molecules.

NPM1 at the tipping point is involved in leukemogenesis

To identify the key signatures leading to the binary cell fate choice at the tipping point, we analyzed the differential expression of genes on the two evolutionary routes and then defined the intersection between the low-expressing gene set 2 in the GMP lineage and the high-expressing gene set 4 in the GMP-like lineage as tipping point signatures (TPS), with a total of 30 genes. (Figure S4A-B). We characterized the TPS pattern in the GMP lineage and GMP-like lineage and found that the KEGG_RIBOSOME pathway associated with tumorigenesis also conformed to this pattern (Figure S4C-D).

Furthermore, CRISPR/Cas9 gene editing was used to verify and screen out genes that affect the differentiation and stemness of Myc-overexpressing HSPCs at the tipping point of malignant transformation. To identify specific genes promoting the differentiation of HSPCs, we ranked the sgRNAs by differences in the mean fluorescence intensity (MFI) of Mac-1 or c-Kit staining and the proportion of Mac-1+ or c-Kit+ cells. Based on all detection indicators, Npm1 and Phgdh were finally selected as key molecules for inducing HSPCs malignant transformation during leukemogenesis (Figure S4E-F). In addition, the characteristic molecular interaction network analysis showed that Npm1 was a key node gene in the GMP-like cells, as both the differentially activated module gene and differentially expressed gene in the GMP-like cells suggested that Npm1 might be a key player in leukemogenesis (Figure 4A). As expected, Npm1 mRNA was gradually upregulated in the transformation route but gradually downregulated in the differentiation route, consistent with our TPS pattern (Figure 4B). The abnormal expression of NPM1 has been reported to be involved in human leukemogenesis (26), and consistent with the upregulation of Npm1 expression in GMP-like cells, NPM1 was also significantly upregulated in AML patients compared to normal samples. (Figure 4C). Therefore, we performed further functional studies of Npm1 to verify its effect on the development of leukemia.

To validate the function of Npm1 in leukemogenesis, we disrupted it in Myc-overexpressing HSPCs by CRISPR/Cas9 (Figure S5A-B) and found that Npm1 loss significantly extended the latency of Myc-driven AML in mice (Figure 4D). Consistently, Npm1 deficiency promoted the differentiation of HSPCs in vitro, as indicated by significantly increased expression of the differentiation marker Mac-1 and decreased levels of the stem cell marker c-Kit (Figure 4E). Transcriptomic analysis by RNA-seq also showed that Myc-overexpressing HSPCs promoted cell differentiation while decreasing stemness, with more GMP gene signatures and fewer GMP-like signatures after the deletion of Npm1. (Figure 4F-G and Supplementary Table 6). KEGG also showed that cell differentiation and apoptosis pathways were enriched in Npm1-deficient tumor-initiating cells (Figure 4H). GSEA indicated that the MYC_TARGETS_V1 and MRNA_SPLICING upregulated genes were significantly negatively enriched in HSPCs with sgNpm1 compared to those with sgScramble, which corroborates our previous analysis of global molecular regulatory events behind the evolution of leukemia from T0 to T3(Figure 4I). Consistent with our data, the GMP gene signature was significantly positively enriched in NPM1-deficient human leukemic cells compared to control shREN cells, and the GMP-like gene signature was significantly negatively enriched, suggesting that NPM1 might be a driver of leukemia (Figure 4J).

Phgdh, although rarely reported in AML, was another top gene in the tipping point of GMP-like cells. Its expression was upregulated along the transformation route but downregulated in the differentiation route (Figure S5C). Furthermore, Phgdh disruption promoted the differentiation of HSPCs (Figure S5D). We also observed that Myc-overexpressing HSPCs reduced the GMP-like gene signature, cell stemness and the cell cycle after Phgdh deletion (Figure S5E-F and Supplementary Table 6). KEGG also showed that lineage differentiation pathways were significantly enriched in Phgdh-deficient preleukemia cells, and GSEA showed that the MYC_TARGETS_V1 and MRNA_SPLICING upregulated genes were also significantly negatively enriched in HSPCs with sgPhgdh (Figure S5G-H). High expression of PHGDH was associated with poor prognosis in AML patients (Figure S5I). As outlined above, our data suggested that the AML tipping point signatures, including Npm1 and Phdgh, might be a critical driver for leukemia initiation at the tipping point.

RNA splicing creates explosive heterogeneity of leukemic cells

Once the PLCs went through the tipping point and embarked on the tumorigenesis route, they gradually gained malignancy and heterogeneity, which could be visualized on the URD map (Figure S6A). All these PLCs could be grouped into five major subpopulations, which were named after their top marker genes (Figure 5A-B and Supplementary Table 7). On the map were all the T0 cells at the “root”, followed by T1 and T2 cells, sequentially, at the “stem”. T3 cells were located at the crown and ended with multiple “branches” (Figure 5C). The expression of the human AML signature genes gradually increased along with the leukemogenesis “tree” (Figure S6B). The ratios of these subpopulations significantly varied among the time points, and Hsp90aa1^hi cells did not emerge until T3, suggesting that this subpopulation might be more aggressive than others (Figure S6C). As expected, we found that these subpopulations displayed distinguishable gene signatures in the TCGA_LAML and TARGET_AML databases, and the hazard ratio calculated that the HSP90AA1^hi subtype was significantly associated with high risk in human AML patients (Figure S6D-G).

We further explored the potential molecular mechanisms underlying the heterogeneity, especially those promoting aggressiveness in the Hsp90aa1^hi subpopulation, and found that Hsp90aa1^hi cells expressed significantly higher levels of splicing regulatory factors than any other subpopulation (Figure 5D), suggesting that the splicing machinery generates leukemia intraheterogeneity and promotes its aggressiveness.

To validate the intraheterogeneity in human AML, we analyzed the single-cell transcriptomes of bone marrow cells from five AML patients. The leukemic cells were recognized by referred copy number variations (Figure S7A). These AML cells could be divided into four distinct subpopulations (Figure 5E), and the ratios of these subpopulations significantly varied among these 5 patients (Figure S7B). Of note, the HSP90AA1^hi subpopulation shared the same marker genes, HSP90AA1 and PDLIM1, with the Hsp90aa1^hi subpopulation in mouse AML. Similarly, the TYROBP^hi subpopulation expressed the same top marker genes, TYROBP and SAT1, as the Tyrobp^hisubpopulation in mouse AML (Figure 5F). Importantly, through unbiased hierarchy clustering according to their transcriptomes, we found that the human HSP90AA1^hi subpopulation was most similar to the mouse Hsp90aa1^hi subpopulation, the human TYROBP^hi subpopulation was most similar to the mouse Tyrobp^hisubpopulation, and the human MPO^hi and SPINK2^hi subpopulations were close to the mouse Hes^hi subpopulation (Figure 5G). Furthermore, the HSP90AA1^hi subpopulation was the most aggressive among all cells, as indicated by the highest cell cycle score, while the TYROBP^hi subpopulation had the slowest cell cycle in human AML (Figure 5H). Similar to those in mouse AML, the HSP90AA1^hi cells expressed the highest levels of splicing factors, while the TYROBP^hicells expressed the lowest (Figure 5I). KEGG analysis showed that splicing pathways were also significantly enriched in the HSP90AA1^hi cell population (Figure S7C). The similarity of the intraheterogeneity and its associated molecular features between human and mouse AML suggested that our mouse AML model faithfully recapitulated the molecular and cellular characteristics of human disease and that abnormal RNA splicing generally promoted leukemic heterogeneity in various AML.

Exon 6 skipping of Tmem134, predominantly in the Hsp90aa1^hi subpopulation, promotes aggressiveness in AML

To deeply characterize abnormal RNA splicing events underlying leukemia heterogeneity, we performed single-cell RNA-seq analysis of 178 leukemic cells in T3 by Smart-seq2 instead of the 10x Genomics platform, which would be able to analyze the full length and more transcripts (27, 28). Based on their gene expression patterns, these cells were partitioned into 5 populations, and the tSNE plot showed that the Hsp90aa1^hi subpopulation was separate from others (Figure 5J and S8A). Indeed, the high expression of the signature genes of the Hsp90aa1^hi cell population was also significantly associated with poor prognosis in human AML patients (Figure 5K).

Consistent with those analyzed by 10x Genomics, the signature genes of Hsp90aa1^hi cells were enriched in the RNA splicing pathway and DNA metabolic process pathway (Figure S8B). Furthermore, we found that all five major subtype events of alternative RNA splicing, including exon skipping, intron retention, alternative 5’ splice site, alternative 3’ splice site, and mutually exclusive exons, were increased in Hsp90aa1^hi cells compared with Hsp90aa1^lo cells (Figure S8C). These data validated the RNA splicing in the Hsp90aa1^hi subpopulation, which displayed more aggressiveness.

Among the alternatively spliced genes in Hsp90aa1^hi cells was transmembrane protein 134 (Tmem134), a putative transmembrane protein with two transmembrane domains on the C-terminus(29). The full-length Tmem134, named Tmem134α, and exon 6 skipped Tmem134, named Tmem134β (Figure 6A). TMEM134 is highly conserved in mice and humans, and importantly, the same exon skipping also occurred in human AML (Figure 6B and S8D). We found that exon 6 of Tmem134 was skipped more in Hsp90aa1^hicells than in other cells (Figure 6C). Importantly, exon 6 skipping of TEME134 (TMEM134 PSI low) was significantly associated with poor prognosis in human AML (Figure 6D). In addition, HSPCs with TMEM134α overexpression hindered leukemogenesis in our Myc-induced leukemia murine model compared to the empty vector, while TMEM134β significantly promoted leukemogenesis (Figure 6E and S8E). We observed that compared to TMEM134α-overexpressing tumor-initiating cells expanded rapidly in vivo with a strong growth advantage (Figure S8F). Similarly, TMEM134β promoted Myc-overexpressing HSPC growth, while TMEM134α inhibited HSPC growth ex vivo (Figure 6F).

Furthermore, we tested the potential effects of exon 6 skipping of TMEM134 in Myc-induced AML cells ex vivo. To exclude the influence of endogenous TMEM134 protein, we generated a Tmem134 knockout Myc AML monoclonal cell line (i.e., M#8) with CRISPR/Cas9 and then introduced TMEM134α or TMEM134β cDNA into M8# by retrovirus (Figure S8G-I). The results showed that without endogenous TMEM134, exogenous TMEM134α significantly reduced the growth of leukemic cells compared to the empty vector, while TMEM134β significantly promoted the growth of leukemic cells (Figure 6G). In addition, we used two different strategies by abrogating the splicing site of exon 6 (sgExon6) or deleting the whole exon 6 (sgIntron5 combined with sgIntron6) for enforced exon 6 skipping of Tmem134 in leukemia cells. Both sgExon6 and sgIntron5+6 significantly increased the growth of AML cells (Figure 6H and S8J). Furthermore, we introduced synonymous mutated TMEM134α/β cDNA into the THP-1-cell line and HL-60 cell line (two kinds of human leukemia cell lines) with endogenous TMEM134 knockout. Consistently, TMEM134α significantly repressed the growth of human leukemia cells, and TMEM134β was promoted (Figure 6I and S8K-N). Then, we tested the in vivo effects of TMEM134 exon 6 skipping on leukemia maintenance by transplanting leukemic cells harboring TMEM134α or TMEM134β into recipient mice. We found that recipients with TMEM134α had significantly extended survival (Figure S8O). Despite the short period of tumor development, we have not yet had time to observe the in vivo growth advantage of TMEM134β-overexpressing leukemia cells after the second transplant. Thus, exon 6 skipping of TMEM134, which likely abrogated the tumor suppression function of full-length TMEM134, drove aggressiveness in AML.

Consistently, transcriptomic analysis by RNA-seq showed that TMEM134α significantly repressed the expression of cell cycle genes in AML (Figure S8P). EdU incorporation assays also showed that the proliferation ratios of cells with sgExon6 or sgIntron5+6 were significantly higher than those of the control cells (Figure S8Q). In both human and mouse AML, TMEM134α was significantly associated with a decreased cell cycle (Figure S8R). The genes upregulated by TMEM134α were expressed at significantly lower levels in Hsp90aa1^hicells than in Hsp90aa1^lo cells, and in contrast, the genes downregulated by TMEM134α were expressed at significantly higher levels in Hsp90aa1^hicells than in Hsp90aa1^lo cells (Figure 6J and Supplementary Table 8). GSEA indicated that the NUP98-HOXA9 upregulated genes were significantly negatively enriched in leukemia cells with TMEM134α compared to those with empty vector (Figure 6K). These data strongly indicated that exon 6 skipping of TMEM134 was critical for generating leukemia intraheterogeneity by promoting the aggressiveness of a subpopulation.

Despite the long latency of tumorigenesis and potential opportunities for early diagnosis and treatment, our understanding of its underlying molecular mechanisms is very limited. In this study, we dissected the trajectory of leukemia initiation and progression at single-cell resolution. The molecular features of this trajectory partially recapitulate those in patients. Stepwise tumorigenesis is characterized by gradual deterioration of differentiation block and uncontrolled proliferation (7, 30). At the very beginning, the cells of origin of cancer, HSPCs for AML, with the driver mutations, would go through a tipping point, a critical time for an irreversible choice between “normal” differentiation and malignant transformation. Furthermore, we found that heterogenicity, a hallmark of cancer, is largely acquired at the late stage of tumorigenesis. Our study describes the landscape of the molecular and cellular features of leukemogenesis driven by Myc overexpression. It would be interesting to test whether these features, including the tipping point and the explosive heterogenicity, would be applied to AML driven by other mutations and other cancers driven by MYC (31).

During the tumorigenesis trajectory, an unexpected observation is gradually increased RNA splicing, although the expression level of Myc, the initial driver, remains the same. Many of these abnormally regulated RNA splicing factors are shared in both mouse and human AML. The stage-specific splicing signatures display increased prognostic value over time. MYC can directly regulate the expression of many splicing genes (23, 32). RNA splicing is frequently involved in normal hematopoiesis and hematopoietic malignancies (33–35). It would be important to explore how these splicing genes are stepwise upregulated along the leukemogenesis trajectory without changes in the Myc level in future studies.

Comprehensively profiling abnormal RNA splicing events and investigating the functions of individual abnormally spliced genes in tumorigenesis are essential to better understand the roles of RNA splicing in the process. We show that exon 6 skipping of TMEM134, which is conserved in human and mouse AML, is significantly associated with poor prognosis. Furthermore, we demonstrate that TMEM134α and TMEM134β have opposite functions in malignancy. Interestingly, exon 6 encodes one of the two transmembrane domains of TMEM134 (29). Thus, this AML-associated exon skipping might disrupt the topology of the resulting protein. More work is needed to elucidate the mechanism of TMEM134 exon 6 skipping for its function in AML.

Ethical Approval

All mouse experiments were performed in compliance with the Guide for the Care and Use of Laboratory Animals of Sichuan University and were approved by the Animal Care and Use Committee of Sichuan University.

Competing interests

The authors declare that they have no competing interests.

Author’s contributions

B. Wu, X. Chen and X. Pan designed and performed the experiments, analyzed the data and wrote the manuscript. X. Deng, S. Li, Z. Wang, J. Wang, X. Hang, D. Liao, J. Xu, and M. Chen performed the experiments and analyzed the data. C. Zhao, Z. Xue, Y. Wang, T. Niu, J. Lin provided resources and designed the experiments. L. Chen, C. Chen and Y. Liu conceived the project, designed the experiments, analyzed the data and wrote the manuscript.

Funding

We thank Dr. Yuquan Wei for his generous support and advice, the members of the CC-LY laboratory for their technical support and suggestions, and the Core Facilities of West China Hospital. This work was supported by the National Key R&D Program of China (2017YFA0505600), the National Natural Science Foundation of China (81770157, 81522003, 81670182, 81570150), and the Sichuan Science and Technology Program (2020YFQ0059, 2020ZYD002, 2018JZ0077, 2017TJPT0005).

Availability of data and materials

All scRNA-seq and RNA-seq data in this study are deposited in NCBI Gene Expression Omnibus under accession number GSE142645. The private code of GSE142645 is sbcpaygmdrghjav.

The analysis code can be found at GitHub (https://github.com/pangxueyu233/Super-RNA-splicing-entropy-drives-stepwise-leukemogenesis).

E. R. Fearon, B. Vogelstein, A genetic model for colorectal tumorigenesis. Cell 61, 759–767 (1990).
C. Lengauer, K. W. Kinzler, B. Vogelstein, Genetic instabilities in human cancers. Nature 396, 643–649 (1998).
S. D. Markowitz, M. M. Bertagnolli, Molecular origins of cancer: Molecular basis of colorectal cancer. N Engl J Med 361, 2449–2460 (2009).
P. C. Nowell, The clonal evolution of tumor cell populations. Science (New York, N.Y.) 194, 23–28 (1976).
M. Greaves, C. C. Maley, Clonal evolution in cancer. Nature 481, 306–313 (2012).
H. Dohner, D. J. Weisdorf, C. D. Bloomfield, Acute Myeloid Leukemia. N Engl J Med 373, 1136–1152 (2015).
L. M. Kelly, D. G. Gilliland, Genetics of myeloid leukemias. Annu Rev Genomics Hum Genet 3, 179–198 (2002).
C. Kandoth et al., Mutational landscape and significance across 12 major cancer types. Nature 502, 333–339 (2013).
A. G. Bick et al., Inherited causes of clonal haematopoiesis in 97,691 whole genomes. Nature 586, 763–768 (2020).
E. L. Bao et al., Inherited myeloproliferative neoplasm risk affects haematopoietic stem cells. Nature 586, 769–775 (2020).
C. Terao et al., Chromosomal alterations among age-related haematopoietic clones in Japan. Nature 584, 130–135 (2020).
P. R. Loh et al., Insights into clonal haematopoiesis from 8,342 mosaic chromosomal alterations. Nature 559, 350–355 (2018).
Y. Liu et al., Deletions linked to TP53 loss drive cancer through p53-independent mechanisms. Nature 531, 471–475 (2016).
S. Picelli et al., Full-length RNA-seq from single cells using Smart-seq2. Nature protocols 9, 171–181 (2014).
X. Han et al., Mapping the Mouse Cell Atlas by Microwell-Seq. Cell 173, 1307 (2018).
Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature 562, 367–372 (2018).
S. Nestorowa et al., A single-cell resolution map of mouse hematopoietic stem and progenitor cell differentiation. Blood 128, e20-31 (2016).
A. G. Kotini et al., Stage-Specific Human Induced Pluripotent Stem Cells Map the Progression of Myeloid Transformation to Transplantable Leukemia. Cell stem cell 20, 315–328.e317 (2017).
P. J. Valk et al., Prognostically useful gene-expression profiles in acute myeloid leukemia. The New England journal of medicine 350, 1617–1628 (2004).
A. P. Patel et al., Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science (New York, N.Y.) 344, 1396–1401 (2014).
H. Luo et al., c-Myc rapidly induces acute myeloid leukemia in mice without evidence of lymphoma-associated antiapoptotic mutations. Blood 106, 2452–2461 (2005).
G. S. Gulati et al., Single-cell transcriptional diversity is a hallmark of developmental potential. Science (New York, N.Y.) 367, 405–411 (2020).
C. M. Koh et al., MYC regulates the core pre-mRNA splicing machinery as an essential step in lymphomagenesis. Nature 523, 96–100 (2015).
X. Wang et al., Sequential fate-switches in stem-like cells drive the tumorigenic trajectory from human neural stem cells to malignant glioma. Cell Res 31, 684–702 (2021).
G. La Manno et al., RNA velocity of single cells. Nature 560, 494–498 (2018).
C. Thiede et al., Prevalence and prognostic impact of NPM1 mutations in 1485 adult patients with acute myeloid leukemia (AML). Blood 107, 4011–4020 (2006).
S. Picelli et al., Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nature methods 10, 1096–1098 (2013).
X. Wang, Y. He, Q. Zhang, X. Ren, Z. Zhang, Direct Comparative Analyses of 10X Genomics Chromium and Smart-seq2. Genomics Proteomics Bioinformatics, (2021).
P. Talaty, A. Emery, K. Holthusen, D. N. Everly, Jr., Identification of transmembrane protein 134 as a novel LMP1-binding protein by using bimolecular fluorescence complementation and an enhanced retroviral mutagen. J Virol 86, 11345–11355 (2012).
R. Pardal, M. F. Clarke, S. J. Morrison, Applying the principles of stem-cell biology to cancer. Nat Rev Cancer 3, 895–902 (2003).
Y. Liu et al., Convergence of oncogenic cooperation at single-cell and single-gene levels drives leukemic transformation. Nature communications 12, 6323 (2021).
C. J. David, M. Chen, M. Assanah, P. Canoll, J. L. Manley, HnRNP proteins controlled by c-Myc deregulate pyruvate kinase mRNA splicing in cancer. Nature 463, 364–368 (2010).
L. Chen et al., Genetic Drivers of Epigenetic and Transcriptional Variation in Human Immune Cells. Cell 167, 1398–1414 e1324 (2016).
L. Chen et al., Transcriptional diversity during lineage commitment of human blood progenitors. Science 345, 1251033 (2014).
F. Wang et al., Single-cell architecture and functional requirement of alternative splicing during hematopoietic stem cell formation. Sci Adv 8, eabg5369 (2022).

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Dissecting the tumorigenesis trajectory by single-cell RNA-seq reveals critical roles of RNA splicing

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Background

Methods

Mice

Murine AML model and AML cell line

Cas9 AML cell line

Retroviral constructs

CRISPRs

Cell growth curve

Cell proliferation and cell differentiation

RNA-seq analyses

SMART-seq2 library preparation and sequencing

Single-cell collection, library preparation and sequencing

Single-cell gene expression quantification

Batch effect reduction, dimension reduction and unsupervised clustering

Trajectory construction and pseudotime calculation

GSVA, gene regulatory network analysis and molecular interaction network analysis

Generation of gene signatures and visualizations

Splicing events calculation

Cell cycle score calculation

Bulk data processing and signature identification in MDS and AML patients

TCGA data analysis

Identification of malignant cells with CNV estimation

RNA Velocity Analysis

Results

Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1