Single-cell transcriptome analyses reveal critical roles of RNA splicing during leukemia progression

JournalFeeds

Authors:
Baohong Wu , et al.

Citation: Wu B, Chen X, Pan X, Deng X, Li S, Wang Z, et al. (2023) Single-cell transcriptome analyses reveal critical roles of RNA splicing during leukemia progression. PLoS Biol 21(5):
e3002088.

https://doi.org/10.1371/journal.pbio.3002088

Academic Editor: Connie J. Eaves, B.C. Cancer Agency, CANADA

Received: October 5, 2022; Accepted: March 21, 2023; Published: May 2, 2023

Copyright: © 2023 Wu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All scRNA-seq and RNA-seq data in this study are deposited in NCBI Gene Expression Omnibus under accession number GSE142645. The analysis code can be found at GitHub (https://github.com/pangxueyu233/Super-RNA-splicing-entropy-drives-stepwise-leukemogenesis). All flow cytometry data for S1 Fig are deposited in the FlowRepository database under accession number FR-FCM-Z64U and FR-FCM-Z64T.

Funding: This work was supported by the National Natural Science Foundation of China (82130007 to YL, 82170171 to CC), the Sichuan Science and Technology Program (2021YFS0027 to LC, 2020YFQ0059 to CC, 2018JZ0077 to YL), the 1.3.5. Project for Disciplines of Excellence, West China Hospital, Sichuan University (ZYJC21009 to YL and ZYGD22012 to YL), the Sichuan University Postdoctoral Interdisciplinary Innovation Fund to XC, the Post-Doctor Research Project, West China Hospital, Sichuan University (2023HXBH019 to BW and 2023HXBH033 to XC), the China Postdoctoral Science Foundation (2022M722272 to JX). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Abbreviations:
AML,
acute myeloid leukemia; BM,
bone marrow; CH,
clonal hematopoiesis; CNV,
copy number variation; DEG,
differentially expressed gene; GSEA,
gene set enrichment analysis; GSVA,
gene set variation analysis; HSPC,
hematopoietic stem and progenitor cell; LSC,
leukemic stem cell; MDS,
myelodysplastic syndrome; MFI,
mean fluorescence intensity; MPN,
myeloproliferative neoplasm; PLC,
preleukemic and leukemic cell; PSI,
percent-splice-in; scRNA-seq,
single-cell RNA sequencing; SF,
splicing factor; TPS,
tipping point signature; UMI,
unique molecular identifier; WBM,
whole bone marrow

Introduction

Tumorigenesis, a process to accumulate enough genetic alterations for normal cells to be transformed into malignant cells, can take decades in patients [15]. 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 and MDS [912]. 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.

Results

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

Since MYC is frequently amplified and/or overexpressed in human AML (13%) and in MDS which may develop AML later and was previously called as preleukemia, we constructed an Myc-induced murine AML model (S1A–S1C Fig) [13]. 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 (S1D Fig). All the recipient mice developed full-blown AML approximately 8 weeks after transplantation with minimal variation (S1E Fig). The leukemic cells could be identified by the expression of GFP, which was linked with Myc. 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 (S1F and S1G Fig). Peripheral blood smears showed step-by-step leukocytosis with increased blasts (S1H Fig). The histological analysis showed that the resulting mouse AML recapitulated the features of human disease (S1I Fig). During leukemogenesis, we harvested the initiating Myc-overexpressing HSPCs (pretransplanted, T0), preleukemic BM cells at 2 weeks (T1), 4 weeks (T2) after transplantation, and full-blown AML BM cells (T3, 8 weeks after transplantation) for 10× 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 (S2A Fig).

Accordingly, we calculated the UMAP plot of all these preleukemic and leukemic cells (PLCs) together with normal BM cells (Fig 1A). The normal cell subpopulations were clustered by their distinct gene signatures, and the remaining GFP-positive cells were defined as the PLC subtype (S2B Fig). 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 (Fig 1B). The ratios of normal populations gradually decreased, while those of PLCs increased over time (S2C Fig).

thumbnail

Fig 1. Serial scRNA-seq analyses revealing the landscape of progressive leukemogenesis.

(A) UMAP plot of normal cells and PLCs for all 9 samples at 5 time points (normal whole bone marrow (2×), T0 (before transplant, 1×), T1 (2 weeks after transplant, 2×), T2 (4 weeks after transplant, 2×), and T3 (8 weeks after transplant, 2×)) of leukemogenesis. Cells are colored according to Myc-GFP expression (upper left corner). (B) UMAP plot of normal cells and PLCs during leukemogenesis ordered by time point. (C) Density scatter plots showing G2 M scores and S scores in PLCs during leukemogenesis. (D) Density scatter plots showing differentiation and stemness signature scores in PLCs during leukemogenesis. (E) Bar graph showing the number of UMI counts of Myc-GFP-positive and Myc-GFP-negative cells at 3 time points; p values were calculated by the Wilcoxon signed-rank test. (F) The heatmap displaying the normalized expression levels of high-risk MDS, MDS_AML, and AML signatures in PLCs at each time point. The underlying data for Fig 1C, 1D, 1E and 1F can be found in S1 Data. AML, acute myeloid leukemia; MDS, myelodysplastic syndrome; PLC, preleukemic and leukemic cell; scRNA-seq, single-cell RNA sequencing; UMI, unique molecular identifier.


https://doi.org/10.1371/journal.pbio.3002088.g001

Since it has been proposed that uncontrolled proliferation and blocked differentiation are 2 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 T0 to T3 (Fig 1C). PLCs also acquired progressively increased differentiation blocks, as indicated by reduced differentiation scores and enhanced stemness scores from T0 to T3 (Fig 1D). We also found that the number of unique molecular identifiers (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 [14] (Figs 1E and S2D and S2E). 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 [15]. Of note, the gene signature of T1 PLCs was highly expressed in patients with high-risk MDS, that of T2 PLCs was highly expressed in AML patients who developed from MDS and that of T3 PLCs was highly expressed in de novo AML (Fig 1F). 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 T0 to T3, 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 (Fig 2A).

thumbnail

Fig 2. The Myc targets, not Myc, triggered a progressively deteriorating RNA splicing abnormality during leukemogenesis.

(A) Heatmap showing GSVA scores of hallmark pathways at each time point during leukemogenesis. (B) The top 10 KEGG pathways of Myc_targets detached in PLCs. (C) The box plot showed the relative expression levels of SFs in PLC during leukemogenesis; p values were calculated by the Wilcoxon signed-rank test. (D) Heatmap showing the relative expression levels of 4 subtypes of MYC targets involved in high SFs (rows) during leukemogenesis (columns) in murine AML. (E) Kaplan–Meier survival curves of TARGET-AML patients with low or high SFs expression at T0, T1, T2, and T3. The p value was calculated by the log-rank test. The underlying data for Fig 2C and 2E can be found in S1 Data. AML, acute myeloid leukemia; GSVA, gene set variation analysis; PLC, preleukemic and leukemic cell; SF, splicing factor.


https://doi.org/10.1371/journal.pbio.3002088.g002

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 T0 to T3 (S3A Fig). However, we observed that the expression levels of Myc targets constitutively increased from T0 to T3 (S3B Fig). 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 (S3C Fig).

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 [16] (Fig 2B). The expression level of splicing factors gradually increased from T0 to T3 (Fig 2C), and stage-specific splicing regulatory factors (SFs) were identified by gene expression (Figs 2D and S3D). 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 T1 to T3 (Fig 2E). 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 [17]. 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 (Fig 3A). Of note, the progenitor-like and GMP-like AML patients’ signatures are mostly expressing in our GMP-like preleukemic cells [18] (S4A Fig). GFP+ cells contained quite a lot of HSC population at the start T0 that were reduced to almost zero at the end time point T3 (S4B Fig). Consistently, 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 (Fig 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 (Fig 3C).

thumbnail

Fig 3. A tipping point for normal differentiation and malignant transformation during leukemogenesis.

(A) UMAP plot of the tipping point for Myc-GFP-positive HSPCs, GMPs, and GMP-like cells. (B) The UMAP plots showing the expression levels of GMP classical markers, Klf2 and Cd177 (top), and GMP-like markers, Mif and Ncl (bottom). (C) The differentiation trajectory from HSPCs to GMP cells (top) and the transformation trajectory from HSPCs to GMP-like cells (bottom) calculated by Slingshot. (D) The RNA velocity at the tipping point. (E) The box plots showed the expression levels of the AML and LSC signatures in HSPCs, GMPs, and GMP-like cells; p values were calculated by the Wilcoxon signed-rank test. (F) The box plots showed the expression levels of the SF signature in HSPCs, GMPs, and GMP-like cells; p values were calculated by the Wilcoxon signed-rank test. The underlying data for Fig 3C, 3E and 3F can be found in S1 Data. AML, acute myeloid leukemia; HSPC, hematopoietic stem and progenitor cell; LSC, leukemic stem cell; SF, splicing factor.


https://doi.org/10.1371/journal.pbio.3002088.g003

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 [19]. 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 (Fig 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 (Fig 3E) compared to HSPCs. However, we observed a significant decrease in the expression of RNA SFs from HSPCs to GMPs but a significant increase from HSPCs to GMP-like cells (Fig 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 2 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 (S4C and S4D Fig). 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 (S4E and S4F Fig).

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 (S4G and S4H Fig). 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 (DEG) in the GMP-like cells suggested that Npm1 might be a key player in leukemogenesis (Fig 4A). As expected, Npm1 mRNA was gradually up-regulated in the transformation route and malignant progression, but gradually down-regulated in the differentiation route, consistent with our TPS pattern (Figs 4B and S4I). The abnormal expression of NPM1 has been reported to be involved in human leukemogenesis [20], and consistent with the up-regulation of Npm1 expression in GMP-like cells, NPM1 was also significantly up-regulated in AML patients compared to normal samples (Fig 4C). Therefore, we performed further functional studies of Npm1 to verify its effect on the development of leukemia.

thumbnail

Fig 4. Npm1 acted as a tipping point signature to promote AML.

(A) The molecular interaction networks showing the node genes in GMP-like cells. Red represents DAMGs and DEGs, simultaneously; orange represents DEGs; green represents DAMGs; gray represents connecting Steiner genes. (B) The relative expression of Npm1 along the differentiation (left) and transforming (right) trajectories. (C) The box plot showed the expression levels of NPM1 in TCGA-LAML and normal samples. ***p.adj < 0.001; p values were calculated by Wilcoxon signed-rank test. (D) The Kaplan–Meier survival curve of recipient mice transplanted with Myc-GFP HSPCs harboring sgScramble or sgNpm1. **p < 0.01 (log-rank test). (E) The bar graphs showed the MFI of differentiated cells (left) and stem cells (right) at 2 days, 4 days, and 6 days in sgNpm1 samples compared to scramble samples. *p < 0.05, **p < 0.01, *** p < 0.001; P values were calculated by unpaired parametric t test. (F) Scatter plot showing differentiation and stemness signature scores in sgNpm1 and sgScramble samples. Measured by RNA-seq. (G) The box plot showed the expression levels of the GMP-like signature (left) and GMP signature (right) in sgNpm1 and sgScramble samples. *p < 0.05, **p < 0.01, p values were calculated by Chi-squared test. (H) The KEGG pathways enriched in sgScramble samples compared to sgNpm1 samples (top). KEGG pathways enriched in sgNpm1 samples compared to sgScramble samples (bottom). (I) GSEA showing the negative enrichment of HALLMARK_MYC_TARGETS_V1 and HALLMARK_MRNA_SPLICING in sgNpm1 cells compared to sgScramble cells. (J) GSEA showing the negative and positive enrichments of GMP-like and GMP signatures in shNPM1 K562 cells compared to shREN cells. The underlying data for Fig 4C–4J can be found in S1 Data. AML, acute myeloid leukemia; DAMG, differentially activated module gene; DEG, differentially expressed gene; GSEA, gene set enrichment analysis; HSPC, hematopoietic stem and progenitor cell; MFI, mean fluorescence intensity.


https://doi.org/10.1371/journal.pbio.3002088.g004

To validate the function of Npm1 in leukemogenesis, we disrupted it in Myc-overexpressing HSPCs by CRISPR/Cas9 (S5A–S5C Fig) and found that Npm1 loss significantly extended the latency of Myc-driven AML in mice (Fig 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 (Fig 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 (Fig 4F and 4G). KEGG also showed that cell differentiation and apoptosis pathways were enriched in Npm1-deficient tumor-initiating cells (Fig 4H). GSEA indicated that the MYC_TARGETS_V1 and MRNA_SPLICING up-regulated 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 (Fig 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 (Fig 4J).

Phgdh, although rarely reported in AML, was another top gene in the tipping point of GMP-like cells. Its expression was up-regulated along the transformation route but down-regulated in the differentiation route (S5D Fig). Furthermore, Phgdh disruption promoted the differentiation of HSPCs (S5E Fig). We also observed that Myc-overexpressing HSPCs reduced the GMP-like gene signature, cell stemness and the cell cycle after Phgdh deletion (S5F and S5G Fig). 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 up-regulated genes were also significantly negatively enriched in HSPCs with sgPhgdh (S5H and S5I Fig). High expression of PHGDH was associated with poor prognosis in AML patients (S5J Fig). We also have observed up-regulated expression of NPM1 and PHGDH in most malignancies compared to normal samples (S5K and S5L Fig). 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 (S6A Fig). All these PLCs could be grouped into 5 major subpopulations, which were named after their top marker genes (Fig 5A and 5B). 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” (Fig 5C). The expression of the human AML signature genes gradually increased along with the leukemogenesis “tree” (S6B Fig). The ratios of these subpopulations significantly varied among the time points, and Hsp90aa1hi cells did not emerge until T3, suggesting that this subpopulation might be more aggressive than others (S6C Fig). Consistently, distinguished expression modules were identified in PLCs at different stages T0, T1, T2, and T3, suggesting that the molecular features of PLCs are changed during leukemogenesis (S6D Fig). The Hsp90aa1hi population had the highest expression of stemness score and the Npghi population had the highest expression of differentiation score (S6E Fig). Moreover, the 5 clusters represented different differentiation stages: the Hsp90aa1hi and Ifitm3 hi populations were more stemness and the Npghi population was well differentiated (S6F Fig). 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 HSP90AA1hi subtype was significantly associated with high risk in human AML patients and HSP90AA1 expressions were up-regulated in even most malignancies (S6G–S6K Fig). Consistently, we also found a positive correlation between the expression levels of Npm1 or Phgdh gene and the Hsp90aa1hi subpopulation signatures in tipping point cells (S6L Fig).

thumbnail

Fig 5. Explosive intraleukemia heterogeneity during leukemogenesis.

(A) The force-directed layout maps showing leukemogenesis trajectories, colored by cell subtype. (B) The dot plot showed the expression of cell type-specific genes (column) across cell populations in PLCs. (C) The force-directed layout maps showing leukemogenesis, colored by time series. (D) The box plot showed the expression levels of SFs in each subpopulation; p values were calculated by the Wilcoxon signed-rank test. (E) UMAP plot of 165,33 single cells from bone marrow of 5 AML patients, colored by cell population. (F) The dot plot showed the expression of cell type-specific genes (column) across cell populations in AML patients. (G) Hierarchy clustering showing the orthologous murine (mmu) and human (hsa) PLC subtypes at the single-cell level. (H) The box plot showed the expression levels of cell cycle signatures in PLC subtypes in AML patients; p values were calculated by the Wilcoxon signed-rank test. (I) The box plot showed the expression levels of SFs in PLC subtypes in AML patients; p values were calculated by the Wilcoxon signed-rank test. (J) The t-SNE plot showing the Myc-GFP-positive cells in T3, performed by Smart-seq2, colored by cell types. (K) Kaplan–Meier survival curves of TCGA-LAML patients with low or high Hsp90aa1hi cell signatures. The p value was calculated by the log-rank test. The underlying data for Fig 5B, 5D, 5F and 5H–5K can be found in S1 Data. AML, acute myeloid leukemia; PLC, preleukemic and leukemic cell; SF, splicing factor.


https://doi.org/10.1371/journal.pbio.3002088.g005

We further explored the potential molecular mechanisms underlying the heterogeneity, especially those promoting aggressiveness in the Hsp90aa1hi subpopulation, and found that Hsp90aa1hi cells expressed significantly higher levels of splicing factors than any other subpopulation (Fig 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 BM cells from 5 AML patients. The leukemic cells were recognized by referred copy number variations (CNVs) (S7A Fig). These AML cells could be divided into 4 distinct subpopulations (Fig 5E), and the ratios of these subpopulations significantly varied among these 5 patients (S7B Fig). Of note, the HSP90AA1hi subpopulation shared the same marker genes, HSP90AA1 and PDLIM1, with the Hsp90aa1hi subpopulation in mouse AML. Similarly, the TYROBPhi subpopulation expressed the same top marker genes, TYROBP and SAT1, as the Tyrobphi subpopulation in mouse AML (Fig 5F). Importantly, through unbiased hierarchy clustering according to their transcriptomes, we found that the human HSP90AA1hi subpopulation was most similar to the mouse Hsp90aa1hi subpopulation, the human TYROBPhi subpopulation was most similar to the mouse Tyrobphi subpopulation, and the human MPOhi and SPINK2hi subpopulations were close to the mouse Heshi subpopulation (Fig 5G). Furthermore, the HSP90AA1hi subpopulation was the most aggressive among all cells, as indicated by the highest cell cycle score, while the TYROBPhi subpopulation had the slowest cell cycle in human AML (Fig 5H). Similar to those in mouse AML, the HSP90AA1hi cells expressed the highest levels of splicing factors, while the TYROBPhi cells expressed the lowest (Fig 5I). KEGG analysis showed that splicing pathways were also significantly enriched in the HSP90AA1hi cell population (S7C Fig). 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 Hsp90aa1hi 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 10× Genomics platform, which would be able to analyze the full length and more transcripts [21,22]. Based on their gene expression patterns, these cells were partitioned into 5 populations, and the t-SNE plot showed that the Hsp90aa1hi subpopulation was separate from others (Figs 5J and S8A). Indeed, the high expression of the signature genes of the Hsp90aa1hi cell population was also significantly associated with poor prognosis in human AML patients (Fig 5K).

Consistent with those analyzed by 10× Genomics, the signature genes of Hsp90aa1hi cells were enriched in the RNA splicing pathway and DNA metabolic process pathway (S8B Fig). Furthermore, we found that all 5 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 Hsp90aa1hi cells compared with Hsp90aa1lo cells (S8C Fig). These data validated the RNA splicing in the Hsp90aa1hi subpopulation, which displayed more aggressiveness.

Among the alternatively spliced genes in Hsp90aa1hi cells was transmembrane protein 134 (Tmem134), a putative transmembrane protein with 2 transmembrane domains on the C-terminus [23]. The full-length Tmem134, named Tmem134α, and exon 6 skipped Tmem134, named Tmem134β (Fig 6A). TMEM134 is highly conserved in mice and humans, and importantly, the same exon skipping also occurred in human AML (Figs 6B and S8D). Furthermore, we also measured the PSI values of TMEM134 exon 6 skipping in multiple cancers and found that it happens in the most of them more or less, suggesting that it may have the general function in pan-cancer (S8E Fig). We found that exon 6 of Tmem134 was skipped more in Hsp90aa1hi cells than in other cells (Fig 6C). Importantly, exon 6 skipping of TEME134 was also validated in multiple human AML cell lines and 6 AML patients (Fig 6D and 6E). 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 (Figs 6F and S9A). We observed that compared to TMEM134α, TMEM134β-overexpressing tumor-initiating cells proliferated rapidly in vivo with a strong growth advantage (S9B Fig). Similarly, TMEM134β promoted Myc-overexpressing HSPC growth, while TMEM134α inhibited HSPC growth ex vivo (Fig 6G).

thumbnail

Fig 6. The exon 6 skipping of Tmem134, associated with the most aggressive Hsp90aa1hi subpopulation, promoted the cell cycle of leukemic cells.

(A) Schematic diagram showing the frequent occurrence of homologous skipping events of TMEM134 exon 6 in murine and human AML. (B) The percentages of TMEM134 exon 6 skipping occurrences in normal and AML patients. Normal (n = 3), AML patients (n = 8). (C) The bar graph showed the Tmem134 exon 6 skipping ratio in Hsp90aa1lo and Hsp90aa1hi cells. *p < 0.05 (Wilcoxon signed-rank test). (D) Semiquantitative PCR showing of the relative mRNA level of TMEM134α or TMEM134β in human AML cell lines. (E) Semiquantitative PCR showing of the relative mRNA level of TMEM134α or TMEM134β in 6 human AML patients. Overview of AML patients including gender, age, blast percent was shown above and genetic alterations was listed in S3 Table. (F) The survival curves of recipient mice transplanted with Myc-overexpressing c-Kit+ HSPCs overexpressing vector, TMEM134α or TMEM134β. **p < 0.01 (log-rank test). (G) Growth curves of c-Kit+ HSPCs overexpressing vector, TMEM134α or TMEM134β under Myc at the indicated time points. **p < 0.01, ***p < 0.001. P values were calculated by unpaired parametric t tests. (H) The growth curves of endogenous Tmem134 KO mouse AML cells overexpressing vector, TMEM134α or TMEM134β at the indicated time points. *p < 0.05, **p < 0.01. P values were calculated by unpaired parametric t tests. (I) Growth curves of Cas9+ AML cells transduced with sgScramble, sgExon6, or sgIntron6+7 at the indicated time points. **p < 0.01, ***p < 0.001. P values were calculated by unpaired parametric t tests. (J) The growth curves of endogenous TMEM134 KO THP-1 cells overexpressing vector, TMEM134α or TMEM134β at the indicated time points. **p < 0.01, ***p < 0.001, p values were calculated by unpaired parametric t test. (K) The box plot showed the relative expression levels of TMEM134α up- and down-regulated genes compared to the vector in the Hsp90aa1lo and Hsp90aa1hi subpopulations. Each bar represents the mean ± SD. ****p < 0.0001, p values were calculated using an unpaired parametric t test. (L) GSEA showing that the TAKDDA_TARGETS_of_NUP98_HOXA9_FUSION_UP gene set was negatively enriched in TMEM134α AML cells compared to vector AML cells. The underlying data for Fig 6B, 6C and 6F–6L can be found in S1 Data. AML, acute myeloid leukemia; GSEA, gene set enrichment analysis; HSPC, hematopoietic stem and progenitor cell.


https://doi.org/10.1371/journal.pbio.3002088.g006

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 (S9C–S9G Fig). 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 (Fig 6H). In addition, we used 2 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 (S9H Fig). Both sgExon6 and sgIntron5+6 significantly increased the growth of AML cells (Fig 6I). Furthermore, we introduced synonymous mutated TMEM134α/β cDNA into the THP-1 and HL-60 cell line (2 kinds of human leukemia cell lines) with endogenous TMEM134 knockout (S9I–S9K Fig). Consistently, TMEM134α significantly repressed the growth of human leukemia cells, and TMEM134β was promoted (Figs 6J and S9L). 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 (S9M Fig). 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, EdU incorporation assays showed that the proliferation ratios of cells with sgExon6 or sgIntron5+6 were significantly higher than those of the control cells (S9N Fig). In both human and mouse AML, TMEM134α was significantly associated with a decreased cell cycle (S9O Fig). Transcriptomic analysis by RNA-seq also showed that TMEM134α significantly repressed the expression of cell cycle genes in AML (S9P Fig). The genes up-regulated by TMEM134α were expressed at significantly lower levels in Hsp90aa1hi cells than in Hsp90aa1lo cells, and in contrast, the genes down-regulated by TMEM134α were expressed at significantly higher levels in Hsp90aa1hi cells than in Hsp90aa1lo cells (Fig 6K). GSEA indicated that the NUP98-HOXA9 up-regulated genes were significantly negatively enriched in leukemia cells with TMEM134α compared to those with empty vector (Figs 6L and S9Q). These data strongly indicated that exon 6 skipping of TMEM134 was critical for generating leukemia intraheterogeneity by promoting the aggressiveness of a subpopulation.

Discussion

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,24]. 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 [25].

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 [16,26]. RNA splicing is frequently involved in normal hematopoiesis and hematopoietic malignancies [2729]. It would be important to explore how these splicing genes are stepwise up-regulated 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 2 transmembrane domains of TMEM134 [23]. 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.

Methods

Retroviral constructs

In Fig 6F and 6G, the cDNA elements overexpressing TMEM134α and TMEM134β were cloned into pMSCV-IRES-GFP-Myc (pMIG-Myc). In Figs 6H and 6J and S9L and S9M, the cDNA elements overexpressing MYC, TMEM134α, and TMEM134β (synthesized by Sangon Biotec, Shanghai) were cloned into pMSCV-IRES-GFP (pMIG) as previously reported [30] 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 S1 Table. 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 S2 Table.

SMART-seq2 library preparation and sequencing

SMART-seq2 was performed on AML cells at T3, following a modified SMART-seq2 protocol published previously [31]. 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.

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 [34]. 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:

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 [15] and AML patient data GSE1159 [35] 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 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.

Supporting information

S1 Fig. Identification of leukemogenesis landscape.

Related to Fig 1. (A) Copy number alteration of MYC in TCGA AML database. Red, Amplification; Pink, Gain; Gray, No alterations. (B) The plots showing the expression levels of MYC in normal samples and AML patients in TARGET_AML and TCGA AML cohorts. **** p.adj < 0.0001.p values were calculated by Wilcoxon test. Created with BioRender.com. (C) The plots showing the expression levels of MYC in normal samples and MDS patients generated from GSE107400. **** p.adj < 0.0001.p values were calculated by Wilcoxon test. (D) The overview of study design. (E) Kaplan–Meier survival curve of mice with Myc-GFP-induced murine AML model, (n = 6). (F) Representative flow cytometry profile of GFP-positive cells from bone marrow in T3 time point. (G) Flow cytometry profiles showing the percentage of lymphocyte and myeloid cells from peripheral blood in T1, T2, and T3 of leukemogenesis. X-axis is myeloid marker Mac-1; Y-axis is lymphocyte markers B220/CD3. (H) Blood smear of peripheral blood during T1 through T3 in Myc-GFP leukemic mouse. Scale bar: 10 μm. (I) Hematoxylin–eosin images of bone marrow (top), spleen (middle), and liver (bottom) during T1 through T3 in Myc-GFP leukemic mouse. Scale bar: 20 μm. The underlying data for S1B, S1C and S1E Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s001

(TIF)

S2 Fig. The characterization of murine AML during leukemogenesis.

Related to Fig 1. (A) The UMAP plot of single cells for all 9 samples at 5 time points, colored by samples. (B) The heatmap showing the expression levels of cell type-specific genes among all cells, and classic markers were labeled on the right. (C) The line charts showing the dynamics percentage of each cell type during leukemogenesis. (D) Bar graph showing the number of genes in 3 times point; p values were calculated by Wilcoxon signed-rank test. (E) Bar graph showing the RNA concentration of Myc-GFP-positive cells and Myc-GFP-negative cells at T1 time point. ***p < 0.001, p values were calculated using an unpaired parametric t test. The underlying data for S2B–S2E Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s002

(TIF)

S3 Fig. MYC/MYC targets during leukemogenesis and stage-specific SFs in human MDS to AML.

Related to Fig 2. (A) The box plot showing the expression levels of Myc at each time point during leukemogenesis; p values were calculated by Wilcoxon signed-rank test. (B) The box plot showing the expression levels of Myc_target at each time point during leukemogenesis; p values were calculated by Wilcoxon signed-rank test. (C) The Kaplan–Meier survival curves of TARGET-AML patients with low or high MYC/MYC targets; p value was calculated by log-rank test. (D) The heatmap showing the relative expression levels of 4 subtypes of MYC targets involved in splicing factors (columns), during leukemogenesis (rows), in human AML. The underlying data for S3A–S3C Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s003

(TIF)

S4 Fig. Differential expression between GMP and GMP like.

Related to Fig 3. (A) The heatmap showing the expression levels of normal HSC progenitor, GMP, myeloid and AML HSC-like, progenitor-like, GMP-like makers in the Galen and colleagues, Cell 2019 (rows) in tipping point HSPC, GMP, and GMP-like cells (columns) [18]. (B) The UMAP plot of single cells at tipping point stages split by samples and colored by cell annotation. (C) Heatmaps showing gene expression dynamics over pseudo time in GMP lineage (left) and GMP-like (right) lineage. Genes (row) are clustered and cells (column) are ordered according to the pseudo time. Gene set1: up-regulate genes in GMP lineage, Gene set2: down-regulate genes in GMP lineage, Gene set3: down-regulate genes in GMP-like lineage, Gene set4: up-regulate genes in GMP-like lineage. (D) Venn diagram showing the overlap among gene set1, gene set2, gene set3, and gene set4. (E) The tipping point signature expression trends along GMP lineage (left) and GMP-like lineage (right). Single cell colored by cluster annotation. (F) The Ribosome signature expression trends along GMP lineage (left) and GMP-like lineage (right). (G) Relative proportions of Mac-1+ (top) and of c-Kit+ (bottom) cell populations, measured by flow cytometry, in HSPCs infected with CRISPRs targeting tipping point genes or scramble. (H) Relative MFIs of Mac-1 (top) and of c-Kit (bottom) staining in HSPCs infected with CRISPRs targeting tipping point genes or scramble. (I) The box plot showing the expression levels of Npm1/Phgdh in PLCs at each time point during leukemogenesis. ****p.adj < 0.0001.p values were calculated by Wilcoxon test. The underlying data for S4G–S4I Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s004

(TIF)

S5 Fig. Phgdh acted as another tipping point signature.

Related to Fig 4. (A) The schematic diagram showing the in vitro flow assay in Cas9 expressing HSPCs infected with Myc-GFP and sgRNA-mCherry. (B) T7 endonuclease I assay showed the sgNpm1/Phgdh efficiency in c-Kit+ cells derived from Cas9 expressing mouse. (C) Western blotting showing the protein levels of NPM1 and PHGDH in HSPCs edited with sgScr or sgNpm1/Phgdh. (D) The relative expressions of Phgdh along with the differentiation (top) and transforming (bottom) trajectories. (E) The bar graphs showed the MFI of stem cells (left) and differentiated cells (right) at 2 days, 4 days, and 6 days, in sgPhgdh samples, compared to scramble samples. *p < 0.05, **p < 0.01, ***p < 0.001, p values were calculated using an unpaired parametric t test. (F) The box plot showing the expression levels of GMP-like signature (left) and GMP signature (right) in sgPhgdh and sgScramble samples. *p < 0.05, **p < 0.05, p values were calculated by likelihood ratio test. (G) The scatter plot showing differentiation/stemness signature scores (left) and G2M/S scores (right) in sgPhgdh and sgScramble samples. Measured by RNA-seq. (H) The KEGG pathways enriched in sgScramble samples, compared to sgPhgdh cells (top). The KEGG pathways enriched in sgPhgdh samples, compared to sgScramble cells (bottom). (I) GSEA showing the negative enrichment of HALLMARK_MYC_TARGETS_V1 and HALLMARK_MRNA_SPLICING in sgPhgdh cells, comparing to sgSCr cells. (J) Kaplan–Meier curve showing the survival of AML patients in TCGA stratified by the expression of PHGDH. P value was calculated by log-rank test. (K) The box plots showing the expression levels of NPM1 in the normal GTEx samples and TCGA cohorts. *p.adj < 0.05, ****p.adj < 0.0001, n.s. not significant, p values were calculated by Wilcoxon test. (L) The box plots showing the expression levels of PHGDH in the normal GTEx samples and TCGA cohorts. *p.adj < 0.05, ****p.adj < 0.0001, n.s. not significant, p values were calculated by Wilcoxon test. The underlying data for S5E–S5L Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s005

(TIF)

S6 Fig. The molecular characteristics of intertumoral heterogeneity.

Related to Fig 5. (A, B) The force-directed layout maps showing leukemogenesis trajectories, colored by pseudo time (A) and the expressions of human AML signature genes (B). (C) The bar graph showing the proportion of each cell type during leukemogenesis. (D) The heatmap showing the dynamics molecular features during leukemia progression. The key genes and top enriched gene sets were labeled. (E) The box plot showing the stemness scores and differentiation scores in PLCs. ****p.adj < 0.0001, p values were calculated by Wilcoxon test. (F) The heatmap showing the expression levels of normal cell subtypes’ signatures in PLCs. (G) Heatmap summarized the expression of 5 malignant cell subtypes from T3 in TCGA-LAML database annotated by clinical information including mutation counts, sex, blast percentages, WBC, diagnosis ages, FAB, and subgroups of patients. (H) Hazard ratio calculated using the expression of 5 malignant cell subtypes from T3 in TCGA-LAML database. Hazard ratios >1 indicate an increased risk of dying. While hazard ratios <1 indicate a beneficial prognosis for the patient; p values of each individual factor based on the multivariate analysis is depicted on the right of the figure with the values: *p < 0.05, **p < 0.01, ***p < 0.001. (I) Heatmap summarized the expression of 5 malignant cell subtypes in TARGET-AML database annotated by clinical information including mutation counts, sex, blast percentages, WBC, diagnosis ages, FAB, and subgroups of patients. (J) Hazard ratio calculated using the expression of 5 malignant cell subtypes from T3 in TARGET-AML database. Hazard ratios >1 indicate an increased risk of dying. While hazard ratios <1 indicate a beneficial prognosis for the patient; p values of each individual factor based on the multivariate analysis is depicted on the right of the figure with the values: *p < 0.05, **p < 0.01, ***p < 0.001. (K) The box plots showing the expression levels of HSP90AA1 in the normal GTEx samples and TCGA cohorts. *p.adj < 0.05, ****p.adj < 0.0001, n.s. not significant, p values were calculated by Wilcoxon test. (L) The scatter plots showing the correlation between the expression levels of Npm1 or Phgdh and the Hsp90aa1hi signatures. The underlying data for S6E and S6K Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s006

(TIF)

S7 Fig. The similarity of the intra-heterogeneity and its associated molecular features between human and murine AML.

Related to Fig 5. (A) The heatmaps showed the results of inferCNV with hierarchical clustering in human AML. (B) The bar graph showing the proportion of 5 AML patients in each leukemia subtypes. (C) Top enriched KEGG pathways for differentially expressed genes in malignant cells. Node size represents gene ratio; node color represents p.adjust. p.adjust, adjusted p value.

https://doi.org/10.1371/journal.pbio.3002088.s007

(TIF)

S8 Fig. Alternative splicing of TMEM134 is common in murine and human malignant cells.

Related to Fig 6. (A) The t-SNE plot of T3 Myc-GFP-positive cells analyzed by Smart-seq2, colored by cell types. (B) The regulatory networks of signature genes in Hsp90aa1hi cells. (C) The bar graph showing the numbers of differentially alternative splicing events between Hsp90aa1lo and Hsp90aa1hi cells. (D) The bar graph showing the skipping ratio of Tmem134 exon6 in HSPC and AML. (E) The PSI values of exon6 skipping of TMEM134 in TCGA pan-cancer samples. The underlying data for S8C–S8E Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s008

(TIF)

S9 Fig. Alternative splicing of TMEM134 promoted cell cycle.

Related to Fig 6. (A) Semiquantitative PCR showed the relative mRNA expression of TMEM134α or TMEM134β in c-Kit+ HSPCs harboring vector, TMEM134α and TMEM134β. (B) Percentage of GFP-positive cells in peripheral blood collected from recipient mice with vector, TMEM134α and TMEM134β after 2- or 4-week transplantation. (C) T7 endonuclease I assay showed sgTmem134 efficiency in Myc-induced murine AML cells. (D) Sequence alignment result showing 2 bp insertion in one of the Tmem134-/- AML clones. (E) Western blotting showing of TMEM134 protein levels in Tmem134 wild-type (Myc AML) or knockout leukemic cells (M8#). (F) Immunofluorescence staining of TMEM134 protein in Tmem134 wild-type (Myc AML) leukemic cells or M8#. (G) Semiquantitative PCR showed the relative mRNA expression of TMEM134α or TMEM134β in endogenous Tmem134 KO mouse AML cells M8# harboring vector, TMEM134α and TMEM134β. (H) Semiquantitative PCR showed the relative mRNA expression of Tmem134α or Tmem134β in Myc-induced leukemic cells harboring sgScramble, sgExon6, or sgIntron5+6. (I) T7 endonuclease I assay showed sgTMEM134 efficiency in THP-1 cell line. (J) Sequence alignment results showing 2 bp insertion and 21 bp deletion in one of the Tmem134-/- THP-1 clones and 71 bp deletion in one of the Tmem134-/- HL-60 clones. (K) The sequence of synonymous mutated (PAM site) TMEM134α/β cDNA rescued in Tmem134 KO cell lines. (L) The growth curves of endogenous TMEM134 KO HL-60 cells with overexpressing vector, TMEM134α or TMEM134β at the indicated time points. **p < 0.01, ***p < 0.001. P values were calculated by unpaired parametric t test. (M) The Kaplan–Meier survival curves of recipient mice transplanted with leukemic cells overexpressing vector, TMEM134α or TMEM134β. **p < 0.01(log-rank test). (N) The bar graph showing the proliferative ability of Myc-GFP AML cells transduced with sgScramble, sgExon6, or sgIntron5+6, with EdU incorporation assay. ***p < 0.001, p values were calculated using an unpaired parametric t test. (O) The bar graphs showed G1 score between Hsp90aa1lo and Hsp90aa1hi cells in murine AML (left), human AML (middle), and TCGA-LAML (right). *p < 0.05, p values were calculated by Wilcoxon signed-rank test. (P) The heatmap showing GSVA enrichments of cell cycle related pathways (rows) in vector, TMEM134α or TMEM134β AML cells (columns). (Q) The heatmap showing the expression levels of genes in TAKEDA_TARGETS_OF_NUP98_HOXA9_FUSION_8D_UP pathway. The underlying data for S9B and S9L–S9O Fig can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.3002088.s009

(TIF)

Source link