The bone marrow microenvironment at single-cell resolution

Abstract

The bone marrow microenvironment has a key role in regulating haematopoiesis, but its molecular complexity and response to stress are incompletely understood. Here we map the transcriptional landscape of mouse bone marrow vascular, perivascular and osteoblast cell populations at single-cell resolution, both at homeostasis and under conditions of stress-induced haematopoiesis. This analysis revealed previously unappreciated levels of cellular heterogeneity within the bone marrow niche and resolved cellular sources of pro-haematopoietic growth factors, chemokines and membrane-bound ligands. Our studies demonstrate a considerable transcriptional remodelling of niche elements under stress conditions, including an adipocytic skewing of perivascular cells. Among the stress-induced changes, we observed that vascular Notch delta-like ligands (encoded by Dll1 and Dll4) were downregulated. In the absence of vascular Dll4, haematopoietic stem cells prematurely induced a myeloid transcriptional program. These findings refine our understanding of the cellular architecture of the bone marrow niche, reveal a dynamic and heterogeneous molecular landscape that is highly sensitive to stress and illustrate the utility of single-cell transcriptomic data in evaluating the regulation of haematopoiesis by discrete niche populations.

Extended Data Figure 1:

RNA-seq analysis of the bone marrow microenvironment populations. a, Representative two-photon imaging of tdTomato⁺ vascular cells (VE-Cad—tdTomato⁺), perivascular cells (LEPR—tdTomato⁺) and osteoblasts (COL2.3—tdTomato⁺). b, Representative flow cytometry of VE-Cad-tdTomato⁺, LEPR-tdTomato⁺ and COL2.3-tdTomato⁺ populations. c, Principal component analysis of vascular (VE-Cad-tdTomato⁺, n = 4) (red), perivascular LEPR⁺ (LEPR-tdTomato⁺, n = 4) (purple) and osteoblast (COL2.3-tdTomato⁺, n = 4) (blue) populations, based on the expression of the 1,000 most-variable genes in bulk RNA-seq. d, Relative expression levels of COL2.3⁺, LEPR⁺and VE-Cad⁺ signature genes across the three subpopulations of the bone marrow niche in bulk RNA-seq. Normalization and statistical analysis were performed using the DESeq2 R package. e, Normalized expression levels of the population-specific markers VE-Cad (Cdh5), LEPR (Lepr) and COL2.3 (Col1a1) for all scRNA-seq clusters. n = 9,622 cells. The data are mean ± s.e.m. Experiments were repeated independently on more than 10 (a, b) biological samples with similar results.

Extended Data Figure 2:

Analysis of VE-Cad⁺, LEPR⁺ and COL2.3⁺populations. a, Schematic workflow of independent and integrated analysis of VE-Cad⁺, LEPR⁺ and COL2.3⁺ scRNA-seq data. b, t-SNE representation of VE-Cad⁺ populations only. Cluster C1 corresponds to arterial cluster V1 (Ly6a^high). Cluster C2 corresponds to sinusoidal cluster V2 (Stab2^high). Cluster C3 is the cycling cluster. c, Normalized expression of arterial, sinusoidal and cycling markers (Ly6a, Stab2 and Mki67, respectively) (n = 4,551 cells). d, Gene signatures of VE-Cad⁺subpopulations in the bone marrow, based on the relative expression levels of the ten most-significant markers for each cluster. e, t-SNE representation of the LEPR⁺ population only. Cluster C1 corresponds to the adipocytic-primed cluster P1 (Mgp^high) and encompasses cluster P2. Cluster C2 corresponds to P3 (Wif1^high), and C3 to P4 (Spp1^high). f, Normalized expression of P1, P3 and P4 markers (Mgp, Wif1 and Spp1, respectively). n = 3,907 cells. g, Gene signatures of LEPR⁺subpopulations in the bone marrow, based on the relative expression levels of the ten most-significant markers for each cluster. h, t-SNE representation of the COL2.3⁺ population only. Cluster C1 corresponds to cluster O1 (Col16a1^high), cluster C2 to O2 (Fbn1^high) and C3 to O3 (Bglap^high). i, Normalized expression of O1, O2 and O3 markers (Col16a1, Fbn1 and Bglap, respectively). n = 1,114 cells. j, Gene signatures of COL2.3⁺ subpopulations in the bone marrow, based on the relative expression levels of the ten most-significant markers for each cluster. Cluster C4 represents arterial vascular cells (Cdh5, Kdr and Ly6e); C5 is glial-like cells (Fabp7, Mpz and Endrb); and C6 is myocyte-like cells (Pgf, Pln and Acta2). The data shown in c, f, i are mean ± s.e.m. MAST with Bonferroni correction (d, g, j).

Extended Data Figure 3:

Characterization of VE-Cad⁺ subpopulations. a, Relative average scRNA-seq expression levels of the previously described arterial and sinusoidal gene signatures at a steady state, within the VE-Cad⁺ clusters V1 and V2. b, Average scRNA-seq expression levels (left) (n = 4,669 cells) and bone marrow immunofluorescence (right) of arterial expression of SCA-1 (Ly6a), CD102 (Icam2) and PODXL (Podxl) (n = 3 mice); LAMA1 staining (blue) labels all bone vessels; yellow arrowheads indicate arterial vessels. Dashed lines mark bone marrow (bm), compact bone (cb), growth plate (gp) and metaphyseal (mp) bone regions. c, Bone marrow immunofluorescence of arterioles co-stained with SCA-1 and PODXL. n = 3 mice. d, Average scRNA-seq expression levels (left) and bone marrow immunofluorescence (right) of sinusoidal VEGFR3 (Flt4) (red) and CD54 (Icam1) (green) markers. n = 3 mice. LAMA1 staining (blue) labels all bone vessels. e, Average scRNA-seq expression levels and representative flow cytometry analysis of the arterial subpopulation (V1) using SCA-1 and scRNA-seq-identified LY6C (Ly6c1) and CD34 (Cd34) from VE-Cad–tdTomato bone marrow (n = 3 mice). Cells were pre-gated on DAPI⁻tdTomato^high cells. The data in b, d, eare mean ± s.e.m.

Extended Data Figure 4:

Characterization of perivascular LepR⁺ subpopulations. a, Relative average scRNA-seq expression levels of the adipocytic(Hp, Lpl, Adipoq, Slc1a5, Cd302, Gas6 and Apoe) and osteo-associated (Ibsp, Spp1, Alpl, Wif1, Bglap, Sp7 and Runx2)genes, as well as markers used for characterization of LEPR⁺ cell (Esm1, Vcam1, Cd200 and Cd63). b, ESM1 (green) bone marrow immunofluorescence of LEPR-tdTomato femur. LAMA1 staining (blue) labels all bone vessels. Yellow arrowhead, LEPR⁺ESM1⁺ cells; white arrowheads indicate LEPR⁺ESM1⁻ cells. c, CD200 (green) and CD63 (red) bone marrow immunofluorescence of LEPR-tdTomato femur. Nuclei, DAPI (blue). Yellow arrowhead, LEPR⁺, CD200⁺, CD63⁺ cells; white arrowhead, LEPR⁺, CD200⁻, CD63⁻ cells. d Human mesenchymal stem cell (hMSC) gene signature module score, overlaid on t-SNE representation (n=9,622 cells). e, Flow cytometry representation of VCAM1^highCD63^low and VCAM1^lowCd63^high cells of tdTomato⁺ subpopulations in LEPR-tdT bone marrow. Cells were pre-gated on DAPI⁻ tdTomato^high cells. f, Fibroblastic colony-forming unit activity of sorted total LEPR⁺ (purple), LEPR⁺VCAM1^lowCD63^high (maroon), and LEPR⁺VCAM1^highCD63^low (yellow) cells from bone marrow of LEPR-tdTomato mice (n=8).The data are mean ± s.d. (f). N.S., not significant, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001 (Student’ t-test, two-tailed (f). Data are representative of two (b, c) or three (e, f) independent experiments.

Extended Data Figure 5:

Characterization of Col2.3⁺ subpopulations and cycling cells. a, Relative average scRNA-seq expression levels of O1 (Edil3, Mmp14, Ostn, Col121, Angptl2 and Col16a1), O2- (Sox9, Comp, Chad and Col10a1), and O3- (Col11a2, Col1a2, Sparc, Bglap2 and Bglap3) associated genes. b-d, Average scRNA-seq expression levels (n = 9,622 cells) and bone marrow immunofluorescence of MMP14 (green) (b) (n = 3 mice), CD9 (green) (c) (n = 3 mice) and CAR3 (green) (d) (n = 3 mice) in Col2.3-tdTomato femur, with arrows indicating co-staining with tdTomato (red). Nuclei, DAPI (blue). Arrowhead, Col2.3⁺MMP14⁺ (b), COL2.3⁺CD9⁺ (c), and COL2.3⁺CAR3⁺ (d). e, Expression levels of Mki67 in all identified subpopulations (n = 9,622 cells). f, Enriched Gene ontology biological processes terms most-strongly associated with cycling cluster (C), color-coded by significance of enrichment and size on the basis of the fraction of overlapping genes. n = 9,622 cells. g, Contribution of VE-cad-tdTomato⁺, LEPR-tdTomato⁺, and COL2.3-tdTomato⁺ cells to the cycling cluster (C) at steady state n = 70 cells. The data b-e are mean ± s.e.m

Extended Data Figure 6:

Effect of treatment with 5-FU on subsets of the bone marrow niche. a, Representative hematoxylin and eosin-stained sections of BM on day 5 following control or 5-FU treatment (n = 3 animals). b, Frequency and numbers of BM LSK cells on day 5 following control (CNTRL) (n = 4) or 5-FU treatment (n = 5). c, Absolute numbers of BM niche cells, vascular VEcad⁺ (CNTRL: n = 4; 5-FU: n = 10), perivascular LepR⁺ (CNTRL: n = 2; 5-FU: n = 5), and Col2.3⁺ osteoblasts (CNTRL: n = 4; 5-FU: n = 5) from control (PBS) and 5-FU treated animals. d, Gene signatures of niche LepR⁺ sub-populations, including cluster P5, based on the average relative expression levels of the 10 most significant markers for each cluster exclusively within the LepR⁺ subset. MAST with Bonferroni correction. e, Relative expression levels of upregulated adipogenesis-associated genes and downregulated osteogenesis-associated genes in LepR⁺ subpopulations in response to 5-FU treatment. f, LepR⁺ cells enriched pathways in response to 5-FU treatment (n = 17,374 cells). Fisher exact test. g, Contribution of VEcad-tdT⁺, Lepr-tdT⁺, and Col2.3-tdT⁺ cells to cycling cluster (C) following 5-FU treatment (n = 418 cells). h, Expression levels of Mki67 in all identified subpopulations at steady state and following 5-FU (n = 17,374 cells). Data are mean ± s.d. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001 (Student’s t-test, two-tailed) (b-c). The data (h) are mean ± standard error.

Extended Data Figure 7:. Validation of scRNA-seq in VE-Cad⁺ and LEPR⁺cells following treatment.

a, b, Heat map and hierarchical clustering of mean normalized expression values of VE-Cad⁺ (a) and LEPR⁺ (b) control-treated and 5-FU-treated samples for 5-FU-modulated genes (top 50 differentially expressed) in two independent scRNA-seq experiments. c, d, log-transformed fold changes of differentially expressed genes (up- and downregulated > 1.2×, adjusted P value < 0.001) in 5-FU-treated versus control-treated VE-Cad⁺ (n = 697 genes) cells (c) and LEPR⁺ (n = 829 genes) cells (d), for two independent experiments. The trend line (dashed) and the confidence interval (grey shading) were calculated using the linear model. MAST with Bonferroni correction. VE-Cad⁺ control-treated, n = 5,796 cells; VE-Cad⁺ 5-FU-treated, n = 1,481 cells; LEPR⁺ control-treated, n = 6,128 cells; LEPR⁺ 5-FU-treated, n = 4,867 cells.

Extended Data Figure 8:. Analysis of Dll4-mCherry, Dll1-mCherry, and Jag1-mCherry reporter animals.

a, Reverse-transcription PCR of Dll4, Dll1 and Jag1 in vascular cells (red), perivascular LEPR⁺ cells (purple) and osteoblasts (blue), normalized to Gapdh (n = 4 mice). b, Low-magnification immunofluorescence images of thymus sections from Dll4-mCherry, Dll1-mCherry and Jag1-mCherry mice. c, Representative flow cytometry, measuring mCherry fluorescence in total bone marrow of Dll4-mCherry, Dll1-mCherry and JAG1-mCherry mice. Indicated values represent percentages of the complete CD144⁺ and mCherry⁺CD144⁻ populations. Cells were pre-gated on DAPI⁻ cells. d, Representative mCherry levels in DAPI⁻CD45^lowTER119^lowCD144⁺ bone marrow endothelial cells from Dll4-mCherry (red) (n = 4), Dll1-mCherry(blue) (n = 3), Jag1-mCherry (black) (n = 4) and control (grey) (n = 5) mice. e, f, Representative immunofluorescence metaphysis and diaphysis of Dll4-mCherry (e) and Dll1-mCherry (f) bone marrow (n = 3 mice). mCherry (red) and LAMA1 (blue). g, h, Representative two-photon images of bone marrow from intact (left) or dextran-injected (right) Dll4-mCherry(g) and Dll1-mCherry (h) mice (n = 3 mice). i, Normalized counts of key differentially expressed genes from bulk RNA-seq performed on CD144⁻DLL1⁺ cells (purple) (from n = 2 mice) and CD144⁺DLL1⁺ cells (black) (from n = 2 mice). j, Representative flow cytometry histogram measuring mCherry fluorescence in NK1.1⁺ population from Dll1-mCherry (pink) and control (black) mice (n = 3 mice). The data are mean ± s.d. N.S., not significant, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, Student’s t-test, two-tailed. Data are representative of two (a, e–h, j) or three (b–d) independent experiments.

Extended Data Figure 9:. Gene expression program of myeloid differentiation is enhanced in VE-Cad-Dll4^i3COIN HSPCs.

a, b, Representative percentage of bone marrow progenitors in VE-Cad^creER-Dll4^i3COIN and littermate-control mice for common lymphoid progenitors gate (CLP gate) (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 11) (a) and granulocyte–monocyte progenitors gate (GMP gate) (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 13) (b). c, Frequencies of CD3⁺ T cells (control, n = 12; VE-Cad^creER-Dll4^i3COIN, n = 11). d, Total numbers of mature haematopoietic subsets in tamoxifen-treated VE-Cad^creER-Dll4^i3COIN and littermate-control mice, including B220⁺ B cells (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 10), CD3⁺ T cells (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 10) and CD11b⁺GR1⁺ myeloid cells (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 10). e, Absolute numbers of thymocytes from VE-Cad^creER-Dll4^i3COIN mice (n = 11) and littermate-control mice (n = 9). f, Representative flow cytometry analysis of thymic subsets in tamoxifen-treated VE-Cad^creER-Dll4^i3COIN and littermate-control mice. g, Frequencies (control, n = 8; VE-Cad^creER-Dll4^i3COIN, n = 5) and absolute numbers (control, n = 6; VE-Cad^creER-Dll4^i3COIN, n = 4) of early thymic progenitors in thymi from VE-Cad^creER-Dll4^i3COIN and littermate-control mice. h, Percentage (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 10) of HSCs, MPP2 cells, MPP3–4 cells, MPP4 cells and LSK cells from the bone marrow of VE-Cad^creER-Dll4^i3COIN and littermate-control mice. i, Total numbers of bone marrow HSCs, MPP2 cells, MPP3 cells, MPP4 cells and LSK cells from VE-Cad^creER-Dll4^i3COIN and littermate-control mice (control, n = 10; VE-Cad^creER-Dll4^i3COIN, n = 10). j, Representative immunofluorescence of early progenitors (Lin⁻CD48⁻CD150⁺) adjacent to the DLL4-producing vascular endothelium in Dll4-mCherry. Lin cocktail, CD11b, GR1, CD41, TER119 and B220. Arrowhead, Lin⁻CD150⁺ progenitors. n = 3 mice. k, scRNA-seq t-SNE visualization of the LSK compartment (n = 21,116 cells), colour-coded by genotype. l, m, Distribution of enrichment scores for myeloid progenitor (l) and HSC (m) gene signatures within the scRNA-seq-profiled HSPC populations from bone marrow of tamoxifen-treated VE-Cad^creER-Dll4^i3COIN mice (pooled n = 2) and littermate-control mice (pooled n = 2). n = 21,116 cells. The data are mean ± s.d. NS, not significant, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, Student’s t-test, two-tailed (a–d, g–i) or Wilcoxon rank-sum test (l, m). Data are representative of four (a–e) or two (g) independent experiments.

Extended Data Figure 10:. Deletion of endothelial Dll1 does not affect early lineage priming of haematopoietic progenitors.

a–f, Flow cytometry analysis of bone marrow progenitors in VE-Cad-Dll1^fl/fl and littermate-control mice, showing equivalent frequencies of common lymphoid progenitors (control, n = 10; VE-Cad-Dll1^fl/fl, n = 8) (a), granulocyte–monocyte progenitors (control, n = 7; VE-Cad-Dll1^fl/fl, n = 9) (b) and MPP4 cells (control, n = 10; VE-Cad-Dll1^fl/fl, n = 10) (c), of B220⁺ B cells (control, n = 10; VE-Cad-Dll1^fl/fl, n = 10) (d, f), CD3⁺ T cells (control, n = 10; VE-Cad-Dll1^fl/fl, n = 10) (d, f) and CD11b⁺GR1⁺ monocytic–granulocytic subset (control, n = 8; VE-Cad-Dll1^fl/fl, n = 9) (e, f). g, Absolute numbers of thymocytes from VE-Cad-Dll1^fl/fl and littermate-control mice (control, n = 6; VE-Cad-Dll1^fl/fl, n = 10). h, Representative flow cytometry analysis of thymic subsets in VE-Cad-Dll1^fl/fl and littermate-control mice. i, Frequencies (control, n = 4; VE-Cad-Dll1^fl/fl, n = 6) and absolute numbers (control, n = 4; VE-Cad-Dll1^fl/fl, n = 6) of early thymic progenitors from thymi of VE-Cad-Dll1^fl/fl and littermate-control mice. The data are mean ± s.d. N.S., not significant, Student’s t-test, two-tailed. Data are representative of three independent experiments.

Figure 1:. Single-cell RNA-sequencing analysis of the BM microenvironment at steady state.

a, Schematic representation of BM niche: osteoblasts (Col2.3⁺, blue), perivascular LepR⁺ (LepR⁺, purple), and vascular cells (VEcad⁺, red). b. Expression levels of genes associated with niche-specific Cre strains, including Cdh5 (VE-cadherin), Lepr (Leptin receptor) and Col1a1 (Col2.3) overlaid on tSNE representation (n=9,622 cells). Dashed lines generally encompass the cells profiled in Col2.3⁺, LepR⁺, VEcad⁺ libraries. c, Gene signatures of VEcad⁺, LepR⁺, and Col2.3⁺ sub-populations based on the relative expression levels of the 10 most significant markers for each of the 10 clusters, displaying 100 randomly-selected cells per cluster. MAST with Bonferroni correction. d, tSNE visualization of color-coded clustering of BM niche (n=17,374 cells). Dashed lines generally encompass the examined VEcad⁺ (red), LepR⁺ (purple), Col2.3⁺ (blue) BM niche populations and the cycling cells (sea green).

Figure 2:. Expression of pro-hematopoietic factors by the BM microenvironment.

a-c, Gene signatures of BM niche sub-populations based on the average relative expression levels of the 10 most significant markers for each cluster exclusively within the parent niche subsets (VEcad⁺ (a), LepR⁺(b), or Col2.3⁺(c)) at steady state. MAST with Bonferroni correction. d, Reconstructed cell differentiation trajectory of LepR⁺ and Col2.3⁺ clusters colored by subpopulation identity. e, Expression levels of indicated genes in LepR⁺ and Col2.3⁺ cells with respect to their pseudotime coordinates. Black lines depict the LOESS fit of the normalized expression values. f, Expression levels of pro-hematopoietic factors in the non-cycling clusters (n=9,622 cells). The data are mean ± standard error.

Figure 3:. Single-cell transcriptome profiling of the BM niche in response to chemotherapy.

a, VEcad⁺ (red), LepR⁺ (purple), and Col2.3⁺ (blue) cells at steady state (dark) or following 5-FU treatment (light) overlaid on tSNE representation (n=17,374 cells). b, tSNE plot of color-coded clusters of all examined (control and 5-FU treated) BM niche cells (n=17,374 cells). Dashed lines generally encompass the examined VEcad⁺ (red), LepR⁺(purple), Col2.3⁺ (blue) BM niche populations and the cycling cells (sea green). c, Frequencies of each sub-population within VEcad⁺, LepR⁺, and Col2.3⁺ compartments under the indicated condition. d, Relative expression of pro-hematopoietic factors examined in Fig. 2e in response to 5-FU treatment (n=17,374 cells). Adjusted *P≤ 0.05, **P≤ 0.01, ***P< 0.001 (MAST with Bonferroni correction). The data are mean ± standard error.

Figure 4:. Vascular endothelial cells are the major source of Dll4 and Dll1 in the BM.

a, Map of BAC clones used to drive transgenes mDll4-mCherry, mDll1-mCherry, and mJag1-mCherry. b, Representative flow cytometry histogram measuring mCherry fluorescence in DAPI⁻CD45/Ter119^lowCD144⁺ endothelial cells from BM of mDll4-mCherry, mDll1-mCherry, and mJAG1-mCherry. Indicated values represent percentages of the complete CD45/Ter119 population. c-d, Representative two photon images of BM from dextran-injected mDll4-mCherry (c) and mDll1-mCherry (d). (c-d) Scale bars are 100 μm and 25 μm, as indicated in the images. Data are representative of three independent experiments (b-d).

Figure 5:. Dll4 expressed by vascular endothelial cells prevents myeloid skewing of hematopoietic progenitors.

a-c, Representative flow cytometry (a,b) and frequencies (c) of mature hematopoietic subsets in VEc-Dll4^i3COIN and littermate control animals, including B220⁺ B cells (CNTRL: n=12; VEc-Dll4^i3COIN: n=11) and CD11b⁺Gr1⁺ cells (CNTRL: n=12; VEc-Dll4^i3COIN: n=11). d, LSK compartment scRNA-seq tSNE visualization (n=21,116 cells), color-coded by the identified HSPC populations. e, Log-transformed normalized expression levels of myeloid-associated genes (Calr, Ctsg, Elane, and Mpo) and Mki67 across HSPC populations in control and VEc-Dll4^i3COIN LSK cells (n=21,116 cells). Violin plots depict kernel density estimation to show the distribution of expression values. Data are mean ± s.d. *P≤ 0.05, **P≤ 0.01, ***P≤ 0.001 (Student t-test, two-tailed (c), MAST test with Bonferroni correction (e)). Data are representative of four (a-c) independent experiments.