Identification of phenotypically, functionally, and anatomically distinct stromal niche populations in human bone marrow based on single-cell RNA sequencing

Hematopoiesis is regulated by the bone marrow (BM) stroma. However, cellular identities and functions of the different BM stromal elements in humans remain poorly defined. Based on single-cell RNA sequencing (scRNAseq), we systematically characterized the human non-hematopoietic BM stromal compartment and we investigated stromal cell regulation principles based on the RNA velocity analysis using scVelo and studied the interactions between the human BM stromal cells and hematopoietic cells based on ligand-receptor (LR) expression using CellPhoneDB. scRNAseq led to the identification of six transcriptionally and functionally distinct stromal cell populations. Stromal cell differentiation hierarchy was recapitulated based on RNA velocity analysis and in vitro proliferation capacities and differentiation potentials. Potential key factors that might govern the transition from stem and progenitor cells to fate-committed cells were identified. In situ localization analysis demonstrated that different stromal cells were localized in different niches in the bone marrow. In silico cell-cell communication analysis further predicted that different stromal cell types might regulate hematopoiesis through distinct mechanisms. These findings provide the basis for a comprehensive understanding of the cellular complexity of the human BM microenvironment and the intricate stroma-hematopoiesis crosstalk mechanisms, thus refining our current view on human hematopoietic niche organization.

identities and functions of the different BM stromal elements in humans remain poorly defined. Based on single-cell RNA sequencing (scRNAseq), we systematically characterized the human nonhematopoietic BM stromal compartment and we investigated stromal cell regulation principles based on the RNA velocity analysis using scVelo and studied the interactions between the human BM stromal cells and hematopoietic cells based on ligand-receptor (LR) expression using Cell-PhoneDB. scRNAseq led to the identification of six transcriptionally and functionally distinct stromal cell populations. Stromal cell differentiation hierarchy was recapitulated based on RNA velocity analysis and in vitro proliferation capacities and differentiation potentials. Potential key factors that might govern the transition from stem and progenitor cells to fate-committed cells were identified. In situ localization analysis demonstrated that different stromal cells were localized in different niches in the bone marrow. In silico cell-cell communication analysis further predicted that different stromal cell types might regulate hematopoiesis through distinct mechanisms. These findings provide the basis for a comprehensive understanding of the cellular complexity of the human BM microenvironment and the intricate stroma-hematopoiesis crosstalk mechanisms, thus refining our current view on human hematopoietic niche organization.

Editor's evaluation
The manuscript by Li and coworkers is a landmark characterization of sorted human nonhematopoietic bone marrow cells by scRNA-seq, which predicts their potential lineage relationships and possible interactions with mature and immature hematopoietic cells. Transcriptionally-different stromal cell subsets are identified convincingly, and their lineage relationships, cell-cell interactions and possible specialized functions are predicted from in-silico studies, paving the way for future necessary functional validation studies. This resource significantly adds to the current understanding of human non-hematopoietic bone marrow stromal cells and their hematopoietic regulatory functions.

Introduction
Bone marrow (BM) is the principal site of hematopoiesis in adult humans. In the BM, hematopoietic stem cells (HSCs) and their progenies are contained in specialized microenvironments which regulate HSC maintenance and differentiation.
In contrast, important HME elements in human BM remain poorly defined, which is due to the fact that the identity and function of the stromal stem/progenitor cells have been difficult to investigate (Li et al., 2016). Nevertheless, a number of recent studies have provided the first important insights into the complexity of the human BM and the potential diverse functional roles of BM stromal cells (Chan et al., 2018;de Jong et al., 2021;Triana et al., 2021;Wang et al., 2021).
Based on single-cell RNA sequencing (scRNAseq) technology, we herein investigated the human non-hematopoietic BM cell compartment aiming to resolve the composition of the human HME at the highest possible resolution, to identify potential novel marrow stromal subsets and cellular hierarchies as well as to establish functional relationships between stromal and hematopoietic elements.

Materials and methods
Human bone marrow mononuclear cell isolation Human bone marrow (BM) cells were collected at the Hematology Department, Skåne University Hospital Lund, Sweden, from consenting healthy donors by aspiration from the iliac crest as described previously . Bone marrow mononuclear cells (BM-MNC) were isolated by density gradient centrifugation (LSM 1077 Lymphocyte, PAA, Pasching, Austria). The plasma layer supernatant containing adipocytes was discarded. The use of human samples was approved by the Regional Ethics Review Board in Lund, Sweden.

Flow cytometry and fluorescence activated cell sorting (FACS)
For sorting of BM-MNCs for scRNAseq analysis, freshly isolated BM-MNCs were incubated in blocking buffer [DPBS w/o Ca2+,Mg2+,3.3 mg/ml human normal immunoglobulin (Gammanorm, Octapharm, Stockholm, Sweden), 1% FBS (Invitrogen)], followed by staining with monoclonal antibodies against CD271, CD235a, and CD45. Sorting gates were set according to the corresponding fluorescenceminus-one (FMO) controls and cells were sorted on a FACS Aria II or Aria III (BD Bioscience, Erembodegem, Belgium). Dead cells were excluded by 7-Amino-actinomycin (7-AAD, Sigma) staining and doublets were excluded by gating on FSC-H versus FSC-W and SSC-H versus SSC-W. Cells were directly sorted into PBS with 0.04% BSA for scRNAseq analysis. Sorting of BM-MNCs for CFU-F or stromal culture was performed using DAPI for dead cell exclusion and cells were directly sorted into stromal cell culture medium [StemMACS MSC Expansion Media, human (Miltenyi Biotec, Bergisch Gladbach, Germany)].

CFU-F (colony-forming unit, fibroblast) assay and generation of cultured stromal cells (cSCs)
Unsorted and FACS-sorted primary BM-MNC were cultured at plating densities of 5×10 4 cells/cm 2 cells and 10-50 cells/cm 2 , respectively. Colonies (≥40 cells) were counted after 14 days (1% Crystal Violet, Sigma). Assays were set up in duplicates or triplicates. For single-cell CFU-F assays, cells were sorted into 96-well plates, cultured in stromal cell culture medium, and colonies were counted after 3 weeks. Thereafter, cells were harvested and culture-expanded for use in subsequent experiments. Medium was changed weekly and cultured stromal cells were passaged at 80% confluency after trypsinization (0.05% trypsin/EDTA, Invitrogen, Carlsbad, USA).
After coverslip mounting, all slides were scanned with a Zeiss 780 Confocal Laser Scanning Microscope. 3D orthographic cross-section views and intensity profiles were generated with the ZEN Microscopy Software (Zeiss).

Immunofluorescence image processing
All image data were analyzed and processed using arivis Vision4D. The sequential regions were aligned and overlayed based on nuclei positions. Sequentially used fluorophores were assigned pseudo colors. The obtained five color stacks were cropped to an overlapping region and a maximum intensity projection of the three-layer z-stack was generated.
Single-cell RNA sequencing and data analysis Single-cell RNA sequencing Single-cell RNA sequencing (scRNAseq) was performed by the Single-Cell Genomics Platform at the Center for Translational Genomics at Lund University. FACS-sorted BM-MNC populations analyzed were CD45 low/-CD235acells from four donors (two males, two females, ages 19, 22, 52, and 53 years) and CD45 low/-CD235a -CD271 + from five donors (three males, two females, ages 21, 25, 32, 58, and 61 years) (Supplementary file 1). Donors younger than 35 years or older than 50 were considered as young donors and old donors, respectively. FACS-sorted BM-MNC cells were encapsulated into emulsion droplets as single cells using the Chromium Controller (10 X Genomics). scRNAseq libraries were constructed using the Chromium Single Cell Gene Expression 3' v3 Reagent Kit according to the manufacturer's instructions. Reverse transcription and library preparation was performed on a C1000 Touch Thermal Cycler (Bio-Rad). Amplified cDNA and final libraries concentrations were measured with a Qubit 4 Fluorometer (Invitrogen) using the dsDNA High Sensitivity Assay (Invitrogen). cDNA and library traces were evaluated on a TapeStation (Agilent) using High Sensitivity D5000 and D1000 Screen Tapes (Agilent). Individual libraries were diluted to 1.5 nM and pooled for sequencing. Pools were sequenced on a NovaSeq 6000 (Illumina) aiming for a sequencing depth of 50,000 reads per cell.
Pre-processing of scRNA-seq data, cell filtering, and normalization scRNAseq data of both CD45 low/-CD235aand CD45 low/-CD235a -CD271 + populations were processed and combined for further data analysis. Sequencing data were demultiplexed and UMI-collapsed using the CellRanger toolkit (version 3.1, 10 X Genomics). Common QC metrics including cell filtration, downsampling, gene filtration, mitochondria and ribosomal reads exclusion were used to identify low-quality cells. Reads were mapped against the human GRCh38 genome paired with the gencode. v31 gene level data using the kallisto/bustools package (version 0.46.0 resp. 0.39.3). Resulting count metrics from the kallisto/bustools were analyzed using Python scanpy/velocyto. Mitochondrial and ribosomal reads were excluded and cells with less than 1000 (diploid) detected UMIs were removed from the analysis. Expression was normalized using the scanpy. pp. downsample_ counts function to 1000 UMIs. Finally, genes not expressed in at least 10 cells were also removed from the further analysis.

Dimensionality reduction
Dimensionality reduction was performed using gene expression data for a subset of 3000 variable genes with the scanpy. pp. highly_ variable_ genes function. UMAP was generated with the default setting of Scanpy version 1.6.0 using the top 50 principal components (PC).

Clustering and visualization
Graph-based clustering of the PCA-reduced data was performed using the Louvain method implemented in scanpy. The resolution parameter was set to 11, which resulted in 103 preliminary clusters that were merged based on the correlation of mean expression for each cluster. Clusters with correlations greater than 0.95 were merged, resulting in the final clusters displayed in Figure 1B. Cell type annotation was performed based on the published marker genes expression (Supplementary file 2; Baryawno et al., 2019).

Differential expression of gene signature
For each cluster, we used the scanpy. tl. rank_ genes_ groups function with default values to identify genes that had a significant fold change of higher than 2 in comparison to the other clusters.

RNA velocity analysis
Bone marrow stromal cell dynamics were analyzed using the scVelo package (version 0.2.2) in Scanpy 1.6.0 (Svensson and Pachter, 2018;Bergen et al., 2020). The data were processed using default parameters following preprocessing as described in scVelo package. Analysis of cellular trajectory inference by RNA velocity was performed using dynamical modeling, which is a generalization of the original RNA velocity method and which allows for characterization of multiple transcriptional states. Latent time analysis, pseudotime analysis and velocity confidence were computed using default parameters. Extrapolated states were then projected on the UMAP embedding produced during the initial analysis step. The clustering and visualization were repeated using the same parameters as above for the selected stromal clusters forming the continuum.
For velocity analysis, the samples were pre-processed using functions for detection of the minimum number of counts, filtering and normalization using scv. pp. filter_ and_ normalise and followed by scv. pp. moments function. The gene-specific velocities were then calculated using scv. tl. velocity with mode set to stochastic, and visualized using scv. pl. velocity_ embedding function. In addition, we used scv. tl. latent_ time function to infer a shared latent time from splicing dynamics and plotted the genes along a time axis sorted by expression along dynamics using scv. pl. heatmap function. Pseudotime trajectory and velocity confidence were computed using velocity_pseudotime, and scv. tl. velocity_ confidence functions.

Cell-cell communication analysis
Cell-cell interaction analysis was conducted with CellPhoneDB2 (Efremova et al., 2020), which is a Python-based analytical tool for calculating the interaction between ligands and receptors between different cell populations. The mean and p value of each ligand-receptor interaction between different clusters were calculated by CellPhoneDB2. To illustrate the strength of specific pathways between different clusters, several common pathways with high mean values were selected for visualization. Dot plots were generated with R studio (Version 1.3.1093). The list of ligand-receptor interaction pairs is shown in Supplementary file 6.

Code availability
The code generated to analyze the single-cell datasets in this study is available through GitHub at https://github.com/Hongzhe2022/MSC_BM_scripts, copy archived at Stefan, 2022.

Single-cell RNAseq identifies distinct cell populations in the human bone marrow microenvironment
We explored the cellular composition of the normal BM stroma by scRNAseq profiling of the nonhematopoietic cell-containing CD45 low/-CD235apopulation from healthy donors ( Figure 1A). The well-known low frequency of the BM stromal stem/progenitor cells (around 0.001-0.01% of BM cells) (Pittenger et al., 1999) imposed challenges to detect CXCL12-expressing stromal cells in the CD45 low/-CD235apopulation ( Figure 1-figure supplement 1A). In order to also capture and be able to analyze extremely rare stromal cell subsets, we therefore also sorted CD45 low/-CD235a -CD271 + cells,  Tormin et al., 2011;Li et al., 2014). As CD45 low/-CD235a -CD271 + cells represent a subset of CD45 low/-CD235acells, we combined the two datasets which allowed for a detailed analysis of the structural and developmental organization of the BM stroma at the highest possible resolution.
In total, our dataset was comprised of 25,067 cells (median of 856 genes and 2937 UMIs per cell, Of note, following the exclusion of hematopoietic cells and the other non-stromal cell types, we identified nine stromal clusters (clusters 3,5,6,8,16,23,29,37,38,Figure 1B and C), which only became detectable after CD45 low/-CD235a -CD271 + cells were included in the analysis. These findings thus validated our CD271 + cell enrichment strategy and proved that a high-resolution transcriptomic analysis of the thus far poorly defined human BM stromal compartment would not have been achieved without CD45 low/-CD235a -CD271 + enrichment.

Cellular heterogeneity of stromal progenitors
The nine stromal progenitor populations shared a unique gene expression profile, including both well-established stromal markers as well as novel markers (Supplementary file 2; Supplementary file 4; Supplementary file 7; Li et al., 2014). These markers clearly distinguished stromal cells from the other non-hematopoietic cell types as demonstrated by the differential expression of several stromal genes ( Figure 1D, and Figure 1-figure supplement 1H).
CD235a -CD271 + cells (77.24%)] from a total of nine healthy donors. Color legend indicates cluster numbers and annotations. Basal, basal cell-like cluster; B prog., B cell progenitors; SC, stromal cells; PC, plasma cells; HSPC, hematopoietic stem and progenitor cells; Mk, megakaryocytes; DC, dendritic cells; GC, granulocytes; Ery., erythroid cells; NK, natural killer cells; EC, endothelial cells; Neuronal, neuronal cell-containing cluster. (C) UMAP display of non-hematopoietic clusters (n=9686 cells). The circled area in the overview (top left corner) indicates the non-hematopoietic clusters selected from (B). (D) Heatmap of representative differentially expressed genes for each of the non-hematopoietic clusters in (C). For stromal marker identification, the top 100 significant differentially expressed genes from each stromal cluster were selected to identify overlapping genes. Three genes shared by all nine clusters were identified in this comparison (NNMT, IFITM3, and DCN). Nine more genes were identified when OC clusters 23 and 29 (identified as OCs in Figure 2A) were removed from the comparison, as OCs showed considerably different gene expression profiles as compared with other stromal clusters. Cluster numbers and corresponding cell types are indicated under the heatmap. The scale bar indicates gene expression levels. EC, endothelial cells. A blow-up heatmap for clusters 38, 29, 23, 6, and 37 is shown under the main heatmap for better visualization.
The online version of this article includes the following figure supplement(s) for figure 1:   Ambrosi et al., 2017). Osteogenic differentiation markers were also expressed but at considerably lower levels compared to adipogenic genes ( Figure 2D). Thus, cluster 5 was termed as highly adipocytic gene-expressing progenitors (HAGEPs). Cluster 16 on the other hand, showed similar expression levels of adipogenic and osteogenic markers and was therefore annotated as balanced progenitors (Figure 2A-D). We observed a gradual increase of more mature osteogenic markers such as RUNX1, CDH11, EBF1, and EBF3 from cluster 5 to cluster 16 and to cluster 38 ( Figure 2D). Accordingly, cluster 38 was denoted as pre-osteoblasts. Clusters 29 and 23 were characterized by prominent expression of osteochondrogenic markers including BGLAP, CHAD, SPP1(OPN), RUNX2, CDH11, CDH2, IBSP (BSP), and SPARC ( Figure 2B and E, Figure 2-figure supplement 1A). Hence, these two clusters were assigned as osteochondrogenic progenitors (OCs). While both adipogenic and osteogenic marker expression was detected in clusters 6, 37 and 8, they also expressed several hematopoiesis-related markers such as SRGN, CD52, CD37, CD48, and PTPRC ( Figure 2B and Figure 2-figure supplement 1A and D). Furthermore, S100 genes, which encode a group of calcium-binding cytosolic proteins, were expressed almost exclusively in cluster 8 ( Figure 2B and F), and the previously reported murine fibroblast markers, S100A4 and S100 A6 (Słomnicki andLeśniak, 2010, Baryawno et al., 2019) were highly expressed in all three clusters. Based on this expression pattern, clusters 6, 37, and 8 were subsequently annotated as pre-fibroblasts.

Inferred trajectories of stromal cell gene expression reconstruct the temporal sequence of transcriptomic events during stromal cell differentiation
To resolve the differentiation dynamics of stromal subsets and infer the directionality of individual cells during differentiation, RNA velocity analysis was performed using the scVelo toolkit based on the changes in the spliced/unspliced ratio of mRNA counts (Bergen et al., 2020). Figure 3A and A' illustrate the transcriptional dynamics as indicated by streamlines with arrows pointing from the MSSC cluster to different progenitor clusters. RNA velocity analysis predicted that the MSSC cluster has two main developmental directions ( Figure 3A') with the first originating from MSSCs (cluster 3) and moving towards HAGEPs (cluster 5), followed by OCs (cluster 29) and terminating at the preosteoblast cluster (cluster 38). The second main differentiation path, which also originated from MSSCs (cluster 3), was directed towards balanced progenitors and ended at either the pre-osteoblasts or the pre-fibroblasts stage. A small fraction within HAGEPs also showed development potential towards balanced progenitors and pre-osteoblasts ( Figure 3A and A'). Latent time analysis, which represents the cell's internal clock and approximates the real time experienced by cells as they differentiate, also identified MSSCs as the root cell population and placed pre-osteoblasts close to the terminal state of differentiation ( Figure 3B).
To provide a more precise approximation of the main differentiation stages, we performed latent time and pseudotime analysis on only those clusters that formed a continuum of cellular differentiation states in UMAP, that is MSSCs, HAGEPs, balanced progenitors, OCs (cluster 29), pre-osteoblasts and one of the pre-fibroblast clusters (cluster 6) ( Figure 3C). Thereby, disturbances of the data analysis by larger gaps in cellular state transitions between the clusters within the continuum and the outlying OCs (cluster 23) and pre-fibroblasts (clusters 8 and 37) were avoided ( Figure 3A). The results of these analyses using two different methods which are based on different mathematical algorithms further corroborated our initial conclusion that the MSSC cluster was at the apex of a differentiation hierarchy with downstream differentiation paths into HAGEPs, balanced progenitors, OCs and pre-osteoblasts plots of adipogenic-(C), osteogenic-(D), osteochondrogenic-(E), and fibroblastic (F) markers in different stromal clusters. Dashed lines separate MSSCs, HAGEPs, balanced progenitors, and pre-osteoblast (left; cluster 3, 5, 16, and 38), OCs (middle; cluster 29, 23), and pre-fibroblasts (right; cluster 6, 37, 8). The boxes in the violin plots indicate the lower, median and upper quartiles.
The online version of this article includes the following figure supplement(s) for figure 2:   Figure 3-figure supplement 1A-B). Furthermore, generally high velocity confidences indicated that RNA velocity calculations were reliable (Figure 3-figure supplement 1C).
In addition to the identification of stromal differentiation directions and hierarchies, scVelo allowed us to detect potential key driver genes that govern cellular state transitions. Using likelihood-based computation, we instructed the program to identify the top 300 genes whose transcriptional dynamics correlated with cellular state transitions along the inferred latent time in the stromal continuum ( Figure 3C, Supplementary file 5). A group of genes, including CSTA, B2M, CXCL12, and COL1A2 was found to drive the gradual and subtle cellular state transitions within the MSSC cluster ( Figure 3C). Known adipogenic genes such as CEBPD (Cao et al., 1991) together with LEPR, CTSC, LMNA, KLF6, and EGR1 were identified to play important roles in the transition from MSSCs to HAGEPs ( Figure 3C). As exemplified by CEBPD, EGR1, and LMNA ( Figure 3-figure supplement 1D-E), gene transcriptional activity and expression level gradually increased from the MSSCs to HAGEPs. Furthermore, both novel drivers (NFKBIA, RBM3, TUBA1A, and DDX24), as well as reported genes such as ANXA1 (Headland et al., 2015), were identified to serve as critical genes in the transition from HAGEPs to OCs ( Figure 3C and Figure 3-figure supplement 1F). VCAN expression has been demonstrated to correlate with bone formation in rats (Nakamura et al., 2005). Accordingly, VCAN was identified as one of the crucial genes for the transition from MSSCs to balanced progenitors and even to the pre-osteoblasts. Other potential novel driver genes along this path included CHL1 and TXNIP (Figure 3-figure supplement 1G) and NEAT1 and EIF5 (not shown). These genes demonstrated elevated expression levels in balanced progenitors and pre-osteoblasts.
In summary, RNA velocity analysis of the transcriptional dynamics allowed us to postulate stromal cell differentiation pathways and identify regulatory factors, including putative driver genes, that are potential key candidates to govern stromal cell fate commitment.
Surface molecule profiling identifies marker combinations for the prospective isolation of functionally distinct stromal progenitor populations Based on the expression patterns of stromal markers, osteochondrogenic genes, and fibroblastic markers, the nine stromal subsets could be further assigned to three main groups (Figure 2 and Figure 4-figure supplement 1A). Group A cells (clusters 3, 5, 16, 38), which included MSSCs, HAGEPs, balanced progenitors, and pre-osteoblasts, showed high expression of stromal markers, adipogenic and osteogenic differentiation genes. Group B cells (clusters 29 and 23), annotated as OCs (Figure 2), demonstrated exclusive and specific expressions of osteochondro-lineage markers. Group C cells (clusters 6, 37, and 8), assigned as pre-fibroblasts (Figure 2), exhibited characteristic expression of a group of S100 genes, such as S100A4, S100A6, S100A8, and S100A9. These clear differences in gene expression between the groups pointed to potential functional differences. We therefore went on to identify suitable cell surface markers and marker combinations, respectively, that would allow for the prospective isolation and functional investigation of the three groups of cells.
As shown in Figure 4A, surface marker gene expression differed between groups A, B, and C. Group B cells showed exclusive expression of NCAM1 and CD9, which were considerably lower expressed by group A and C cells ( Figure 4A). Group C cells expressed several surface molecules which are conventionally considered as hematopoietic markers (CD52, CD37, PTPRC, and CD48). Thus, the lack of expression of NCAM1, CD9 and group C-specific antigens could be used to identify stromal cells colored by inferred latent time. Inferred latent time is represented by a color scale from 0 (the earliest latent time) to 1 (the latest latent time). (C) Left: UMAP display of stromal cell clusters that form a continuum of different cellular states (n=6376 cells). Colors correspond to different stromal clusters as in Figure 2A. Right: Heatmap constructed by the top 300 likelihood-ranked genes demonstrates gene expression dynamics along latent time. Colors on top of the heatmap correspond with cluster colors in UMAP (left). Key genes are highlighted by different colors on the right. Gene name colors correspond to different developmental stage transitions. Red: genes responsible for MSSC intra-cluster transition; orange: genes responsible for the transition to HAGEPs; light green: genes responsible for the transition to OCs; dark green: genes responsible for the transition to balanced progenitors and pre-osteoblasts.
Taken together, the data-driven FACS gating strategy allowed us to isolate distinct stromal subsets for functional investigation. Both MSSCs and HAGEPs demonstrated considerably higher in vitro colony-forming capacities in comparison with balanced progenitors, OCs, pre-osteoblasts, and prefibroblasts. While MSSCs and HAGEPs exhibited full multi-differentiation potentials, pre-fibroblasts had only limited differentiation potentials. OCs showed enhanced osteoblastic and chondrogenic differentiation capacities and compromised adipogenic potential, which is consistent with their differentiation marker expression profiles (Figure 2).

In-situ localization of stromal cell subpopulations
Based on the markers identified above, we used CD271 (NGFR), CD81, NCAM1, and CD45 to identify stromal cells (CD271 + ) and to investigate the in-situ localization of MSSCs (CD271 + CD81 ++ ) as well as OCs (CD271 + NCAM1 + ) using bone marrow biopsies from hematologically normal donors. As shown in Figure 4E and Figure 4-figure supplement 1D-E, CD271-expressing cells were found in perivascular regions, endosteal regions, and throughout the stroma. We observed that perivascular cells surrounding the capillary endothelium and larger vessels primarily expressed CD271 alone. The bone-lining cells proximal to the surface of trabecular bone showed expression of both CD271 and NCAM1, confirming the osteochondrogenic nature of OCs, which correlated very well with the predicted endosteal localizations of murine osteoblasts and chondrocytes (Baccin et al., 2020). In contrast to the CD271/NCAM1 double-positive cells, CD271/CD81 double-positive cells were localized in bone marrow stromal regions, including both peri-adipocytic and perivascular regions with a considerable fraction of CD271/CD81 cells being located in close proximity or encircling adipocytes. Co-staining with CD45 confirmed that perivascular, periadipocytic, and bone-lining cells were CD45 negative. The results from the stromal cell in situ localization analysis by scanning microscopy were confirmed by detailed confocal microscopy analysis as illustrated in Figure 4-figure supplements 2-16. These results thus demonstrated the distinct anatomical localizations of MSSCs and OCs, which confirms and extends our previous findings on differently localized putative hematopoietic niche cells (Tormin et al., 2011) and point to possible divergent physiological functions of the different stromal subsets.

In silico prediction of cell-cell interactions in the human bone marrow microenvironment
It has been proposed that different niches exist for different types of hematopoietic stem and progenitor cells. We therefore went on to study the interactions between the different human BM stromal cell populations identified herein and hematopoietic cells based on ligand-receptor (LR) expression using CellPhoneDB (Efremova et al., 2020) complemented by analysis of additional published LR pairs for erythroid cells (Kleven et al., 2018).
As shown in Figure 5A and Figure 5-figure supplement 1A, we identified a wide range of interactions between stromal cells and hematopoietic cells. Generally, cells belonging to the three stromal groups (A-C) were predicted to have multiple interactions with a large number of hematopoietic cell types, including HSPCs and more mature hematopoietic cells ( Figure 5A and  . Interestingly, only a few LR pairs were identified between stromal cells and plasma cells ( Figure 5A).
To study the functional relationships between stromal and hematopoietic cells in more detail and to identify subset-specific interactions, we plotted selected ligand-receptor pairs that are known to be relevant for hematopoietic-stromal crosstalk ( Figure 5B). Selected LR pairs, including essential cytokines and established hematopoiesis-supporting molecules were analyzed for effects from stromal cells on hematopoiesis and vice versa, including stromal subset-specific interactions, intra-, and interstromal cell interactions, and other potentially important interactions ( Figure 5B, Figure 5-figure supplement 1B, Figure 5-figure supplement 3B).

Stromal cell regulation of hematopoiesis
Regarding the regulatory and supportive effects of stromal cells on hematopoiesis, we found that group A and C stromal subsets had the potential to interact with essentially all hematopoietic cell types through similar pathways. CXCL12 was highly expressed by group A and C stromal cells, especially the MSSC population ( Figure 6A), and CXCL12-CXCR4 crosstalk was detected in almost all hematopoietic clusters except for erythroid progenitors ( Figure 5B, Figure 5-figure supplement 1B, and Figure 6A). Furthermore, large dot sizes and high statistical power indicated that CXCL12mediated interactions represent one of the major interaction mechanisms between these stromal cells and hematopoietic cells ( Figure 5B, Figure 5-figure supplement 1B). In contrast, CXCL12 interactions with CXCR3 and DPP4 were less significant ( Figures 5B and 6A).
KITLG (SCF)-involving stroma-hematopoiesis interactions were also identified. However, these interactions were not as strong as those involving CXCL12 ( Figure 5-figure supplement 2A). Interestingly, MDK, a neurite growth factor expressed by stromal cells, showed interactions with the HSPCs, lymphoid lineage progenitors (B, T, and NK), and monocytes through LRP1 or SORL1 ( Figure 5-figure supplement 2B), which has thus far only been described in mouse fetal liver (Gao et al., 2022). Furthermore, other cytokines known to be important for the maintenance of CD34 + HSPC cells include FLT3-ligand (FLT3LG) and thrombopoietin (THPO). However, whereas expression of FLT3 was detected in several hematopoietic cell types including HSPCs, dendritic cells and monocytes ( Figure 5-figure supplement 2C), THPO and MPL expression were barely detected ( Figure 5figure supplement 2D). JAG1 expression was highly enriched in stromal cells while NOTCH1, 2 and 4 expressions were detected in hematopoietic cells, suggesting that activation of Notch signaling is involved in stromal-hematopoietic interplay ( Figure 5-figure supplement 2E). All three stromal groups were furthermore predicted to interact with T cells through the IL7-IL7R RL combination, indicating a T cell supporting mechanism by stromal cells (Figures 5B and 6B, Figure 5-figure  supplement 1B).
In addition to soluble factor-mediated interactions, in silico analysis also identified direct cell-to-cell communication between group A and C stromal cells and hematopoietic cells mediated by VCAM1 ( Figure 6C), FN1, PLAUR and TNC, some of which have also been identified in a previously reported study on murine cells (Mende et al., 2019). These results indicated that group A and C cells have of selected ligand-receptor pairs between MSSCs and OCs, respectively, and different hematopoietic cell types from this dataset (B) and more defined hematopoietic progenitors from another study (C). The y-axis label indicates the pair of interacting cell clusters ('cluster X_cluster Y' indicates cluster X interaction with cluster Y by ligand-receptor expression). The x-axis indicates the interacting receptor/ligand (R/L) or ligand/receptor pairs, respectively ('molecule L_molecule R', molecule L interacts with molecule R). The means of the average expression level of the interacting molecules are indicated by circle size (adjacent lower scale bar). P values are indicated by circle color and correspond to the upper scale bar. Cluster numbers of the hematopoietic cell types are listed in the table. The dashed line separates MSSC from OC interactions. DC-I, dendritic cell progenitor cluster I; GMP, granulocyte-macrophage progenitor cluster I; HSC-I, hematopoietic stem cell cluster I; HSC-II, hematopoietic stem cell cluster II; Ly-I, lymphoid progenitor cluster I; MEP-I, megakaryocyte-erythroid progenitor cluster I; MPP-I, multipotent progenitor cluster I.
The online version of this article includes the following figure supplement(s) for figure 5:     the potential to support hematopoietic cells via direct cell-to-cell contact as well as hematopoiesissupporting cytokines.
In contrast to the strong contribution of CXCL12-mediated interactions between group A and C stromal cells and hematopoietic cells, SPP1-mediated interactions were predominant in the crosstalk between group B stromal cells (OCs) and hematopoietic cells (except for megakaryocytes), indicating a stromal group-specific interaction ( Figure 5B and Figure 5-figure supplement 1B). Interacting partners of SPP1-expressing OCs were erythroid cells, lymphoid progenitors (B, T, NK progenitors) and CD34-enriched HSPC cluster which expressed SPP1 binding partners such as CD44 and PTGER4 ( Figures 5B and 6D and Figure 5-figure supplement 1B). SPP1-involving integrin alpha 4 beta 1 (a4b1) and alpha 9 beta 1 (a9b1) complexes could also play important roles in the interaction of OCs with a number of hematopoietic cells. Consistent with this, FACS analysis demonstrated that SPP1 expression was higher in OCs in comparison with MSSCs ( Figure 5-figure supplement 2F). Moreover, in situ staining demonstrated that SPP1-expressing cells were exclusively colocalized endosteally with NCAM1 + OCs ( Figure 6E), pointing to a localization-specific regulation. Importantly, the CXCL12-and SPP1-mediated stromal cell type-specific interactions detected in our dataset were furthermore confirmed when performing R/L analysis of our stroma cell data with well-defined human BM hematopoietic stem and progenitor clusters obtained from a recently published study (Sommarin et al., 2021; Figure 5C).
Taken together, these data demonstrated that the most predominant interactions were CXCL12and SPP1-mediated interactions for group A/C stromal cells and OCs, respectively, indicating that regulation of hematopoiesis is stromal group and location-specific.

Erythroid cell-specific interactions
Unlike the majority of the hematopoietic cells in our dataset, the interplay between erythroid cells and group A and C stromal cells was not dependent on CXCL12-involving interactions as indicated by lacking expression of CXCL12 receptors in erythroid clusters ( Figure 6A). CellPhoneDB analysis demonstrated that while group A and C stromal cells had the potential to interact with the erythroid cells through VCAM1-, FN1-, TNC-and PLAUR-involving complexes, the interactions between OCs (group B) and erythroid cells were mainly mediated by . In addition to that, we found that TF (transferrin) was exclusively expressed by all stromal clusters and its corresponding receptors, TFRC and TFR2 were highly expressed by erythroid progenitors ( Figure 6F).

Hematopoietic regulation of stromal cells
Aside from the supportive and regulatory role of stromal cells for hematopoiesis, it is well described that regulatory signals originating from hematopoietic cell types can affect stromal cells (Baksh et al., 2005). Accordingly, we found that the crucial stromal cell growth factors PDGFs (PDGFA, PDGFB, and PDGFC) were highly expressed by megakaryocytes while their corresponding receptors (PDGFRA and PDGFRB) were predominantly expressed by stromal cells (Figure 6G), suggesting interactions through these axes ( Figure 5B, Figure 5-figure supplement 1B). Our analysis also predicted that hematopoietic cells control stromal cells via the TGFB1-EGFR and TGFB1-TGFBR1/2/3 axes ( Figure 5B, Together, these data indicate that the regulatory effects between stromal cells and hematopoietic cells are bidirectional, which is consistent with the increased fibroblastic colony size in the presence of CD45 + hematopoietic cells in CFU-F assays ( Figure 5-figure supplement 2I). (pink), and NCAM1 (cyan) and scanned with an OlympusVS120 slide scanner. Left upper image: all markers are shown; left lower image: just NCAM1 channel is shown. Middle panel: confocal laser scanning analysis of BM biopsies co-stained with CD271 (red), NCAM1 (green), SPP1 (white), and DAPI (blue) presented as 3D orthographic cross-section view. Right panel: Single staining channel data for CD271 (red), NCAM1 (green), SPP1 (white), and corresponding blow-ups for indicated areas. NCAM1, SPP1 and CD271(NGFR). White scale bars represent 50 μm and yellow scale bars represent 25 μm. White dashed lines indicate the trabecular bone surface lining regions.

Figure 6 continued
Inter-and intra-stromal signaling As shown in Figure 5-figure supplement 3A, CFU-F frequencies were significantly reduced when stromal cells were cultured under single-cell conditions in comparison with bulk-cultured cells, suggesting the presence of inter-and intra-stromal signaling among stromal cells. Indeed, multiple stimulatory pathways were operative between different stromal subsets as well as within the same stromal cell population ( Figure 5-figure supplement 3B). Among them, FGF receptor-mediated signaling and various collagen-involving pathways were identified that reflected interactions between different stromal cells ( Figure 5-figure supplement 3C). In addition, CXCL12-, MDK-, GRN-, NRP1-, and BST2-mediated interactions were also found to significantly contribute to the regulatory mechanisms of stromal cells ( Figure 5-figure supplement 3B-D).
In summary, our data showed all three stromal groups demonstrate the potential to interact with a wide range of hematopoietic cells. According to this analysis, group A and C cells were predicted to communicate with hematopoietic cells mainly through CXCL12-, VCAM1-, FN-, and MDK-involving pathways, whereas group B cells have the potential to interact with hematopoiesis mainly via SPP1mediated pathways. Furthermore, regulation was predicted to be bidirectional, that is stromal cells are regulated by hematopoietic cells through PDGF-, OSM-, and TGFB1-mediated signaling. Additionally, multiple intra-and inter-stromal pathways were identified. Our in silico analysis of scRNAseq data of both stromal and hematopoietic clusters thus predicted the presence of a complicated interaction network in human bone marrow microenvironment.

Cell-cell interaction between endothelial cluster and hematopoietic clusters
As endothelial cells are an important niche component, we also compared the interaction pattern between endothelial cells and hematopoietic cells with that of MSSCs. As illustrated in Figure 7A, the endothelial cluster had a broad interaction spectrum with all hematopoietic clusters. In addition to the endothelial-specific SELL-SELPLG and ICAM3-aLb2 complex pathways, it was predicted that the endothelial cluster interacted with hematopoietic clusters mainly through HLA molecules-, TNFRSF family members-and LGALS9-mediated interactions. However, endothelial cluster interaction seemed to be less specific compared with stromal cells. Unlike stromal cells, CXCL12-mediated interactions were not identified as the main signaling pathways in the crosstalk between endothelial cells and hematopoietic cells, indicating that endothelial cells and stromal cells communicate with hematopoietic cells through different mechanisms.

Discussion
The BM microenvironment plays a critical role in regulating hematopoiesis. Recent landmark studies have illustrated the molecular complexity of the murine BM microenvironment (Baryawno et al., 2019;Tikhonova et al., 2019;Baccin et al., 2020). However, the exact definition of the cell populations that form the BM stroma in humans remains elusive. In this study, based on single-cell RNA sequencing technology, we gained important insights into the stromal components of this important organ including their developmental and functional roles in hematopoiesis.
Several recent studies using scRNAseq approaches provided the first important insights into the complexity of the human BM and the potential functional roles of BM stromal cells (de Jong et al., 2021;Triana et al., 2021;Wang et al., 2021). A first overview of the cellular composition of the nonhematopoietic cells in human BM was provided by a recently published single-cell proteo-genomic reference map (Triana et al., 2021) and a single-cell transcriptomic map of non-hematopoietic cells (de Jong et al., 2021). However, due to the low frequencies of the stromal stem/progenitor cells, it has been difficult to investigate the potential heterogeneity of the human BM stromal cell compartment even when employing single cell approaches. Enrichment strategies represent a possible approach for a detailed analysis of the transcriptional diversity of BM stromal cells, as recently reported by Wang et al., 2021. However, in their study, CD271 + cells from femoral shafts from old patients with osteoarthritis and osteoporosis were used, which might not reflect transcriptomic profiles of normal stroma cells from hematologically active bone marrow (Wang et al., 2021). Because of the extremely low frequency of stromal cells in human BM, we chose a sorting strategy that also included CD45 low cells to ensure that no stromal cells were excluded from the analysis. Furthermore, we also added sorted CD271 + cells, which are highly and exclusively enriched for stromal stem/progenitor cells (Tormin et al., 2011). Therefore, despite the expected presence of hematopoietic cells in our dataset, this approach ensured that stromal cell resolution was substantially increased.
Following exclusion of the hematopoietic clusters, we identified primarily nine stromal clusters which could further be assigned to six stromal cell types, that is MSSCs, HAGEPs, balanced progenitors, OCs, pre-osteoblasts, and pre-fibroblasts. Of note, corresponding murine stromal cell types have been reported for most of these cell types (MSSCs, HAGEPs, OCs, pre-osteoblasts, and prefibroblasts), whereas the cluster annotated as balanced progenitors has not been described either in mice or in humans (Baryawno et al., 2019;Tikhonova et al., 2019).
Additionally, we identified three non-hematopoietic cell clusters which contained KRT5-(cluster 0) and neuronal gene-expressing cells (cluster 39), respectively, as well as endothelial cells. KRT5 is a keratin family member and is widely used as a basal cell marker. However, KRT5 expression has also been reported for several immune cell types, such as plasmacytoid dendritic cells and NK cells. As cells in our analysis did not express typical immune cell genes, further characterization of this cluster needs to be performed in future studies. Multiple genes indicating the presence of different cell types were identified in cluster 39. The genes identified with the highest fold changes included several neuronal markers (NEUROD1, CHGB, ELAVL2, ELAVL3, ELAVL4, STMN2, INSM1, ZIC2, NNAT) (Supplementary file 4). Due to the heterogeneity of this cell composition, this cluster was therefore annotated as neuronal cell-containing cluster.
Thus, this analysis allowed us to dissect the human bone marrow non-hematopoietic cell compartment at a thus far unreached molecular level, to resolve the stromal cell heterogeneity and generate a detailed transcriptional fingerprint of distinct stromal populations.
Furthermore, we were able to define phenotypically and functionally distinct stromal subsets, which were predicted to be organized in a hierarchical manner according to RNA velocity analysis. We identified MSSCs as the most primitive population, which is in accordance with published murine and human stromal stem cell studies (Zhou et al., 2014;Baryawno et al., 2019;Wolock et al., 2019). Pre-fibroblasts, pre-osteoblasts, and OCs were placed downstream of MSSCs in the developmental hierarchy, thus confirming the reliability of the velocity analyses. By using a likelihood-based computation, we also detected potential key driver genes that are candidates to govern the transition between different cellular states. Identified genes included both novel and established transcription factors and certainly, these data provide the basis for future important mechanistic studies of the dynamic processes associated with cellular state transition and fate choices.
Self-renewing human skeletal stem cells (hSSC) with osteogenic and chondrogenic potential have been previously identified in human fetal and adult bones based on the phenotype PDPN + CD146 -CD73 + CD164 + (Chan et al., 2018). However, the identity of a tri-lineage multipotent stromal progenitor in adult human bone marrow was not reported. Herein, we were able to provide a phenotypical definition that allowed for the faithful identification and prospective isolation of molecularly defined stromal stem and progenitors from hematopoietically active adult bone marrow. These cells have thus far remained incompletely characterized and our data thus complement and extend reported phenotype definitions of MSCs from bone marrow and skeletal tissues (Jones et al., 2002;Quirici et al., 2002;Gronthos et al., 2007;Sacchetti et al., 2007;Delorme et al., 2008;Li et al., 2016;Chan et al., 2018). The phenotypic MSSCs defined in this study also showed expression of previously reported stromal surface markers such as CD51 (Pinho and Frenette, 2019) and CD146 (Tormin et al., 2011;data not shown). Besides these surface markers, scRNAseq also identified VCAN as a potentially important stromal gene based on its high expression in stromal clusters. However, as an extracellular matrix protein, FACS analysis of cellular VCAN expression can only be achieved based on its intracellular expression after fixation and permeabilization. Additionally, VCAN is also expressed by monocytes (cluster 36). Therefore, VCAN is not an optimal marker to isolate viable stromal cells.
Furthermore, we performed a preliminary analysis which indicated possible age-and gender-related stromal differences as also reported for murine cells (Singh et al., 2016;Maryanovich et al., 2018;Aguilar-Navarro et al., 2020). However, this potentially interesting finding needs to be confirmed by additional experimentation in forthcoming studies.
In mice, elegant studies using genetic approaches involving in-vivo labelling of cell types and ablation of candidate niche cells and niche factors have provided detailed insight into the functional roles of stromal cells in regulating hematopoiesis (Pinho and Frenette, 2019). Comparable studies in human are impossible to perform and analysis of human stromal cells function has been limited to in-vitro models. We therefore investigated cellular interactions based on ligand-receptor interactions and found a broad and complicated crosstalk network, which is likely to reflect -at least in partthe in-situ situation as data were generated with directly isolated cells. We found that stromal cells generally showed hematopoiesis supporting potential, which is consistent with previously published in-vitro coculture data Li et al., 2020) and which we confirmed in co-cultures with BM CD34 + cells and MSSC-or OC-derived stromal feeder cells ( Figure 5-figure supplement 3E). Furthermore, whereas all stromal cell groups showed the potential to interact with almost all hematopoietic cells, cell-cell interaction mechanisms differed considerably between different groups of stromal cell groups, that is between group A/C cells which included the MSSCs and group B cells (OCs). Whereas the former demonstrated the potential to interact with hematopoiesis through cytokines such as CXCL12, IL7, KITLG, MDK, and TF as well as adhesion molecules providing direct cellto-cell contact (e.g. VCAM1), SPP1-mediated interactions were predicted to be the main pathways for the latter. Although we have not formally examined the expression of the receptors for CXCL12 and SPP1 on hematopoietic cells in the current study, expressions of CXCR4 (CXCL12 receptor) and CD44 (SPP1 receptor) on human BM hematopoietic cells has been demonstrated previously (Ishii et al., 1999;AbuSamra et al., 2017).
These results thus indicated that hematopoietic cells were maintained by different stromal populations through diverse but nevertheless stromal cell-specific pathways. Taking furthermore into account that SPP1-expressing OCs were located endosteally and that CD271 single positive stromal cells including MSSCs were localized in the perivascular and stromal regions, respectively, suggested the possibility that different stromal cells provide specialized niches for hematopoietic cells in different locations. Our data thus confirm and extend previous reports describing that endosteal and perivascular niches exist in human BM (Tormin et al., 2011;Pinho et al., 2013), and conceptually support a model where different hematopoietic cell types are differentially regulated by distinct niche milieus which are composed of specific cytokine-producing stromal cells. Finally, in silico ligand-receptor interaction predictions suggested that endothelial cells and stromal cells showed the potential to communicate with hematopoiesis through different mechanisms, which is certainly an interesting finding that needs to be further addressed in future studies.
In this study, we used human BM aspiration samples in which bone-lining or bone-attaching cells are most likely underrepresented. Our experiments could therefore ideally be complemented by studies on fresh bone marrow biopsy samples from healthy individuals, which however were not available to us, and which might be difficult to acquire because of ethical considerations. In our study, tri-lineage differentiation capacities of stromal clusters were evaluated using standard in-vitro assays, which may not necessarily reflect the bona fide stromal cell potentials. Thus, these data need to be complemented in follow-up studies by in vivo differentiation studies, for example using humanized BM ossicles. The potential key genes that govern the stromal cell fate commitment were predicted using likelihoodbased computation and the stroma-hematopoiesis crosstalk mechanisms were inferred based on ligand-receptor interactions. The identified key genes and interactions will need to be validated in future experiments, for example using in vitro engineered human BM models or in vivo humanized BM ossicles (Dupard et al., 2020). While the CellPhoneDB analysis predicted the communication networks between stromal and hematopoietic cells, this approach infers potential interactions using transcriptomics data without considering the spatial proximity of the cells and functional relevance of the interaction. A more comprehensive overview of cellular communication is anticipated when combined with spatial localization analysis and functional assays in future studies. Despite these limitations, we believe that our results provide the basis for a better understanding of the cellular complexity of the human BM microenvironment and put a new light on the current concept that specialized niches exist for distinct types of hematopoietic stem and progenitor cells (Morrison and Scadden, 2014).