Single-cell analysis reveals altered tumor microenvironments of relapse- and remission-associated pediatric acute myeloid leukemia

Abstract

Acute myeloid leukemia (AML) microenvironment exhibits cellular and molecular differences among various subtypes. Here, we utilize single-cell RNA sequencing (scRNA-seq) to analyze pediatric AML bone marrow (BM) samples from diagnosis (Dx), end of induction (EOI), and relapse timepoints. Analysis of Dx, EOI scRNA-seq, and TARGET AML RNA-seq datasets reveals an AML blasts-associated 7-gene signature (CLEC11A, PRAME, AZU1, NREP, ARMH1, C1QBP, TRH), which we validate on independent datasets. The analysis reveals distinct clusters of Dx relapse- and continuous complete remission (CCR)-associated AML-blasts with differential expression of genes associated with survival. At Dx, relapse-associated samples have more exhausted T cells while CCR-associated samples have more inflammatory M1 macrophages. Post-therapy EOI residual blasts overexpress fatty acid oxidation, tumor growth, and stemness genes. Also, a post-therapy T-cell cluster associated with relapse samples exhibits downregulation of MHC Class I and T-cell regulatory genes. Altogether, this study deeply characterizes pediatric AML relapse- and CCR-associated samples to provide insights into the BM microenvironment landscape.

Fig. 1. Single-cell transcriptional profiling identifies heterogenous putative AML blast cell clusters.

a UMAP embedding of the four paired Dx and EOI biologically independent samples (three with relapse; patients 3, 5, 6, and one with CCR; patient 14) consisting of 19,350 high quality single cells partitioned into fourteen clusters. These clusters are colored based on expression of canonical cell type(s) gene markers (Fig. 1b). b Dot plot showing expression of different cell type specific marker genes that were used to annotate the cell types. The color intensity represents gene expression with red: high, yellow: medium, cyan: low, and size of dot represents percentage of cells expressing each gene in individual cell types. c Relative proportion of different patient samples in each cell type cluster. Most of the blast cell clusters are made up of cells from one or two patient samples indicating AML blast heterogeneity. On the other hand, most of the immune cell clusters contain cells from multiple patient samples. D denotes Dx samples, E denotes EOI samples. d Split UMAP showing the putative blast cells that are significantly over-represented in Dx samples (Dx enriched: ~95% of total blast cell clusters are made up of Dx samples with remaining ~5% made up of EOI samples; EOI enriched: non-blast cells making ~95% of EOI samples). e Heatmap of the selected set of 44 genes showing significant overexpression (Wilcoxon rank sum test, two-sided; Fold Change >1.2, P < 0.01) in Dx enriched blast cells in comparison to EOI enriched differentiated cells. Relative gene expression is shown with blue and red colors representing low and high expressing genes, respectively. Columns and rows represent the cells and genes, respectively. Patient samples used are shown with colored human icons at the bottom of the figure.

Fig. 2. Development of blast progenitor signature from heterogeneous AML-blasts.

a Heatmap of 20 genes that show significant overexpression in the diagnosis (Dx) TARGET AML-1 samples with high blast enrichment (> 60% blasts and 30–60% blasts) as compared to Dx low (<30%) blasts and end of induction (EOI) samples. The Y-axis represents gene names and X-axis represents the patient’s ID. Color scale: blue and red colors represent low and high expression of genes, respectively. b Feature plots of select genes (CLEC11A, NREP) from 20-genes set showing uniform overexpression in the lassoed blast clusters (UMAPs of Dx, EOI samples from patients 3, 5, 6, 14). Color scale represents high (red), medium (yellow), and low (gray) expression of the select genes. Violin plots of select genes (CLEC11A, NREP) expression in Dx (n = 4) and EOI (n = 4) samples, showing overexpression in the former. c Representative expression profiles of select genes (CLEC11A, NREP) show a progressive downregulation pattern from high blast % to low blast % to EOI samples in TARGET AML-1 dataset (biologically independent samples in >60% blasts group, n = 203; 30–60% blasts group, n = 60; <30% blasts group, n = 14; EOI group, n = 24, two-sided Wilcoxon rank sum test, **P < 0.01, ***P < 0.001). Comparisons of normalized CLEC11A expression in >60% blasts vs. <30% blasts (P = 9.6e−05) and >60% blasts vs. EOI (P = 1.1e−08) groups and NREP expression in >60% blasts vs. <30% blasts (P = 0.0052) and >60% blasts vs EOI (P = 0.00018) groups had significant results. Boxplots show the distribution of expression with the center of the box representing the median, upper and lower bounds representing 75% and 25% percentiles, and upper and lower whiskers extending to the largest value no further than 1.5 times interquartile range from bounds of box. Source data are provided as a Source data file. d Survival correlations showed the higher expression of CLEC11A (biologically independent samples: High: n = 23, Low: n = 122; HR = 1.96, 95% Confidence Interval (CI): 1.10–3.48, Cox HR P = 0.022, one-sided log-rank P = 0.0197) and NREP (biologically independent samples: High: n = 60, Low: n = 85; HR = 1.95, 95% CI: 1.21–3.14, Cox HR P = 0.0064, one-sided log-rank P = 0.0055) is associated with poor survival. One-sided log-rank tests were used to compare survival curves of high and low expression groups. Cox proportional hazards models were used to calculate hazard ratios and Wald tests were used to determine significance. e Receiver operative Curve (ROC) depicting performance of Support Vector Machine (SVM) approach based 7-gene signature (CLEC11A, PRAME, AZU1, NREP, ARMH1, C1QBP, TRH) discriminated AML-blasts from other cells. AUC (Area under curve) is 0.968, and gray shaded area around curve represents 95% CI (0.941–0.989). f UMAP plot depicting the blast cells (lassoed) of the longitudinal samples (Dx, EOI, Rel) from two patients (patient 5 and 6). g Bar plot showing fraction of Dx, EOI, and Rel samples of patients 5 and 6 in each cell type. D denotes Dx samples, E denotes EOI samples, R denotes relapse samples. h UMAP of integrated Dx, EOI AML samples (patients 3, 5, 6, 14), and healthy BM and normal HSCs from external datasets. i. Feature plot of the 7-gene signature module score in UMAPs split by dataset (AML, healthy BMs, normal HSCs). Color scale: purple represents a high module score, i.e., the 7-gene signature is highly expressed by the cells, and yellow represents a low module score. j Performance of SVM classifier as AUC using ROC analysis on independent validation set of Dx, EOI samples from five patients (patients: 16–20). AUC is 0.784, and gray shaded area around curve represents 95% CI (0.739–0.831). k Performance of SVM classifier as AUC on prospective publicly available validation dataset of pediatric AML, healthy BM, and normal HSC samples. AUC is 0.809 (95% CI: 0.734–0.866).

Fig. 3. AML-blast cells from samples with relapse depicted differences in transcriptome profile in comparison to samples with CCR.

a Split UMAP plot of Dx AML blast cells of relapse- and CCR-associated biologically independent samples from patients’ 1–14 (shown with colored human icons at bottom of the figure). Blast cells formed fifteen clusters based on transcriptome profiles and were labeled manually based on top genes. b Stacked bar plots showing the proportion of cells from each patient in different clusters (Red: Relapse-associated, Blue: CCR-associated). D denotes Dx samples. c Heatmap of top differentially expressed genes between select relapse-associated (0, 2, 7, 13) and CCR-associated (1, 3, 4) samples clusters. Color scale: blue and red colors representing low and high relative expression of genes, respectively. The highlighted genes have been discussed in the main text for association with leukemia development and progression by regulating key pathways. d Kaplan–Meier plots show high expression of DEGs such as FLNA and RFLNB/FAM101B in AML-blasts, having more  expression in relapse-associated samples, were associated with poorer OS in TARGET AML dataset. Similar analysis on genes (MPO, TRH) with greater expression in the CCR-associated clusters depicted significant association of high expression with better OS. Number of samples in high and low expression groups are shown on plots. One-sided log-rank tests were used to compare survival curves of high and low expression groups. Cox proportional hazards models were used to calculate hazard ratios and Wald tests were used to determine significance. e Pathways that were significantly (P < 0.01) activated (Z-score >1.5)/inhibited (Z score < −1.5) in samples with relapse vs. samples with CCR. Activation and inhibition of pathways was determined using one-tailed Fisher’s Exact tests. f Upstream regulatory molecules significantly activated (orange) in the blasts enriched with relapse-associated samples.

Fig. 4. Immune microenvironment analysis of Dx samples.

Non-blast cells from biologically independent samples of patients 1–14 (shown with colored human icons at bottom of the figure) were selected to perform a focused analysis. a Split UMAP plot of BME cells based on relapse- and CCR-associated status of samples. b Expression-based dot plot of canonical cell-specific markers used for the annotation of lymphoid, myeloid, and erythroid lineages. The color intensity represents gene expression with red: high, yellow: medium, cyan: low, and size of dot represents percentage of cells expressing each gene in individual cell types. c t-distributed Stochastic Neighbor Embedding (t-SNE) plot depicting fifteen subclusters of T-lymphocytes. The subclusters were annotated based on the canonical cell type markers: Naive T cells (CCR7⁺, LEF1⁺, TCF7⁺), CD4⁺ effector T cells (CD4⁺, CCR6⁺, CXCR6⁺, CCL5⁺), CD8⁺ cytotoxic T cells (CD8A⁺, CD8B⁺, GZMB⁺, GNLY⁺, PRF1⁺) and exhausted T cells (HAVCR2⁺, LAG3⁺, PDCD1⁺, NFATC1⁺, TIGIT⁺, TOX⁺). Clusters enriched with cells from CCR, and cells from relapse samples are lassoed. d Stacked bar plots showing the cluster-wise proportion of cells from each patient (Red: relapse, Blue: CCR). e Proportion of T cells in relapse- and CCR-associated patients (top panel; CCR: n = 8 biologically independent samples, relapse: n = 6 biologically independent samples, two-tailed unpaired t-test, *P < 0.05. % T cells: P = 0.0452, % T cells of CD45+ : P = 0.0470). Overall enrichment of naive and exhausted T cell signatures across relapse- and CCR-associated samples based on ssGSEA score (bottom panel; n = 599 cells from n = 8 biologically independent CCR samples; n = 1773 cells from n = 6 biologically independent relapse samples; two-tailed unpaired t-test, ****P < 0.0001). Bar plots show the mean value, with error bars representing mean + standard error of mean (SEM). Source data are provided as a Source data file. f Heatmap of top DEGs in clusters 3, 5, 6 having more relapse-associated samples and clusters 7, 8, 10 having more CCR-associated samples (the T cell subclusters that were made up of multiple samples were selected for this analysis). Relative gene expression is shown with blue and red colors representing low and high expressing genes. g Pathways that are significantly activated (Z-score >1.5) and inhibited (Z-score < −1.5) with genes differentially expressed in the relapse-associated T cells. Pathways achieved P < 0.01 based on one-tailed Fisher’s Exact test.

Fig. 5. CCR-associated samples at diagnosis depicted enrichment of inflammatory monocytes/macrophages.

We conducted the focused analysis of monocyte/macrophages clusters from biologically independent Dx samples (patients 1–14, shown with colored human icons at bottom of the figure). a UMAP plot depicting monocytes/macrophages clusters (lassoed) in Dx samples. Cluster 1 (Mono/mac 1) and 2 (Mono/mac 2) are enriched with cells from CCR-associated patient samples whereas cluster 6 (mono/mac 6) is enriched with cells from relapse-associated samples. b SsGSEA analysis shows significantly increased expression of M1 macrophages genes (S100A8, S100A9, S100A12, TYROBP, VCAN, CD68, MNDA, CYBB, STAT1) in CCR-associated samples compared to relapse-associated samples (clusters 1, 2, and 6). n = 5810 cells from n = 8 biologically independent CCR samples, and n = 1223 cells from n = 6 biologically independent relapse samples, two-tailed Wilcoxon rank sum test (P = 1.80e−177), ****P < 0.0001. Boxplots show the distribution of enrichment scores with the center of the box representing the median, upper and lower bounds representing 75% and 25% percentiles, and upper and lower whiskers extending to the largest value no further than 1.5 times interquartile range from bounds of box. Source data are provided as a Source data file. c Heatmap of top DEGs from the comparison of relapse- (cluster 6) and CCR-associated dominant (clusters 1, 2) clusters. Relative gene expression is shown with blue and red colors representing low and high expressing genes. The highlighted genes are the ones mentioned in the main text. d Upstream regulatory molecules significantly inhibited (blue) and activated (orange) in the CCR-associated samples enriched monocytes/macrophages clusters in comparison to the relapse-associated cluster. e Split violin plots showing expression levels of genes encoding specific regulators in CCR- (clusters 1, 2) and relapse-associated (cluster 6) samples.

Fig. 6. ScRNA-seq analysis of EOI samples identified post-therapy residual blast cells with distinct transcriptome landscapes.

a Split UMAP plot shows the putative blast cells (colored in cyan) that are over-represented in Dx samples (patients 3, 5, 6, 14) and reduced in the EOI samples (patients 3, 5, 6, 14, 15). b Heatmap of top markers that are differentially expressed between EOI post-therapy residual and Dx blast cells. Relative gene expression is shown with blue and red colors representing low and high expressing genes. The highlighted genes are the ones mentioned in the main text. c Split violin plots of pathways significantly enriched in the EOI residual blast cells compared to Dx blast cells (Dx: n = 4173 cells from n = 4 biologically independent samples, EOI: n = 264 cells from n = 5 biologically independent samples, linear model comparison of means, ***P < 0.001, exact P-values are located in Supplementary Table 11). d Regulators that are significantly activated (orange color) and inhibited (blue color) in EOI residual blasts (n = 264 cells) compared to treatment responsive Dx blasts (n = 4173 cells). The significance of transcriptional regulators was determined using one-tailed Fisher’s exact test. The regulators with P < 0.01 and absolute Z-score of ±1.5 were considered statistically significant.

Fig. 7. EOI non-blast cells analysis exhibits different patterns in samples from those with relapse and those with CCR.

a UMAP plot of annotated non-blast cells clusters split based on relapse- and CCR-association status. b Bar plot showing the proportion of cell types in relapse- and CCR-associated samples at EOI. The Y-axis represents the fraction of each cell type. c Bar plot showing sample/patient fraction (Y-axis) in major cell types (X-axis). d Bar plot of mean scaled ssGSEA scores of mono/mac cells from those with relapse and those with CCR (CCR: n = 2203 cells from n = 2 biologically independent samples, relapse: n = 2416 cells from n = 3 biologically independent samples) show that M1 markers are significantly more enriched in patients with CCR (two-sided t-test, P = 3.8e−05, ***P < 0.001) and M2 markers are more enriched in those with relapse (two-sided t-test, P = 0.046, *P < 0.05). Source data are provided as a Source data file. e The bar plot shows that the ratio of the naive T cell-2 cluster to total T-cells (i.e., naive T cell-2 and naive T cell-1) is significantly higher (P = 0.0232) in samples with relapse (n = 3) compared to samples with CCR (n = 2) (two-sided t-test, *P < 0.05). Bars represent the mean values, with error bars representing the mean + SEM values. Source data are provided as a Source data file. f Bar plot of significantly activated (Z-score >1.5) or inhibited (Z-score < −1.5) pathways in T cell-2 cluster (relapsed-enriched, n = 1538 cells) compared to T cell-1 cluster (n = 2272 cells). The significantly (One-tailed Fisher’s exact test, P < 0.05) activated and inhibited pathways are shown with red and blue colors, respectively.