A Single-Cell Taxonomy Predicts Inflammatory Niche Remodeling to Drive Tissue Failure and Outcome in Human AML

Cancer initiation and drug resistance are orchestrated by an interplay between tumor-initiating cells, residual normal tissue stem and progenitor cells, and ancillary tissue-resident cells, comprising the stromal and immune environment (1). Ultimately, insight into the signaling cascades operative between all these cells in a neoplastic state is required to obtain a comprehensive appreciation of the complexity of tumorigenesis and tumor survival in the context of therapy.

One important component of all human tissues are stromal cells. Bone marrow stromal cells (BMSC) pervade the tissue in extensive networks that are estimated to comprise up to 20% of the marrow's cellular volume and associate with the vast majority of hematopoietic cells (2). They represent a likely heterogeneous population of cells, with a subset considered critical for the maintenance and homeostatic regulation of hematopoietic stem/progenitor cells (HSPC; ref. 3). A subset of stromal cells comprising the HSPC niche has been characterized in mice by expression of the leptin receptor (LEPR; refs. 4, 5) and high expression of CXCL12 [hence dubbed “CXCL12-abundant reticular (CAR) cells”; refs. 6, 7] located in the vicinity of sinusoid vessels (2, 5, 8). These stromal cells are major sources of stem cell factor (SCF) and interleukin-7 and are considered critical regulators of hematopoietic stem cells (HSC), multipotent progenitors, lymphoid progenitors as well as natural killer (NK), B and plasmacytoid dendritic cell (DC) development (4, 9–12).

It is, however, important to emphasize that these insights have been derived mostly from genetic studies in nonhuman species, particularly mice. The biology, heterogeneity, and interactions of LEPR⁺ stromal cells in the human bone marrow (BM) have remained largely elusive, in part because it is difficult to retrieve these cells in sufficient numbers, likely because they reside along fibers of the extracellular matrix, limiting their capture by aspiration.

Insight into the biology of human LEPR⁺ stromal cells, and in particular their interaction with HSPCs, is of likely further relevance to the biology of myeloid neoplasms, including acute myeloid leukemia (AML). AML is caused by genetic events occurring in HSPCs, but mouse modeling and humanized ex vivo modeling have implicated stromal niche alterations in the initiation, maintenance, drug resistance, and progression of AML (13–23). Transcriptional alterations in Lepr⁺ stromal cells, including the downregulation of hematopoietic regulatory factors, have further been documented in mice transplanted with MLL-AF9 leukemic cells, suggesting that a leukemic state may attenuate the supportive function of Lepr⁺ stromal cells for normal hematopoiesis (24).

The relevance of these proposed concepts and mechanisms for human disease have, however, remained uncertain, primarily because we lack comprehensive insights on LEPR⁺ BMSCs in the human leukemic marrow and their potential pathologic interaction with hematopoietic elements.

Here, we performed comprehensive, paired, single-cell transcriptional sequencing and networking analyses of all cells in human normal bone marrow (NBM) and AML aspirates.

The data establish a comprehensive taxonomy of the predicted cellular interactions between LEPR⁺ stromal niches, HSPCs, and adaptive and innate immune cells in the human NBM and the disruption of these signaling pathways defining AML. Induced deterioration of stromal niche support in AML through inflammatory activation of LEPR⁺ cells is predicted to broadly affect tissue signaling and induce repression of normal hematopoiesis via disruption of signaling toward high-output HSCs and committed progenitors but is ultimately associated with a favorable prognosis. The data provide human disease relevance to previous findings in model systems and provide a resource of intercellular signaling in the human BM defining native and neoplastic hematopoiesis, with a particular emphasis on the communication between HSPCs and their stromal niches.

To generate a cellular taxonomy of the human NBM, representing both rare HSPCs and stromal niche populations, allowing assessment of their cellular diversity, and predicted intercellular signaling, we performed single-cell RNA sequencing (scRNA-seq) on viably frozen BM aspirates from four healthy donors for allogeneic transplantation (median age 51; range, 47–53 years; Supplementary Table S1). To ensure a robust representation of all BM cell types in the dataset (enabling accurate assessment of population heterogeneity at high resolution), we used flow cytometry to sort and purify all cells in the aspirates into five fractions: the nonhematopoietic stromal (CD45⁻CD235a⁻CD71⁻CD31⁻) and nonstromal/endothelial (CD45⁻CD235a⁻CD71⁻CD31⁺) fraction, the HSPC fraction (CD45⁺CD34⁺), the myeloid fraction (CD45⁺CD34⁻CD117⁺/CD33⁺), and the nonmyeloid (lymphoid) fraction (CD45⁺CD34⁻CD117⁻CD33⁻; Supplementary Fig. S1A). Fractions were combined into two pools (HSPC with myeloid and nonhematopoietic with nonmyeloid) and subjected to scRNA-seq in separate runs, and the data of the two separate runs were subsequently merged into a single scRNA-seq dataset (Supplementary Fig. S1A).

We acquired high-quality data from 46,740 NBM cells. The purification strategy resulted in a robust representation of low-frequency cell populations such as nonhematopoietic stromal cells and HSCs in the dataset (Fig. 1A).

Using Clustifyr, a package to classify cells from scRNA-seq data using external references (25), the sequenced cells were classified into BMSCs, endothelial cells (EC), CD34⁺ HSPCs, erythroid progenitors, megakaryocytes (MK), DCs, CD14⁺ monocytes (CD14⁺ Mono), CD16⁺ monocytes (CD16⁺ Mono), NK cells, CD4⁺ T cells (CD4 T), CD8⁺ T cells (CD8 T), B cells (B), and plasma cells (Supplementary Fig. S1B). The accuracy of the unsupervised cell annotation was validated by checking the expression of canonical markers for each cell type (Fig. 1B; Supplementary Fig. S1C and S1D). The BMSC cluster, for instance, specifically expressed BMSC markers NGFR (CD271), PRRX1, the perivascular MSC marker LEPR and the HSPC niche factors CXCL12 and KITLG (Fig. 1B; Supplementary Fig. S1C and S1D). The HSPC cluster expressed CD34, AVP (26), and KIT (Fig. 1B; Supplementary Fig. S1C and S1D) and the CD4⁺ T, CD8⁺ T, and NK clusters had enhanced expression of CD3D, CD8B, and NKG7, respectively (Fig. 1B; Supplementary Fig. S1C and S1D).

To dissect the cellular heterogeneity within the HSPC fraction, we further classified this subset using the K nearest-neighbor classification method, resulting in distinct subtypes that were subsequently clustered into HSC/multipotent progenitors (MPP), lymphomyeloid primed progenitors (LMPP), MK and erythroid progenitors (MEP), erythroid progenitors, common lymphoid progenitors (CLP), and granulocyte–monocyte progenitors (GMP)/monoblasts (enriched) subpopulations (Supplementary Fig. S2A). Annotation of each HSPC subpopulation was based on comprehensive transcriptional analysis, including CD34 expression (highly expressed in HSC/MPPs, LMPPs, CLPs but gradually deceased in GMP/monoblasts and erythroid progenitors; Supplementary Fig. S2B), genes indicative of cell-cycle status (cells in the HSC cluster were predominantly retained in a noncycling (G₁) phase; Supplementary Fig. S2B), gene enrichment scores to identify HSPC subsets (ref. 27; Supplementary Fig. S2C) and expression of genes indicative of HSPC type (HLF, AVP, and HES1 for HSCs, myeloid markers LYZ, MPO, and S100A8/A9 for GMPs, B cell lineage marker CD79A for CLPs and erythroid markers CA1, HBB, and HB for MEPs; Supplementary Fig. S2D and S2E; ref. 26). Trajectory analysis was congruent with the view that HSCs are at the basis of a cellular hierarchy giving rise to erythroid, myeloid, and lymphoid progeny (Supplementary Fig. S2F).

scRNA-seq studies in mice have revealed transcriptional differences among HSCs, providing an explanation for their functional heterogeneity. Distinct transcriptomes discern long-term (LT), low-output HSCs that self-renew rather than proliferate or differentiate from short-term (ST), high-output HSCs that provide lineage output toward mature progeny during normal steady-state hematopoiesis (28–30). Whether a similar hierarchy of HSCs, reflected in their transcriptome, exists within the HSC pool in human native hematopoiesis has remained largely unknown.

Subclustering of the human HSC/MPP population distinguished 5 subsets (HSC 0–4; Fig. 1C). Cluster 0 constitutes 23.85% ± 1.18% of the HSC pool (Fig. 1C). The transcriptional signature of this cluster 0 had a strong relationship with the transcriptional signature previously reported for low-output LT-HSCs in mice (28, 31), whereas clusters 1–4 displayed strong transcriptional congruence with signatures of high-output ST-HSCs (Fig. 1D).

The “low-output” HSC cluster 0 was also enriched for HSC signatures previously associated with platelet/MK-biased HSCs (28, 31), previously demonstrated to reside at the apex of the hematopoietic stem cell hierarchy (ref. 32; Fig. 1D).

Transcriptional similarity of cells in cluster 0 with murine LT-HSCs included expression of genes associated with self-renewal and quiescence such as Mllt3, Socs2, Txnip, and Ndn, MHC class II components [Cd74, H2-Eb1(HLA-DRB1)], and transcription regulators (Hlf, Hes1; Fig. 1E; Supplementary Table S2). Significantly differentially expressed genes in clusters 1–4 included CDK6, previously reported to be a critical regulator for the transition of LT-HSCs to ST-HSCs and MPPs (33, 34) and FLT3 (Fig. 1E).

Gene set enrichment analysis (GSEA) demonstrated transcriptional activation of interferon and STAT activation–related inflammatory signaling (TNF signaling via NFκB; interferon gamma response; IL2–STAT5 signaling) in the HSC-0 cluster (Fig. 1F), reminiscent of an inflammatory transcriptional signature defining a subset of CD74^high MHCII^high HSCs resistant to myeloablation in mice (35). Interestingly, human cluster-0 HSCs display concomitant overexpression of genes encoding proteins inhibiting inflammation, including TNFAIP3 (A20; ref. 36) and all members of the NR4A subfamily of nuclear receptors (NR4A1, NR4A2, and NR4A3; Fig. 1E), previously shown to protect LT-HSCs against DNA damage and repress the proliferative inflammatory response of HSCs (37, 38), whereas cell-cycle-related Hallmark gene signatures (E2F targets and G₂–M checkpoint) were significantly enriched in clusters 1–4 (Fig. 1F). Cell trajectory analysis was congruent with the notion that the HSC-0 subset may be at the base of a cellular hierarchy with differentiation trajectories toward HSC clusters 1–4 (Fig. 1G).

Collectively, the data reveal the existence of distinct transcriptional subsets within the HSC/MPP population in human native hematopoiesis, suggesting a linear hierarchy from low-output LT-HSCs to high-output ST-HSC and MPPs, that has previously been demonstrated in mice.

The generation of a comprehensive cellular taxonomy allows the assessment of predicted cellular interactions in the human BM. We analyzed the predicted intercellular communication between all retrieved cell types based on possible ligand–receptor cross-talk with CellChat, which predicts the major signaling routes and how they integrate into cellular function, using network analysis and pattern recognition approaches (39), revealing a complicated cell interaction network in which all cell types were involved (Fig. 2A).

BMSCs were identified as a key “communication hub” in these analyses, with predicted signaling routes to almost all other cell types in the human marrow, with a predicted magnitude of signaling interaction exceeding that of any other cell population in the marrow (Fig. 2A). This predicted signaling included well-established interactions of BMSCs with HSPCs via KITLG (SCF; to HSC/MPPs, LMPPs, MEPs, and erythroid progenitors; ref. 9) and all immune subsets via CXCL12 (10, 11, 40), lymphoid cells via IL7 (to CLPs, CD4⁺ T, CD8⁺ T, and NK; refs. 4, 41) myeloid cells via CSF1 and IL34 (to GMP/monoblasts, DCs, CD14⁺ monocytes, and CD16⁺ monocytes; ref. 42), ECs via VEGF and ANGPTL (43), and predicted autocrine signaling via IGF, BMP, TGFβ, and adiponectin, previously implicated in fate-decision signaling (Fig. 2B).

Collectively, the data describe a cellular taxonomy of the stromal and immune environment of the human NBM and their predicted intercellular signaling, unveiling a predicted central role of BMSCs, not only in the homeostatic regulation of HSCs and HSPC subsets but also all innate and adaptive immune cells, in line with experimental data from mouse studies (9, 12).

BMSCs, as shown in mice, represent a heterogeneous population with pleiotropic function in tissue homeostasis and support of distinct hematopoietic cell types (44, 45). Heterogeneity is, in part, related to anatomic localization (endosteal vs. medullary localized BMSCs), but also within the medullary fraction of BMSCs (likely best represented in BM aspirates), heterogeneity exists in mice (2). Whether such heterogeneity exists in humans has remained largely unknown. The robust representation of stromal cells in our datasets allowed for addressing this question.

Subclustering of the BMSC population distinguished 4 subsets (BMSC-0-3; Fig. 2C). All BMSC clusters expressed the sinusoidal stromal HSPC niche marker LEPR (9), albeit at different levels, as well as preadipocytic markers (LPL, ADIPOQ, and CD36) and HSC niche factors (CXCL12, KITLG, and ANGPT1; Fig. 2D; Supplementary Fig. S3A), whereas expression of osteolineage differentiation markers (BGLAP, RUNX2, and SPP1), chondrocyte differentiation markers (SOX9, ACAN, and COL2A1), fibroblast markers (S100A4 and SEMA3C), and pericyte markers (NES, NG2, and ACTA2) could not be detected or were lowly expressed (Supplementary Fig. S3A).

These findings seem consistent with the notion that the BMSCs captured in our analyses represent the human equivalent of the sinusoidal Lepr-expressing stromal niche cells in mice that display adipogenic differentiation capacity (24, 46) and are essential for HSPC maintenance (4, 9). Other microenvironmental cells, such as ECs, pericytes, and osteoblasts, could not be retrieved from human BM aspirates or only in numbers too low to perform robust analyses (in the case of ECs).

The BMSC-0 subset, comprising 10.5% ± 5.11% of all stromal cells, displayed the highest expression of LEPR and genes encoding critical HSPC regulatory factors (CXCL12, KITLG, ANGPT1, and IL7) in comparison with other clusters (Fig. 2D). Transcriptome-wide comparison of this subset with the other stromal subsets revealed differential expression of 443 genes (144 upregulated and 299 downregulated, P_adj < 0.05; Supplementary Fig. S3B; Supplementary Table S3).

Ligand–receptor analysis revealed that the BMSC-0 subset had the strongest predicted KITLG–KIT, CXCL12–CXCR4, and IL7–IL7 receptor interaction with HSPCs, whereas these interactions were less strong in other clusters (BMSC-1, -2, and -3; Fig. 2E). BMSC-0–derived KITLG was predicted to signal to all HSC subsets (HSC0-4), whereas other interactions, in particular signaling via the CD74 receptor, were predicted to be strongest to the LT-HSC subset of HSCs (Supplementary Fig. S3C).

Taken together, the data indicate that heterogeneity exists within the human BMSC population with a minor subset, characterized by the highest expression of genes encoding critical HSPC maintaining factors, predicted to have the strongest interactions across HSPC subsets, particularly LT-HSCs.

Earlier approaches to establish cellular hierarchies in human AML have focused on hematopoietic cells (47), precluding the assessment of interactions between rare cell populations such as HSPCs and their stromal niches and their potential involvement in AML pathogenesis. We, thus, generated scRNA-seq tissue maps to establish the cellular taxonomy of the BM in AML with mutations in the gene encoding nucleophosmin (NPM1), among the most frequently mutated genes in AML, representing approximately 30% of newly diagnosed patients with AML (48).

BM aspirates from six patients with NPM1-mutant (NPM1m) AML at diagnosis (median age 53.5; range, 39–59 years; Supplementary Table S1) were sorted and subjected to scRNA-seq, according to the strategy described for NBM, resulting in transcriptomes of 49,758 cells comprising the taxonomy of the BM in NPM1m AML (Fig. 3A). The cells representing the clusters annotated within the NBM could largely be retrieved from the AML BM, although disruption of architecture within defined populations was observed, as described below.

The predicted central role of LEPR⁺ BMSCs as intercellular signal traffickers in the human NBM raised the question of whether, and how, LEPR⁺ BMSC characteristics are disrupted in AML and if alterations could contribute to disease pathogenesis.

Direct comparison of the BMSC transcriptome at the population level (outlined in Fig. 3A) in AML to normal controls showed significant (P_adj < 0.05) differential expression of 973 genes (541 genes upregulated and 432 genes downregulated; Fig. 3B; Supplementary Table S4). Among the most upregulated genes in AML were genes associated with inflammatory activation, including NFKBIA, PTGS2, CD44, CXCL2, CXCL8, TGFβ1, and ANXA1, and genes associated with extracellular matrix (ECM) remodeling, including LOLX2, TGFBI, COL5A3, and COL6A3 (Fig. 3B; Supplementary Fig. S4A and S4B; Supplementary Table S4), with significant downregulation of genes encoding critical HSPC regulatory genes such as KITLG, CXCL12, and IL7 and the gene encoding the HSPC niche marker LEPR and preadipocyte marker LPL (Fig. 3B; Supplementary Fig. S4A and S4B). This was reflected in gene signatures (GSEA Hallmark) consistent with inflammatory signaling (TNFα signaling via NFκB; IL2–STAT signaling; inflammatory response) among the top upregulated transcriptional signatures and downregulation of signatures associated with adipogenesis in AML BMSCs (Fig. 3C; Supplementary Fig. S4C). Disruption of metabolic pathways was suggested by gene signatures indicative of glycolysis in AML BMSCs (Hallmark Hypoxia; Hallmark Glycolysis) versus oxidative phosphorylation and fatty acid metabolism in normal BMSCs (Fig. 3C; Supplementary Table S4). In addition, the GO terms connective tissue development and cartilage development were also significantly increased in the AML BMSCs (Fig. 3D; Supplementary Fig. S4C), indicative of remodeling of the extracellular matrix.

To assess how this disruption of transcriptional programs in the overall BMSC population relates to BMSC heterogeneity in AML, we performed subclustering of the BMSC population (Fig. 3E). This revealed near loss of the cluster predicted to have the strongest interaction with HSPCs (BMSC-0; 3.53% ± 4.84% vs. 11.02% ± 3.39% in AML and NBM, respectively; P = 0.044), with a concomitant relative increase of BMSC-2 (46.4% ± 6.81% vs. 7.85% ± 1.24%; P = 0.009; Fig. 3E and F; Supplementary Fig. S4D). Levels of genes encoding critical HSPC maintaining factors were significantly reduced in all BMSC subsets in AML, which was most pronounced for KITLG and IL7, with dramatic impaired expression in BMSC-0 (Fig. 3G).

The BMSC-2 cluster in AML was transcriptionally characterized by the upregulation of genes and transcriptional programs associated with inflammation (Hallmark: Inflammatory Response) and genes and signatures related to ECM remodeling (Fig. 3H and I). This pattern of inflammatory disruption of BMSC architecture and reduced expression of HSPC factors was consistent among all six NPM1m AML samples examined (Supplementary Fig. S4A and S4D).

Stromal inflammation was confirmed in situ by demonstra­ting increased CD44 protein expression, which is a marker of stromal activation and inflammation (49), in CD271⁺(CXCL12⁺) BMSCs in AML by immunofluorescence on biopsies and flow­cytometric assessment of aspirates (Fig. 3J and K).

Trajectory analysis suggested that the relative reduction of BMSC-0 and increase in the inflammatory BMSC-2 subset may reflect a linear hierarchy characterized by a gradual increase in inflammatory activation accompanied by a loss of expression of HSPC regulatory factors (Supplementary Fig. S4E).

Taken together, the data are congruent with the view that in NPM1m AML, BMSCs are remodeled by inflammatory activation, resulting in a dramatic expansion of an inflammatory subset with a concomitant loss of the BMSC subset predicted to support the maintenance of normal HSPCs.

Next, we sought to explore the consequences of BMSC remodeling in AML for the residual normal and leukemic hematopoiesis. AML is characterized by suppression of normal hematopoiesis and expansion of clonal cells, but the cellular and molecular mechanisms promoting these processes have remained incompletely understood. Mouse modeling has suggested important contributions of stromal HSPC niches, but the relevance for human AML has remained uncertain.

The localization of NPM1 mutations at the 3′ terminal coding region of the gene, captured by polyA-RNA-seq, allowed us to distinguish NPM1-mutant AML cells from their nonleukemic hematopoietic counterparts, by assessing the mutational state of NPM1 in all cells in the dataset (Fig. 4A and Methods). Analysis of mutation status confirmed that adaptive immune subsets (CD4⁺/CD8⁺ T cells, plasma cells, and B cells), NK cells, and BMSCs are of nonleukemic origin in AML. The very small fraction of NPM1m cells within these fractions differentially and highly expressed myeloid genes (AZU1, LYZ, MPO, and ELANE; Supplementary Fig. S4F), indicative of blast contamination. NPM1m leukemic cells reside predominantly within the LMPP and GMP subsets of the HSPC fraction (Fig. 4A and B).

The remodeling of stromal niches in human AML, in particular the loss of stromal subsets predicted to have the strongest supportive interaction with HSPCs, and the marked downregulation of KITLG-encoding SCF, are predicted to have major consequences for residual normal hematopoiesis. Depletion of Scf from Lepr⁺ BMSCs in mice results in cytopenia, hypocellularity of the marrow, and a reduction in HSC number (5), indicating that stromal SCF is critical for the maintenance of normal HSCs and hematopoiesis. Lepr⁺ BMSC-derived Scf is also critical for the maintenance of kit⁺-restricted hematopoietic progenitors (in particular CLPs and MEPs and to a lesser extent CMP and GMP; ref. 9).

Ligand–receptor analyses confirmed that BMSC–HSPC interactions deemed critical for HSPC maintenance, survival, and proliferation (such as KITLG–KIT and CXCL12–CXCR4) were attenuated in nonleukemic HSPC subsets in AML in comparison with their counterparts in NBM (Supplementary Fig. S5A). It is relevant to note that these nonleukemic HSPCs in the AML marrow comprise both healthy, unmutated, HSPCs and “preleukemic” HSPCs carrying founder mutations in genes such as DNMT3A and TET2 (Supplementary Table S1), which cannot be reliably detected by 10× scRNA-seq. We thus evaluated the cellular changes in the residual normal/preleukemic HSPC compartment associated with the predicted disruption of niche signaling. The HSPC fraction in AML displayed a relative reduction of erythroid progenitors (17.41% ± 9.77% vs. 6.31% ± 5.74% of the HSPC fraction in NBM and AML, respectively, P = 0.1) and GMP/monoblasts (22.08% ± 3.29% vs. 12.14% ± 7.15%, P = 0.02), as well as near loss of the CLP population (8.85% ± 3.79% vs. 0.15% ± 0.378%, P = 0.019; Fig. 4B), recapitulating observations in mice with depletion of Scf from Lepr⁺ BMSCs (9).

In line with the notion that residual normal HSPCs in AML may be affected by the loss of stromal niche factors, and in particular, SCF propagating survival and proliferation signals, which occur via activation of PI3K–AKT–MTORC signaling (50, 51), GSEA and Hallmark pathway analyses demonstrated enrichment of inflammatory and apoptotic signatures in residual normal/preleukemic subsets in comparison with their counterparts in the human NBM (Hallmark: TNFα via NFκB; Hallmark: apoptosis and Hallmark: P53 pathway) with concomitant depression of signatures indicative of MTORC signaling, active cellular metabolism, and proliferation (Hallmark: TORC1 signaling; Hallmark: E2F targets; Hallmark: G₂–M checkpoint, Hallmark: Mitotic spindle and Hallmark: Oxidative phosphorylation; Fig. 4C and D). Interestingly, these transcriptional programs indicative of cellular stress were activated most extensively in residual normal/preleukemic cells within the HSC/MPP, MEP, and erythroid progenitor populations (Fig. 4D), the progenitor cells known to be most repressed in AML (resulting in anemia and thrombocytopenia).

The CD34^high HSC compartment was largely of nonleukemic origin (with only 3.66% ± 2.88% of cells within this HSC/MPP population detected positive for NPM1 mutation; Fig. 4A), in line with the existing notion that NPM1m AML has low expression of CD34 and may find its origin in transformation of committed progenitors (52). Reclustering of the HSC population showed that, strikingly, the ST, high-output, HSC/MPP subsets (HSC-1–4) were almost completely depleted from the residual normal/preleukemic HSC pool in AML (74.78% ± 1.29% vs. 16.47% ± 16.12% in NBM and AML, respectively, P = 0.005) with relative conservation of the LT, low-output, HSC-0 subset (Fig. 4E). In line with this observation, the gene signatures indicative of low-output HSCs and quiescence were significantly enriched in the residual normal/preleukemic HSCs (Fig. 4F), whereas the high-output HSC marker genes such as FLT3, STMN1, PLAC8, and MPO were downregulated (Fig. 4G; Supplementary Table S5). The remodeling of the HSC compartment was also found in a patient (X320) in which no founder mutations were detected (Supplementary Table S1), suggesting that the relative loss of “high-output” HSCs is not caused (solely) by the presence of founder mutations such as DNMT3A and TET2 in these cells.

The finding that this specific HSC subset, predicted to be the LT-repopulating subset of HSCs resistant to chemotherapy in mice, resists inflammatory stress seems consistent with the long-standing observation that the human BM can reconstitute normal hematopoiesis after chemotherapeutic eradication of leukemic cells.

Importantly, NPM1m cells within the HSPC fraction, predominantly found within the LMPP and GMP subsets (Fig. 4B), showed preservation of transcriptional MTORC1 and survival signaling (Fig. 4H), suggesting relative resistance to the factors driving the cellular stress in comparison with residual normal HSPCs. Moreover, several interactions were predicted between inflammatory molecules encoded by genes differentially expressed in the “inflamed” BMSC subset (BMSC-2) and the LMPP-like, GMP-like and/or monocyte-like NPM1m AML cells, including HGF–CD44, IL6–IL6 receptor, and JAG1–NOTCH1/NOTCH2 interactions (Supplementary Fig. S5A), which have been reported to drive the initiation, proliferation, and/or chemoresistance of AML (17, 53, 54).

Collectively, the scRNA-seq data, in conjunction with previously established relevance of niche factors for normal hematopoiesis in mice, support the view that niche remodeling by inflammatory activation and suppression of HSPC maintenance factors in AML suppress normal hematopoiesis, whereas LT-HSCs and NPM1m leukemic cells are relatively resistant to the deprivation of supportive signaling from stromal niches. In addition, inflamed BMSCs express factors previously associated with leukemia propagation. Niche remodeling is thus predicted to be a driving force in the competitive advantage of mutated cells over their normal and preleukemic counterparts in leukemogenesis.

Additionally, LEPR⁺ BMSC transcriptional remodeling in AML is predicted to affect cellular signaling to lymphoid immune subsets (CD4 T, CD8 T, NK, B, and plasma cells) in the leukemic bone marrow (Supplementary Fig. S5B). This complex cellular taxonomy of the intercellular signaling defining human AML is anticipated to provide a resource to instruct future experimental interrogation of the relevance of these interactions.

The scRNA-seq analyses suggested that the disruption of BMSC architecture in AML may reflect inflammatory activation of LEPR⁺ BMSCs. To provide experimental support for this view, we tested whether secreted inflammatory factors could induce the inflammatory alterations observed in the BMSC compartment in AML. The significantly enriched gene signature “TNFα signaling via NFκB” (Fig 3C) suggested that TNFα may be one such factor. TNFα levels have earlier been demonstrated to be increased in the plasma of patients with AML (55) and TNF is overexpressed in AML at the transcriptional level (Supplementary Fig. S6A). The scRNA-seq data showed that TNF is highly expressed in adaptive and innate immune cells (CD8 T cells and NK cells) and myeloid lineage cells (GMP/monoblasts and CD14⁺ monocyte populations) in NPM1m AML (Supplementary Fig. S6B and S6C), whereas its canonical receptor TNFRSF1A is predominantly expressed in LEPR⁺ BMSCs (Supplementary Fig. S6D), predicting that LEPR⁺ BMSCs could be one of the cell types in AML most affected by TNFα signaling.

To test whether TNFα can induce inflammatory remodeling of LEPR⁺ BMSCs and attenuate normal hematopoiesis, we performed an in vivo experiment, injecting C57BL/6 mice intraperitoneally with recombinant TNFα (Fig. 5A). This resulted in cytopenia (anemia and thrombocytopenia; Fig. 5B) and a significant reduction of the frequency of LEPR⁺ stromal cells within the nonendothelial (CD45⁻Ter119⁻CD31⁻CD144⁻CD51⁺Sca1⁻) niche (Fig. 5C and D). RNA-seq confirmed high levels of Kitlg and Cxcl12 expression, specifically in the LEPR⁺ stromal cell population, in line with the notion that these cells represent HSPC niches (Fig. 5E; ref. 5). TNFα exposure induced overexpression of inflammatory markers, including Cd44 and Cxcl2, with concomitant downregulation of Kitlg and Cxcl12 (Fig. 5E), in line with the notion that TNFα leads to stromal inflammatory remodeling with subsequent reduction of HSPC niche cells, recapitulating findings in human AML. This was further associated with a reduction in BM cellularity and depletion of distinct subsets of HSPCs, in particular MPPs (Lin⁻cKIT⁺Sca1⁺CD48⁻CD150⁻), CMPs (Lin⁻cKIT⁺Sca1⁻CD34⁺CD16⁻), GMPs (Lin⁻cKIT⁺Sca1⁻CD34⁺CD16⁺), and MEPs (Lin⁻cKIT⁺Sca1⁻CD34⁻CD16⁻; Fig. 5F), resembling the reduction in MPPs and committed progenitor fractions in human NPM1m AML.

To begin testing the functional significance of inflammatory BMSC niche remodeling on normal versus leukemic hematopoiesis using an ex vivo coculture system, we exposed human HS-5 stromal cells to TNFα (10 ng/mL) for 24 hours to trigger the inflammatory response, followed by coculturing with either human CD34⁺ HSPCs or primary NPM1m AML cells for 72 hours (Supplementary Fig. S6E). This resulted in a significant reduction of immunophenotypic HSPCs (CD45⁺Lin⁻CD34⁺) and CFU-GEMM (Supplementary Fig. S6E), whereas the number of NPM1m AML cells was not affected by the inflammatory activation of stromal cells induced by TNFα (Supplementary Fig. S6E). In line with these in vitro findings, in vivo administration of TNFα in a well-established MLL-AF9 mouse transplantation model of AML did not result in a reduction of leukemic cells (Supplementary Fig. S6F), while it resulted in a marked decrease in the frequency of LEPR⁺ BMSCs (Supplementary Fig. S6G). The reduction in the frequency of stromal niches in the MLL-TNFα condition was associated with worsening of anemia and a (nonsignificant) reduction in the number of residual normal MPPs (Supplementary Fig. S6H), recapitulating aspects of the negative effect of TNFα on normal hematopoiesis (Fig. 5).

Together, the data fit a model in which overexpression of TNFα, perhaps in part by activation of the innate and adaptive immune system in NPM1m AML, results in suppression of hematopoiesis with relative resistance of clonal leukemic cells. This is associated with inflammatory remodeling of stromal niches predicted to affect the survival and proliferation of HSPC subsets, although direct effects of TNFα on HSPCs have also been demonstrated (56) and cannot be excluded.

We next asked the question of whether the inflammatory remodeling of BMSCs in AML was restricted to NPM1m cases, or rather a biological commonality in AML across different, distinct, genetic subtypes. To answer this question, bulk RNA-seq was conducted on highly purified CD45⁻CD71⁻CD235a⁻CD31⁻CD271⁺ BMSCs isolated from a cohort of 62 newly diagnosed patients with AML (Supplementary Fig. S7A), uniformly treated within an intensive chemotherapy clinical trial (57) and selected to represent the mutational landscape of AML (Supplementary Table S1). The purity of the sorted stromal population was confirmed by excluding the expression of hematopoietic transcripts (including CD45 (PTPRC), CD34, MPO, and GYPA; Supplementary Fig. S7B) and expression of canonical stromal markers (CD271 (NGFR), COL1A1, CD90 (THY1), and PRRX1 (Supplementary Fig. S7B).

Overall, 1,579 genes were differentially expressed in the AML BMSCs in comparison with age-matched normal controls (n = 8; Supplementary Table S6), including 1,413 upregulated and 166 downregulated genes (Supplementary Table S6). GSEAs revealed differentially expressed gene sets, very similar to those identified by scRNA-seq in the NPM1m subset, including enrichment of the signatures indicative of inflammatory activation (TNFα signaling via NFκB; inflammatory response) and hypoxia, as well as depletion of gene sets indicative of adipogenesis and oxidative phos­phorylation; Fig. 6A). Consistent with this, expression of the NFκB-related inflammatory genes/cytokines such as NFKBIA, CD44, CXCL2, CXCL3, CXCL8, CCL2, LIF, and PTGS2 was significantly upregulated in AML and adipogenesis markers (LPL, ADIPOQ and CD36) significantly downregulated (Fig. 6B; Supplementary Table S6). Similarly, the expression of LEPR and genes encoding the HSPC regulatory factors CXCL12, KITLG, and ANGPT1 was significantly reduced in AML BMSCs (Fig. 6B), in line with the scRNA-seq data from NPM1m patients. Notably, inflammatory activation of BMSC in AML (as represented by the TNFα via NFκB score) was strongly associated with the reduction in expression of HSPC niche genes such as LEPR, CXCL12, KITLG, and ANGPT1 (Supplementary Fig. S7C), congruent with the notion that inflammatory activation results in reduction of HSPC niches and downregulation of HSPC factors, as established in our in vivo experiments (Fig. 5C–E).

To obtain insight into the heterogeneity of these signatures among genetic subtypes, scores were allocated to individual samples and categorized among genetic subtypes using gene set variation analysis (GSVA; ref. 58; Fig. 6C), revealing that inflammatory scores were relatively high in the NPM1m patients and tended to be lower in patients carrying either a RUNX1, ASXL1, or TP53 mutation, although heterogeneity existed within these genetic subsets (Fig. 6C). Correlation of BMSC signatures to genetic risk categories as defined by the ELN (ref. 59; shown to be of prognostic relevance in our cohort of 62 patients, as expected; Supplementary Fig. S7D), revealed that patients in the favorable risk category displayed the most pronounced NFκB inflammatory activation and suppression of niche supportive factors (LEPR, KITLG, and ANGPT1; Fig. 6D), whereas BMSCs in the adverse risk group tended to be less perturbed in comparison with NBM.

Together, the data demonstrate that the remodeling of BMSCs, revealed by scRNA-seq in NPM1m patients and defined by activation of inflammatory signaling and reduced expression of HSPC maintenance factors associated with the loss of the HSPC niche subset, is a biological commonality in AML and that the extent of this remodeling may vary within and between genetically defined subsets of patients.

We next asked the question whether inflammatory remodeling of stromal niches is related to clinical outcomes in AML. The loss of stromal HSPC niche cells has recently been associated with better outcomes upon chemotherapeutic treatment in AML in a mouse model (16). In the MLL-AF9 AML model, genetic depletion of stromal niche cells resulted in delayed relapse after cytarabine treatment (16). The finding that the HSPC niche (BMSC-0) subset of stromal cells is depleted in AML (to varying degrees) as a result of inflammatory remodeling prompted us to interrogate the prognostic value of a transcriptional signature reflecting BMSC remodeling in the context of intensive chemotherapy in human AML.

In order to analyze the correlation between BMSC heterogeneity with clinical outcome in AML, we made a “by-proxy” assessment of BMSC heterogeneity from bulk transcriptomes of all 62 patients using CIBERSORTx, a machine learning method to impute gene-expression profiles and provide an estimation of the relative abundance of cellular subsets in a mixed population (60). The BMSC scRNA-seq dataset was used as a reference, and the relative abundance of all four BMSC subsets was computationally retrieved (Fig. 7A).

In line with our scRNA-seq observations in NPM1m patients, the predicted size of the BMSC-0 subset was significantly reduced in this larger AML population, with a concomitant increase in the size of the inflamed cluster (BMSC-2; Fig. 7A). Interestingly, stratifying patients with AML into a “BMSC-0 preserved” group and a “BMSC-0 depleted” group based on the median level revealed that depletion of the BMSC-0 niche subset was significantly associated with better overall survival (5-year OS 57.2% vs. 31.0% of patients; P = 0.041) and reduced risk of relapse (32.0% vs. 65.2% of patients; P = 0.032; Fig. 7B).

Next, we sought to refine this analysis and get better insight into the genes in the remodeled stromal cells that shape this association between niche remodeling and clinical outcomes in AML. To this end, we generated a list of genes significantly differentially expressed in BMSCs in human AML (in comparison with BMSCs in NBM) that intersected with inflammation-associated genes (differentially expressed in AML BMSC cluster 2 vs. other clusters in AML) and HSPC niche genes (differentially expressed genes in NBM BMSC cluster 0 vs. other clusters in NBM), resulting in a list of 189 genes (Fig. 7C; Supplementary Table S7).

By constructing a penalized multivariable Cox regression model through nested cross-validation (61), 13 of 189 genes were identified to be correlated to OS (Fig. 7C and D). The strongest negative association (coefficient: 0.34; overexpressed in the HSPC niche population and related to poor prognosis) was found for KITLG (Fig. 7D), indicative of a positive correlation between the loss of the HSPC niche population and favorable outcome. Other genes overexpressed in the HSPC niche population (BMSC-0) or inflammatory population (BMSC-2) and associated with poor OS include ARPC5L, encoding an actin-related protein involved in cell migration, MEDAG, a positive regulator for adipocyte differentiation, and PTGDS, encoding a glutathione-independent prostaglandin D synthase that catalyzes the conversion of prostaglandin H2 (PGH2) to prostaglandin D2 (PGD2; Fig. 7D). Upregulation of the inflamed cluster (BMSC-2) genes ZBTB21, SPAG9, ELL2, and SMAD7 (a TGFβ inhibitor positively regulated by inflammatory cytokines; refs. 62, 63), on the other hand, were strongly associated with a favorable outcome (Fig. 7D).

Based on this 13-gene signature, a score was calculated that allowed for stratifying patients with AML into a “BMSC-inflammatory remodeling” group (score^low) and a “BMSC niche-preserved” group (score^high) based on the median level (Fig. 7D). The low score, “BMSC-inflammatory remodeling” group had significantly better OS (5-year OS 80% vs. 8.6%; P < 0.0001) and lower relapse probability (5-year relapse probability 13% vs. 78.4%; P < 0.0001; Fig. 7E). Importantly, this significant survival benefit was found throughout all distinct ELN2017 genetic risk categories (Supplementary Fig. S8A).

Finally, we sought to confirm our finding that stromal niche factors, associated with inflammatory remodeling, are associated with outcomes in an independent cohort of patients with AML. Although many of the factors in the prognostic stromal niche gene signature were either not specifically expressed in stromal cells or not detected in published transcriptional datasets of AML, we found the gene KITLG (i) to be an important determinant of the prognostic value of the stromal prognostic gene signature (Fig. 7D and F), (ii) to be specifically expressed in BMSCs in the AML taxonomy (Supplementary Fig. S8B), and (iii) to be detectable in publicly available transcriptional datasets generated from BM aspirates of patients with AML (18–65 years, treated with intensive chemotherapy), namely, the TCGA-LAML (64) and Bohlander (GSE37642; ref. 65) cohorts, making it the ideal candidate gene to confirm our findings in independent cohorts. In these datasets, stratifying patients with AML into a “KITLG-low” group and a “KITLG-high” group, based on the 75th percentile expression (to clearly discriminate patients with high levels of expression; Supplementary Fig. S8C) confirmed a significantly better OS (5-years OS 44.0% vs. 0%, P = 0.017 in TCGA-LAML; and 39.8% vs. 20.8%, P = 0.034 in Bohlander AML; Fig. 7G), confirming findings in the “training” cohort (Fig. 7E).

Taken together, the data support a predictive working model in which inflammatory remodeling of LEPR⁺ stromal niche cells suppresses normal hematopoiesis (via downregulation of HSPC regulatory factors) but in which leukemia relapse-initiating cells may remain critically dependent on residual levels of stromal niche support (in particular SCF/KITLG-cKIT signaling) for their survival in the context of chemotherapeutic treatment (Supplementary Fig. S9A). Strong inflammation-associated loss of stromal niches may take away this support, resulting in the impediment of leukemia (initiating) cells and reduced risk of relapse (Supplementary Fig. S9B). Alternatively, but not mutually exclusive to this, direct engagement of leukemia-relapse-initiating cells by an activated immune system (driving inflammatory alterations in stromal cells) may explain the association between inflammation and reduced relapse risk in these patients. Future experiments are warranted to test these working models based on our data.

Tumor-promoting inflammation is considered an enabling characteristic of tumorigenesis via tumor-promoting effects that immune cells have on neoplastic cells and disease progression (1). The exact mechanisms linking inflammation to oncogenesis, however, remain incompletely understood.

This, to our knowledge, is the first study to examine the predicted interactions between residual tissue-resident stem/progenitor cells versus their neoplastic counterparts within their native niches and immune environment in human cancer at cellular resolution. By exploiting insights in the hematopoietic system, in which stem cells and their niches have been defined at near cellular resolution, we provide experimental support for the view that tumor-associated inflammation can result in the remodeling and deterioration of innate stem/progenitor cell niches. This is predicted to result in the loss of their capacity to support residual normal HSPCs with relative resistance of neoplastic cells. Additionally, cytokines (e.g., TGFβ1, IL6, and JAG1) and extracellular matrix components produced by the inflamed niches may directly restrict nonleukemic HSPC growth (66) while promoting leukemic cell progression (17, 67), thus providing a conceptual basis for tissue repression and competitive advantage of neoplastic cells in AML (Supplementary Fig. S9A).

The data provide human disease relevance to concepts previously postulated by experiments in various murine models, in which leukemic cells alter stromal niches in ways that were proposed to inhibit normal hematopoiesis (by suppression of key HSPC factors; ref. 24) or secretion of inhibitory factors (22, 23), but at the same time remain dependent on these niches for their survival under chemotherapeutic conditions (16, 68).

Our human data implicate inflammatory signaling as a key driver of niche deterioration in AML with potential therapeutic significance.

It is important, however, to note that the cellular taxonomies we established do not include all niche cells that constitute the mammalian hematopoietic system. Specifically, ECs and stromal niche subsets with close anatomic relationship to the (trabecular) bone, as identified in the murine bone marrow (24, 69, 70), are likely absent or underrepresented in human BM aspirates. The limited availability of core biopsies (not a universal diagnostic procedure in AML), as well as the low amount of available tissue and heterogeneity (sampling bias) of core biopsies preclude robust analyses of these types of niche cells presently.

The stromal cells we retrieved from human aspirates show transcriptional resemblance to the LEPR⁺ perivascular stromal HSPC niches identified in mice (5, 71). This subset of stromal cells, thought to largely overlap with so-called CXCL12-abundant reticular (CAR) cells, is considered to be the largest subset of stromal cells in the mammalian BM (2) and has previously been demonstrated to be instrumental for the maintenance of normal HSCs and committed progenitors (72).

In addition to inflammatory remodeling of LEPR⁺ stromal niches, we describe transcriptional heterogeneity in the human nonleukemic HSC pool and report a previously unanticipated loss of a specific subset of HSCs/MPPs in AML that show strong transcriptional congruence with ST, high-output HSCs previously defined in mice, with relative conservation of LT, low-output HSCs. Due to limitations of 10× scRNA-seq in detecting mutations in genes that are lowly expressed or not near the 3′ end of the transcript (such as DNMT3A or TET2), we were unable to faithfully discern cells carrying “founder” mutations from truly normal (unmutated) cells in the nonleukemic HSC pool. The data, however, are consistent with the recent finding that residual healthy and preleukemic cells in AML are predominantly dormant (73), and, importantly, put this finding in the context of normal hematopoiesis, revealing a loss of the transcriptional high-output subset of HSCs in AML. These findings may provide not only a cellular basis for the BM failure characterizing AML but also the long-standing observation that normal hematopoiesis is typically reconstituted after chemotherapeutic eradication of leukemic cells, which indicated that a population of cells with HSC characteristics must be able to survive these conditions. Our data indicate that LT-HSCs in the human BM are relatively resistant to inflammatory stress and the loss of SCF from stromal niches, perhaps through the expression of inflammation-inhibiting pathways such as the NRA family and TNFAIP3 (36–38). Other mechanisms likely play a role, and it may be that the transition from LT-HSCs to ST-HSCs is specifically impaired in human AML by yet-to-be-identified factors in the leukemic BM. The identification of transcriptional human HSC subsets and elucidation of the cellular taxonomy of AML reported here is anticipated to enable and facilitate future investigations in these directions.

The nature of the drivers of stromal inflammation remains largely speculative at this point. Our data demonstrate that TNFα, which is commonly overexpressed in myeloid malignancies, is able to induce inflammatory activation and remodeling of the stromal microenvironment. The AML taxonomy indicated that activated immune cells may be a source of TNFα in AML, reminiscent of recent work in mice demonstrating that, in viral infection, activated BM-resident CD8⁺ T cells induce damage to stromal niches associated with activation of inflammatory pathways and reduction of expression of Cxcl12 and Scf (74). It is, however, reasonable to assume that other sources of TNFα may exist and that multiple inflammatory cytokines are overexpressed in the AML environment that may contribute to stromal remodeling. Future investigations, instructed by the current taxonomies, may shed further light on this. It is further important to note that the negative effects of TNFα on residual normal hematopoiesis in AML are unlikely to be entirely stromal niche-dependent, as TNFα has been demonstrated to exert both positive and negative direct effects on HSPCs (56, 75–78).

Inflammatory remodeling of stromal niches and the associated loss of niche support was associated with a favorable outcome upon chemotherapeutic treatment, supporting the notion that it may pose a therapeutic vulnerability for leukemia-relapse-initiating cells that may remain critically dependent on some residual niche support to ascertain their survival during chemotherapeutic treatment (Supplementary Fig. S9B). Alternatively, or in addition, inflammatory responses driven by activated innate and adaptive immune cells may contribute to the eradication of leukemia-relapse-initiating cells following chemotherapy in niche-independent fashions. This seems consistent with a recent study demonstrating that inflammatory programs in hematopoietic cells result in suppression of adaptive immunity (including T-cell subsets) and poor prognosis (79). Our data indicate that inflammatory activation of stromal niches may be an additional mechanism of T-cell activation-mediated favorable prognosis. Indeed, the observations that NPM1-mutated AML is associated with a low (hematopoietic) inflammation score (and associated activation of T cells; ref. 79), but a relatively high stromal inflammation score in our study (relative to “poor risk” genetic subtypes), would be consistent with this. In light of these data, the importance of clearly defining the cellular source and molecular makeup of signatures when using the relatively undescriptive term “inflammation” should be stressed.

The finding that niche signatures associated with clinical outcomes in patients with AML treated with intensive chemotherapy may have an impact on risk stratification and therapeutic decision-making in AML, in which current prediction models are instructed by parameters from hematopoietic cells, rather than their stromal environment (59).

Finally, the presented data are anticipated to provide an important resource of human BM signaling, and in particular stem cell niche interactions, defining normal and leukemic hematopoiesis, to serve as a platform for discovery and validation of findings from nonhuman model systems.

NBM samples were obtained by hip bone aspiration from healthy donors for allogeneic transplantation.

AML BM aspirates were obtained from newly diagnosed patients with AML (age 18–65) included in the HOVON-132 clinical trial testing the addition of lenalidomide to intensive treatment in younger and middle-aged adults with newly diagnosed AML (57). Patients received two cycles of intensive induction chemotherapy. Cycle 1 included idarubicin at 12 mg/m² (3-hour infusion on days 1, 2, and 3) and cytarabine at a dose of 200 mg/m² (per continuous infusion on days 1–7) with or without lenalidomide. Cycle 2 contained daunorubicin 60 mg/m² per 1-hour infusion on days 1, 3, and 5 plus cytarabine 1,000 mg/m² given intravenously for 3 hours twice per day on days 1 to 6 with or without the addition of lenalidomide. Patients in CR or CR with incomplete hematologic recovery (CRi) after cycle 2 received consolidation with 1 final additional cycle of intensive chemotherapy with mitoxantrone–etoposide (cycle 3), autologous stem cell transplantation (auto-SCT), or allogeneic stem cell transplantation (allo-SCT). The study was approved by the ethics committees of the participating institutions and was conducted in accordance with the Declaration of Helsinki. All patients gave their written informed consent (57). BM specimens were collected by hip bone aspiration at diagnosis.

Mononuclear cell fractions of human BM aspirates were isolated using lymphoprep and viably frozen in PBS supplemented with 40% heat-inactivated fetal calf serum (FCS; Corning) and 10% dimethyl sulfoxide (DMSO; Sigma-Aldrich). Marrow aspirates for scRNA-seq were cryopreserved within 24 hours after collection. All specimens were collected with informed consent, in accordance with the Declaration of Helsinki.

Viably frozen BM aspirates (mononuclear cell fractions) were thawed in a water bath at 37°C and washed with warm Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% FCS as described in the 10X Genomics protocol “Fresh Frozen Human Peripheral Blood.’

For isolation of hematopoietic fractions, 10% of the thawed BM cells (>20 × 10⁶) from each individual were stained for sorting in PBS supplemented with 0.5% FCS at 4°C with the following antibodies: CD45-APC (1:20, clone 2D1; eBioscience), CD34-AF700 (1:50, clone 581; BioLegend), CD117-PE-CF594 (1:50, clone YB5.B8; BD Biosciences), CD33-PE (1:50, clone P67.6; BD Biosciences), CD3-PE-Cy7 (1:50, clone SK7; BioLegend), CD19-APC-Cy7 (1:50, clone HIB19; BioLegend), and CD38-FITC (1:50, clone HIT2; Life Technologies). For the exclusion of dead cells, 7AAD (1:100; Beckman Coulter) was used. The 7AAD⁻CD45⁺CD34⁺ HSPC fraction, 7AAD⁻CD45⁺CD34⁻CD117⁺ and 7AAD⁻CD45⁺CD34⁻CD33⁺ myeloid fraction, and 7AAD⁻CD45⁺CD34⁻CD117⁻CD33⁻ nonmyeloid (lymphoid) fraction were sorted in DMEM supplemented with 10% FCS using a FACSAria III (BD Biosciences) and BD FACSDiva version 5.0.

The nonhematopoietic fraction was sorted from the thawed BM cells as described before (80). 90% of the thawed cells (>200 × 10⁶) were stained with biotinylated antibodies against CD45 (1:50, clone HI30; BioLegend) and CD235a (1:50, clone HIR2, BioLegend) followed by depletion using magnetic antibiotin beads (20 μL per 10⁷ cells; Miltenyi Biotec) in PBS supplemented with 2% FCS and iMag (BD Biosciences). After depletion, the remaining cells were stained for sorting in PBS containing 0.5% FCS at 4°C with the following antibodies: streptavidin-AF488 (1:100, Invitrogen), CD45-BV510 (1:50, clone HI30; BioLegend), CD235a-PE-Cy7 (1:50, clone HI264; BioLegend), CD71-AF700 (1:20, clone MEM-75; Exbio), CD271-PE (1:50, clone ME20.4; BioLegend), CD105-APC (1:50, clone SN6; eBioscience), CD31-APC-Cy7 (1:20, clone WM59; BioLegend), CD34-eFluor610 (1:50, cone 4H11; eBioscience), and CD144-V450 (1:50, clone 55-7H1; BD Biosciences). 7AAD (1:100; Beckman Coulter) was used for dead cell exclusion. FACSAria III (BD Biosciences) and BD FACSDiva version 5.0 were applied for sorting. 7AAD⁻CD45⁻CD235A⁻CD71⁻ nonhematopoietic fraction was sorted in DMEM supplemented with 10% FCS for scRNA-seq and 7AAD⁻CD45⁻CD235A⁻CD71⁻CD31⁻CD274⁺ BMSCs were sorted in TRIzol (Life Technologies) for bulk RNA-seq.

Prior to scRNA-seq, two fractions sorted separately from the same donors were pooled together (nonhematopoietic fraction with lymphoid fraction and myeloid fraction with HSPC fraction) followed by encapsulation for barcoding and cDNA synthesis using the Chromium Single-cell 3′ Reagent kit v3 (10X Genomics). 3′ gene-expression library was constructed according to the manufacturer's recommendations. The quality and quantity of libraries were determined using an Agilent 2100 Bio-analyzer with 2100 Expert version B.02.11.SI811 software and a High-Sensitivity DNA kit. Libraries were sequenced on a NovaSeq 6000 platform (Illumina), paired-end mode, at a sequencing depth of around 45,000 reads per cell, followed by computational alignment using CellRanger (version 3.0.2, 10X Genomics).

Seurat (R package, version 4.0.0; ref. 81) was used for data preprocessing and downstream analysis. For data preprocessing, datasets were first subjected to quality control steps that included removing doublets (a high ratio of RNA counts vs. feature numbers; >6) and filtering out apoptotic cells determined by the high transcriptional output of mitochondrial genes (>5% of total). Subsequently, single-cell data from separate runs of the same donor were merged to generate a complete picture that includes all cell types for each individual. Aiming at generating a tissue map robustly representing the BM taxonomy in healthy individuals and patients with AML, the merged datasets of each individual were integrated using the integration function in Seurat followed by linear dimensional reduction. That included scaling of gene expression across all cells, principal component analysis on the most variable genes (k  =  2,000) and unsupervised clustering using a shared nearest-neighbor modularity optimization-based clustering algorithm (resolution 0.3–1).). To mitigate the effects of cell-cycle heterogeneity, cell-cycle phase scores were calculated based on canonical markers and were regressed as described in the Seurat protocol. The data were visualized using uniform manifold approximation and projection for dimension reduction (UMAP; ref. 82). Cell types were identified by the Clustifyr R package (25) and by reviewing the expression of canonical markers associated with particular cell types. In order to reveal BMSC heterogeneity at a high resolution, BMSCs were separated using the Cellselector function in Seurat and independently preprocessed/analyzed. For HSC analysis, the HSC/MPP population was reclustered in order to study HSC heterogeneity robustly in single samples. Only the samples with >10 cells in the HSC/MPP population were used for the analysis.

Cell trajectory analysis was performed using the Monocle3 R package (https://github.com/cole-trapnell-lab/monocle3). Cellchat and CellphoneDB were applied to predict ligand–receptor interaction between cell types (39, 83). For differential expression gene analysis between samples, a DEseq2 (84)-based pseuDE package was developed for the aggregation of single-cell level gene counts and application of differential expression (DE) analysis (software is available at https://github.com/weversMJW/pseude). For GSEA (85) between groups, variable genes were identified using Seurat's Findmarker function and ranked using the formula −sign(log₂foldchange) × log₁₀(P value). The R package fgsea was used for the analysis with a permutation of 1,000 using predefined gene sets from the Molecular Signatures Database (MSigDB 6.2) as input. Gene enrichment scores for individual cells were calculated using Seurat's AddGeneScore function, which calculates the score by counting the average expression levels of provided genes, subtracted by the aggregated expression of randomly selected control genes.

Identification of cells carrying the NPM1 mutation was performed using an in-house–developed tool. In short, alignment results produced by the Cell Ranger pipeline in the form of a BAM file were used as input. Aligned reads, representing sequenced cDNA molecules, were extracted at the mutation position and screened for the NPM1 mutation taking into account the cellular barcode (CB) and unique molecular identifier (UMI), irrespective of whether the cell is of hematopoietic origin. Following UMI-based consensus sequence building, mutant-carrying reads were assigned to the cell of origin based on the CB. Reads without direct detection of the NPM1 mutation, i.e., the known 4-nucleotide insert is absent, were screened by a second approach. The NPM1 mutation is located close to an intron-exon boundary which can result in improper alignment when the mutation is positioned near the extremities of the cDNA molecule, i.e., the 4-nucleotide insert mutation is partially, abnormally, or not introduced (only missense variants) in the aligned read, which in turn complicates detection of the mutation. To prevent this issue, all reads without detected NPM1 mutation in the first screening round were aligned to the reference genome (hg38) and a modified version thereof including the 4-nucleotide insert mutation, in both cases taking into account the local NPM1 splicing pattern (i.e., splicing from the penultimate to terminal exon). Reads aligning best to the reference genome were labeled wild-type and those aligning better to the modified version were labeled mutant. Using UMI-based consensus building, mutant-carrying and wild-type reads identified in the second round of screening were attributed to cells based on the CB. Subsequent to the two screening rounds, each cell was assigned a number of UMI-controlled wild-type and mutant reads covering the mutation position. Cells with at least one mutant read detected were classified as NPM1 mutants. Cells with at least three wild-type reads detected covering the mutation position are classified as likely NPM1 wild-type. All remaining cells that did not pass these criteria were classified as “nonassignable (NA).” Software is available at https://github.com/RemcoHoogenboezem/annotate_bam_statistics_sc.

RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer's instructions, in combination of isolation with GenElute LPA (Sigma-Aldrich). cDNA was prepared using the SMARTer procedure through the SMART-Seq v4 Ultra Low Input RNA Kit (Clontech) for Illumina Sequencing. Quantity and quality of cDNA production were assessed using the Agilent 2100 Bio-analyzer and the High Sensitivity DNA kit. cDNA libraries were generated using the TruSeq Sample Preparation v2 guide (Illumina) and paired-end sequencing on a NovaSeq 6000 (Illumina). Adaptor sequences and polyT tails were trimmed from unprocessed reads using fqtrim version 0.9.7. (http://ccb.jhu.edu/software/fqtrim/). Read counts and transcripts per million (TPM) were determined with Salmon version 1.2.1 (86). Gene count estimates were determined from Salmon output with the tximport R package version 1.16 (87). These gene count estimates were in turn normalized, prefiltered according to standard practices, and used for determining the differential gene expression between groups of interest through the DESeq2 package (version 1.28; ref. 84) using default parameters. GSEA was performed with GSEA software (version 3.0, Broad Institute) using predefined gene sets from the Molecular Signatures Database (MSigDB 6.2). Gene lists were ranked on the basis of the log2 FC made available through the DESeq2 package. Classic enrichment statistics with 1,000 permutations were used to determine significant enrichment within gene sets. GSVA (R package) was applied for gene enrichment analysis in single samples.

BM biopsies were fixed in 4% paraformaldehyde, decalcified in 150 mmol/L EDTA, washed in 70% ethanol, and embedded in paraffin. Sections of 5-μm BM biopsies were deparaffinized in xyleen and hydrated in a graded series of ethanol. Antigen retrieval was achieved by microwave treatment in TRIS-EDTA buffer (1M Tris, 0.5M EDTA, 0.05% Tween-20, PH = 9.0). Slides were subsequently blocked using BloxALL (Vector Laboratories), to block endogenous peroxidases, PBS supplemented with 0.5% Tween-20, 5% human serum, 5% goat serum, and 5% donkey serum to circumvent nonspecific binding, and were stained overnight at 4°C with primary antibodies. The primary antibodies that were used in this study included rabbit anti-human CD271 (1:100, HPA004765, BioLegend), mice anti-human CXCL12 (1:50, clone 79018, Invitrogen), and rat anti-human CD44 (1:50, clone IM7, Invitrogen). The secondary antibodies included Alexa Fluor 488–labeled donkey anti-mouse (Invitrogen, 1:200), horseradish peroxidase (HRP)-labeled goat anti–rat (1:100, Jackson ImmunoResearch), and cy5-labeled donkey anti-rabbit antibodies (1:800, Jackson ImmunoResearch). A Tyramide superboost kit (Invitrogen) was used to label CD44 in AF555 following the manufacturer's protocol. The stained sections were mounted in ProLong Diamond containing DAPI (Invitrogen). Images were acquired using a Leica SP5 confocal microscope at 40× magnification and were subsequently analyzed using Fiji software.

Frozen BM aspirates were thawed on ice and stained with biotinylated antibodies against CD45 (1:50, clone HI30; BioLegend) and CD235a (1:50, clone HIR2, BioLegend) followed by depletion using magnetic anti-biotin beads (20 μL per 10⁷ cells; Miltenyi Biotec) in PBS supplemented with 2% FCS and iMag (BD Biosciences). After depletion, the remaining cells were stained in PBS containing 0.5% FCS at 4°C with the following antibodies: Streptavidin-AF488 (1:100, Invitrogen), CD45-BV510(1:50, clone HI30; BioLegend), CD235a-PE-Cy7(1:50, clone HI264; BioLegend), CD71-AF700 (1:20, clone MEM-75; Exbio), CD271-PE (1:50, clone ME20.4; BioLegend), CD31-APC-Cy7 (1:20, clone WM59; BioLegend), and CD44-APC (1:50, clone IM7, Sony). 7AAD (1:100; Beckman Coulter) was used for dead cell exclusion. Samples were measured using FACSymphony A5 Cell Analyzer and analyzed with a FlowJo_v10.6.1 program.

Three-month-old C57BL/6 mice were purchased from Charles River and maintained in specific pathogen-free conditions in the Experimental Animal Center of the Erasmus MC (EDC). These mice were intraperitoneally (i.p.) injected with 5 μg recombinant murine TNFα (PeproTech) for 5 consecutive days, followed by peripheral blood collection for complete blood count, and were sacrificed at day 6 by cervical dislocation. For establishment of MLL-AF9 AML mice model, retroviral plasmid containing MLL-AF9-EGFP (a gift from Dr. Stefan Erkeland, ErasmusMC Rotterdam) were transfected into the platinum-E (Plat-E) retroviral packaging cell line using LipoD293 transfection reagent (SignaGen) for retroviral production. The virus was transduced into freshly isolated Lin⁻ BM cells derived from 3-month-old B6.SJL mice (from Charles River) after Lin⁺ cell depletion using mouse Lineage Cell Depletion Kit (Miltenty Biotec). One million transduced cells (containing ±5% GFP⁺ cells) were intravenously injected into lethally irradiated (9.5 Gy) age-matched C57BL/6 mice. When AML had developed (white blood cell counts >15 × 10³/mm³) 30 days after injection, mice were sacrificed and BM cells (%GFP⁺>90%) were collected. For testing the effects of TNFα on AML, 50,000 BM cells derived from the MLL-AF9 AML mice were intravenously injected into nonirradiated 3-month-old C57BL/L mice. After 30 days of transplantation, 5 μg recombinant murine TNFα (PeproTech) was i.p. injected into the mice for 5 consecutive days, followed by peripheral blood collection for complete blood count. Mice were sacrificed on day 11 (counted from the start date of TNFα injection) by cervical dislocation. Animal studies were approved by the Animal Welfare/Ethics Committee of the EDC in accordance with legislation in the Netherlands (approval No. EMC 2067, 2714, 2892, and 3062).

Mice collagenased bone fraction cells were stained with CD45.2-APC-Cy7 (1:200, clone 104, eBiolegend), Ter119-BV510 (1:50, clone TER-119, BioLegend), CD31-PE-Texas-red (1:100, clone MEC13.3, BD Biosciences), CD144-PE-Cy7 (1:200, clone BV13, BioLegend), CD51-PE (1:50, RMV-7, BioLegend), Sca1-PacificBlue (1:100, D7, BioLegend), and LEPR-biotin (1:50, catalog RB01, R&D Systems) followed by staining of streptavidin-APC (1:100, BioLegend) and 7AAD (1:100) in PBS + 0.5% FCS. 7AAD⁻CD45.2⁻Ter119⁻CD31⁻CD114⁻CD51⁺Sac1⁻LEPR⁺ and 7AAD⁻CD45.2⁻Ter119⁻CD31⁻CD114⁻CD51⁺Sca1⁻LEPR⁻ BMSCs were sorted in TRIzol using FACSAria III (BD Biosciences) and BD FACSDiva version 5.0. Data were analyzed using FlowJo software. For HSPC analysis, mice BM mononuclear cells were first stained with lineage (Lin) cocktail containing biotinylated antibodies against Ly-6G/Ly-6C (Gr-1) Biotin RB6-8C5 (BD Biosciences; 1:100), CD11b Biotin M1/70 (BD Biosciences; 1:100), Ter119 Biotin TER-119 (BD Biosciences; 1:100), CD3e Biotin 145-2C11 (BD Biosciences; 1:100), CD4 Biotin GK1.5 (BD Biosciences; 1:100), CD8 Biotin 53-6.7 (BD Biosciences; 1:100) and B220 Biotin RA3-6B2 (BD Biosciences; 1:100). After one washing step, cells were incubated with Streptavidin Pacific Orange (Life Technologies; 1:200), together with a combination of the following antibodies: Sca1 Pacific Blue E13-161.7 (BioLegend; 1:100), CD48 AF700 HM48-1 (BioLegend; 1:100), CD150 PE-Cy7 TC15-12F12.2 (BioLegend; 1:100), CD127 (IL7RA) APC A7R34 (BioLegend; 1:100), cKit PE-CF594 2B8 (BD Biosciences; 1:100), CD16/32 APC-Cy7 2.4G2 (BD Biosciences; 1:100), FLT3 APC A2F10 (1:40), and 7-ADD (1:100).

The STR-authenticated human BMSC cell line HS-5 was purchased from ATCC (CRL 3611) in 2020 and maintained in RPMI-1640 medium (Thermo Fisher Scientific) supplemented with 10% FCS and 1% penicillin and streptomycin. The cells used in the experiments had passage numbers <12. Mycoplasma testing was performed monthly using a Mycoplasma PCR detection kit (Thermo Fisher Scientific). Primary CD34⁺ cells were obtained from umbilical cord blood (Erasmus MC) using a Ficoll gradient protocol and by magnetic-activated cell sorting. Primary AML cells were obtained from patients’ aspirates as described above. For ex vivo coculture, 180,000 HS-5 cells were seeded in each well of a 12-well plate in RPMI-1640 medium supplemented with or without 10 ng/mL recombinant human TNFα (PeproTech). After 24 hours and washing of HS-5 cells with PBS, CB CD34⁺ cells or primary AML cells were seeded on top of the HS-5 cells and covered with GMP serum-free Stem Cell Growth Medium (CellGenix GmbH) supplemented with 50 ng/mL TPO, 50 ng/mL FLT3, 50 ng/mL SCF, and 20 ng/mL IL3 (only for AML cells; all from PeproTech). Cells were cultured for 3 days at 37°C with 5% CO₂.

Counting of CB CD34⁺ cells after 3 days of ex vivo coculture was accurately obtained with flow-count fluorosphere beads (Beckman Coulter), in combination with the following antibodies: Lin-cocktail-FITC (1:100, catalog 22-7778-72, eBioscience), CD45-APC (1:20, clone 2D1; eBioscience), CD34-AF700 (1:50, clone 581; BioLegend). For counting of AML cells, an antibody cocktail containing CD45-APC (1:20, clone 2D1; eBioscience), CD34-AF700 (1:50, clone 581; BioLegend), CD117-PE-CF594 (1:50, clone YB5.B8; BD Biosciences), CD33-PE (1:50, clone P67.6; BD Biosciences) was used in combination with flow-count fluorosphere beads. In both cases, dead cells were excluded based on the DAPI (1:7,500) gate. The data were acquired using an LSRII flow cytometer (BD Biosciences) and analyzed using FlowJo software.

On day 3 of coculture, 2000 mononuclear cells (MNCs) from HS-5-Veh (PBS) condition and in HS-5-TNFα condition were resuspended in IMDM (Thermo Fisher Scientific). This cell suspension was mixed with MethoCult GF H84434 (STEMCELL Technologies), which allows the growth of colonies from all three lineages, and triplicate dishes were plated. The MethoCult plates were kept at 37°C in a 5% CO₂ incubator for 2 weeks until colony counting under a light microscope (Zeiss).

The R package dCVnet (https://github.com/AndrewLawrence/dCVnet) was used to perform Lasso-penalized nested cross-validated Cox regression [k-fold.outer = 5, k-fold.inner = leave-one-out cross-validation (i.e., 1)] to build a predictive model for OS based on the centered and scaled gene expression of the 189 genes from all available AML samples (n = 62). This resulted in the retention of 13 out of 189 genes weighted according to their contribution to the model. Similar to Ng and colleagues (61), we obtained a weighted score by taking the linear combination of the centered and scaled gene-expression levels of the 13 retained genes weighted by the obtained regression coefficients. Score = (ADAMTS4 × −0.066) + (ZBTB21 × −0.079) + (HAS2 × −0.035) + (SPAG9 × −0.104) + (KITLG × 0.345) + (CXCL3 × 0.021) + (MEDAG × 0.224) + (PTGDS × 0.094) + (ARPC5 L × 0.269) + (UHRF1BP1 L × −0.0008) + (SMAD7 × −0.069) + (ANXA2 × 0032) + (ELL2 × −0.072). The median score value was used to dichotomize the AML cohort into low- and high-score groups. These groups were labeled as “BMSC niche-preserved” (low score) and “BMSC niche inflammatory remodeling” (high score), respectively. Distinction into the “favorable,” “intermediate,” and “adverse” groups was based on the ELN2017 genetic risk classification. The log-rank test was used to assess statistical differences between the survival distributions, a P ≤ 0.05 was considered statistically significant. For clinical outcome analysis using publicly available datasets, TCGA-LAML and Bohlander AML (GSE37642) datasets were acquired using R package TCGAbiolinks and GEOquery, respectively. Patients between 18 and 65 years old were selected and stratified into KITLG-low and KITLG-high groups based on the 75% percentile of normalized counts of KITLG. The log-rank test was used to assess statistical differences between the survival distributions, a P ≤ 0.05 was considered statistically significant.

Statistical analysis was performed using Prism 8 (GraphPad Software) and/or R program. Unless otherwise specified, unpaired, two-tailed Student t test (single test), one-way ANOVA (multiple comparisons), or Spearman rho test (correlation analysis) was used to evaluate statistical significance, defined as P-value < 0.05. All results in bar graphs are mean value ± SD.

The RNA-seq data generated in this study were deposited and publicly available in European Genome-Phenome Archive (EGA) at accession number EGAS00001007330. The data analyzed in this study were obtained from the database of Genotypes and Phenotypes (dbGaP) at phs000178 (TCGA-LAML), and Gene-Expression Omnibus at GSE37642. All other data supporting the findings of this study are cited in Methods, in supplementary documents, or are available upon request from the authors.

No disclosures were reported.

L. Chen: Conceptualization, resources, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft. E. Pronk: Formal analysis, methodology. C. van Dijk: Data curation, formal analysis, methodology. Y. Bian: Methodology. J. Feyen: Methodology. T. van Tienhoven: Formal analysis, methodology. M. Yildirim: Investigation, methodology. P. Pisterzi: Methodology. M.M. De Jong: Methodology. A. Bastidas: Methodology. R.M. Hoogenboezem: Data curation, software, methodology. C. Wevers: Software, methodology. E.M. Bindels: Data curation, methodology. B. Löwenberg: Conceptualization, investigation. T. Cupedo: Conceptualization, investigation. M.A. Sanders: Conceptualization, investigation, methodology. M.H.G.P. Raaijmakers: Conceptualization, supervision, funding acquisition, investigation, writing–original draft, project administration.

The authors wish to thank Nathalie Papazian and Michael Vermeulen for technical assistance, Mariette ter Borg for providing CD34⁺ cells; the HOVON/SAKK leukemia working group and all its members and participating sites for conducting the HOVON132 trial; Peter Valk, Patrycja Gradowska, and Jurjen Versluis for providing the genetic and clinical data; the Josephine Nefkens Precision Cancer Treatment Program for infrastructural support and members of the Erasmus MC Department of Hematology for providing scientific discussion and members of the Erasmus MC animal core facility EDC for help with animal care. This work was supported by grants from the Dutch Cancer Society (KWF Kankerbestrijding), Amsterdam, the Netherlands (grant EMCR 10488 and 11092 to M.H.G.P. Raaijmakers).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.