Mapping Distinct Bone Marrow Niche Populations and Their Differentiation Paths

Abstract

The bone marrow microenvironment is composed of heterogeneous cell populations of non-hematopoietic cells with complex phenotypes and undefined trajectories of maturation. Among them, mesenchymal cells maintain the production of stromal, bone, fat, and cartilage cells. Resolving these unique cellular subsets within the bone marrow remains challenging. Here, we used single-cell RNA sequencing of non-hematopoietic bone marrow cells to define specific subpopulations. Furthermore, by combining computational prediction of the cell state hierarchy with the known expression of key transcription factors, we mapped differentiation paths to the osteocyte, chondrocyte, and adipocyte lineages. Finally, we validated our findings using lineage-specific reporter strains and targeted knockdowns. Our analysis reveals differentiation hierarchies for maturing stromal cells, determines key transcription factors along these trajectories, and provides an understanding of the complexity of the bone marrow microenvironment.

Keywords: bone marrow; differentiation; scRNA-seq; stromal cell; transcription factor.

Figure 1.. scRNA-seq Sequencing Reveals Heterogeneous Gene Expression in Bone Marrow Stromal Cells

(A) Schematic of scRNA-seq of non-hematopoietic (CD45⁻/Ter119⁻-), non-endothelial (CD31⁻) mouse bone marrow cells.(B) Heatmap of the most specific significantly enriched genes for each cell cluster.(C) Selected gene sets significantly enriched in the most highly specific genes of each cluster.(D) MA plot for genes significantly differentially expressed (permutation test, false discovery rate [FDR]-corrected p < 0.05) in each cluster versus all other cells. Selected genes of interest are highlighted in black, and all of the significant genes are shown in gray.(E) SPRING plot of single-cell transcriptomes. Each point is one cell, and colors indicate graph-based cluster assignments.

Figure 2.. Characterization of Stromal Subpopulations within the scRNA-seq Data

(A) Heatmap of the five most differentially expressed genes significantly enriched in each cell cluster.(B) Violin plots of previously characterized cluster-specific genes.(C) SPRING plots of stromal cells, colored by expression of the indicated gene. Shown are genes previously shown to play a role in the bone marrow HSC niche.(D) Expression of key lineage specific genes and transcription factors.

Figure 3.. Population Balance Analysis Predicts Early Differentiation Trajectory of Adipocytes, Osteoblasts, and Chondrocytes

(A) Velocyto-calculated RNA velocity vector field overlaid on a SPRING plot.(B) Prediction of start and end cell states using RNA velocity-based Markov process simulation.(C) PBA-predicted differentiation trajectories for each stromal lineage. Colors indicate the ordering of cells from least (black) to most (red) differentiated, with gray cells excluded from the ordering.(D) Heatmaps of dynamically varying genes for each lineage, with cells ordered from least to most differentiated and genes ordered by the clustering of expression patterns. Gene expression was smoothed using a Gaussian kernel along the cell ordering axis.(E) Clustering of gene expression traces for significantly variable genes along the PBA-predicted osteoblast differentiation ordering. Z score normalized, Gaussian-smoothed expression traces were clustered using k-means clustering. Individual gene traces are shown in black, and the average cluster trace is shown in red.

Figure 4.. Validation of Predicted Gene Expression States and their Differentiation Potential

(A) Flow cytometry of tdTOM⁺ cells for each marked Cre reporter, showing the expression of cluster-specific surface markers.(B) Sorted tdTOM⁺ cells were used for qPCR of lineage-specific genes. Error bars represent the SDs of three replicates.(C) Cell surface markers were used to sort cells for qPCR of known lineage-specific genes. Left: the relative gene expression based on the qPCR measurements of sorted cells. Error bars represent the SDs of four replicates. Right: the relative average gene expression in the scRNA-seq clusters. Error bars represent 95% confidence intervals based on 1,000 bootstrap iterations.(D) Cells were sorted based on specific surface markers, and their differentiation potential was determined. Top: representative images for the differentiation of 100 sorted Pre-OP cells. Bottom: the differentiation potential (percentage of wells with colonies) of each FACS population. Replicates are shown as individual points, and error bars represent the SDs of 4 independent trials, with 12 wells per trial.

Figure 5.. Validation of Predicted Lineage-Specific Transcription Factors in Cultured Stromal Cells

(A) SPRING plot quantifying the mapping of cultured stromal cells to their most similar in vivo counterparts.(B) SPRING plots showing the expression of key stromal cell transcription factors.(C) Representative images of adipocytes, osteoblasts, and chondrocytes with varying degrees of differentiation following shRNA knockdown in passage 1 cultured MSCs. Scale bars represent 50 μm.(D) Transcription factor knockdown using shRNA impeded the ability of stromal cells to differentiate into specific lineages. Differentiation potential was measured as an area of staining in wells, relative to controls. Error bars represent the SDs of four independent experiments.