Connecting the Dots: Resolving the Bone Marrow Niche Heterogeneity

Single-cell sequencing approaches have transformed our understanding of stem cell systems, including hematopoiesis and its niche within the bone marrow. Recent reports examined the bone marrow microenvironment at single-cell resolution at steady state, following chemotherapy treatment, leukemic onset, and aging. These rapid advancements significantly informed our understanding of bone marrow niche heterogeneity. However, inconsistent representation and nomenclature among the studies hinder a comprehensive interpretation of this body of work. Here, we review recent reports interrogating bone marrow niche architecture and present an integrated overview of the published datasets.


HEMATOPOIESIS
Hematopoiesis is a continuous process of generating blood and immune cells and is one of the best-studied stem cell systems in modern biology (Weissman, 2000). Since the discovery of hematopoietic stem cells (HSCs), the field has been driven by technological advances such as flow cytometry, mass cytometry, high-resolution microscopy, and now single-cell sequencing approaches. In vitro approaches, such as colony-forming assays, defined the intermediate stages between a rare population of multipotent hematopoietic stem cells and the terminally differentiated cell populations. In vivo, early hematopoietic progenitors were further classified as long-term HSCs (LT-HSCs) capable of self-renewal and unlimited differentiation potential and multipotent progenitors (MPPs), characterized by limited self-renewal capacity (Reya et al., 2001). The MPPs have been further separated into myeloid-biased MPP2 and MPP3 as well as lymphoidbiased MPP4 subsets that differentiate to lineage-restricted common myeloid progenitors (CMPs), granulocyte/macrophage progenitors (GMPs), and common lymphoid progenitors (CLPs) (Pietras et al., 2015). Extensive literature continued to define the complexity of hematopoietic differentiation, incorporating further subpopulations and subdivision of downstream progenitors (Cabezas-Wallscheid et al., 2014;Liggett and Sankaran, 2020).
The last decade has witnessed the introduction and wide-scale adoption of high-throughput approaches, including single-cell RNA sequencing (scRNA-seq). Unlike more traditional flow or mass cytometry methods that are limited by predefined markers, scRNA-seq technologies simultaneously capture thousands of genes in each cell, providing an unbiased characterization of their transcriptional diversity even within phenotypically homogenous populations. Notably, a major advantage of single-cell transcriptional profiling is the identification of novel populations and states as well as shifts in their relevant frequencies (Kowalczyk et al., 2015).

HEMATOPOIETIC DIFFERENTIATION AT SINGLE-CELL RESOLUTION
In one of the earlier large-scale single-cell studies, Paul et al. (2015) characterized 2,730 myeloid progenitor cells. Clustering analysis based on an expectation maximization (EM)-like procedure resulted in 19 subpopulations, including groups corresponding to known erythrocyte, monocyte, and neutrophil progenitors, as well as more transcriptionally heterogeneous clusters. The data showed that the traditional GMP and CMP definitions include several subpopulations with distinct lineage commitment transcriptional profiles. Nestorowa et al. (2016) comprehensively profiled individually sorted 1,656 murine hematopoietic stem and progenitor cells (HSPCs), with an average of 6,558 protein-coding genes detected per cell. The study utilized the diffusion maps dimensionality reduction technique (Haghverdi et al., 2015;Angerer et al., 2016) to visualize and interpret the continuous processes of early hematopoietic differentiation. The cells were assigned to 12 commonly sorted HSPC phenotypes based on surface markers. When displayed on the diffusion map, all populations, except CMPs, localized within defined regions. The relative ordering of populations within the diffusion map was in agreement with the established hematopoietic hierarchy and further showed intermediate transition zones containing transcriptionally similar cells that traditional gating assigned to differing phenotypes. Interestingly, the single-cell analysis revealed that the three broad trajectories defined by the diffusion map begin to diverge immediately following the LT-HSC stage.
The findings were extended to the human system. Velten et al. (2017) combined transcriptomic and functional single-cell data to investigate the dynamics of lineage commitment of individual cells. Healthy human HSPCs were defined by the absence of lineage markers and CD34 expression (Lin − CD34 + ). 1,413 HSPCs from the bone marrow of two individuals were profiled by scRNA-seq and 2,038 cells individually cultured ex vivo were used to determine lineage potential. The transcriptomic and functional datasets were integrated based on surface marker expression. Clustering within the Lin − CD34 + CD38 − compartment that includes HSCs and their immediate progeny, such as MPPs, was unstable and the cells formed a continuously connected unit. In contrast, more differentiated Lin − CD34 + CD38 + progenitors formed clusters corresponding to defined progenitor populations. The scRNA-seq data separated the cells into a continuum of low-primed undifferentiated (CLOUD)-HSPCs and discrete populations of restricted progenitors associated with increased CD38 expression. Analysis of the discrete Lin − CD34 + CD38 + sub-populations identified the major branches of hematopoiesis, including B-cell progenitors, megakaryocyte/erythrocyte committed progenitors, neutrophilprimed progenitors, monocyte/dendritic cell progenitors, and eosinophil/basophil/mast cell progenitors. The authors developed a dimensionality reduction technique STEMNET to view the transformation from HSCs to the lineage-restricted progenitors by placing each cell along a trajectory from the least-primed HSCs in the center toward one of the six restricted progenitor populations in the corners of a simplex. The authors proposed a model in which distinct lineages emerge from CLOUD-HSPCs rather than a series of discrete progenitors. Zheng et al. (2018) profiled 21,306 CD34 + progenitor cells from human cord blood across five biological replicates without enrichment or depletion for specific lineages. The analysis revealed previously defined progenitor populations in addition to the intermediate continuous states. The authors additionally integrated the data with the previous bone marrow dataset from Velten et al. to show the high concordance between the two systems. Pellin et al. (2019) transcriptionally mapped the fates of the early human hematopoietic progenitors using 6,011 CD34 + and 15,401 Lin − single cells profiled with the inDrops protocol. The hierarchical continuum of states that branch out toward cells expressing established lineage signatures was visualized using SPRING (Weinreb et al., 2018a). This two-dimensional force-directed graph layout generates a graph of nodes representing cells connected to their nearest neighbors in high-dimensional gene expression space. The organization of the data broadly segregated the cells into known immunophenotypic subpopulations and detected extensive transcriptional heterogeneity among HSPCs. The analysis suggested a structured hierarchy of fate decisions rather than a single step transition from undifferentiated HSPCs to committed states. This organizational structure was confirmed by inferred transcriptional trajectories as well as by population balance analysis (PBA) (Tusi et al., 2018), which in addition to providing a static continuum description of cell states aims to infer dynamic properties such as fate potential. This work confirmed a continuous hierarchical organizational structure from HSPCs to lineage-committed states not fully defined by traditional immunophenotyping. Importantly, the authors compared the human and mouse HSPCs, revealing a comparable branching structure within the two species despite differences in cell surface markers commonly used to isolate the various subpopulations.
Collectively these studies were able to isolate and define highly granular HSPC populations. They highlighted the transitional phases of hematopoietic differentiation, building on the traditional tree-like model of hematopoiesis toward a more fluid path of stem cell differentiation (Laurenti and Göttgens, 2018;Watcham et al., 2019).

BONE MARROW NICHE
Postnatally, HSCs reside in the bone marrow. While substantial progress has been made in understanding the hematopoietic hierarchy, the extrinsic factors guiding HSC self-renewal and differentiation are less defined. The bone marrow is a complex organ composed of numerous cell types that regulate HSC function through physical and biochemical interactions (Pinho and Frenette, 2019). The discovery that SLAM family surface receptors CD150, CD244, and CD48 were differentially expressed among HSPC populations in a way that can track developmental potential made it possible to image and localize those populations within the bone marrow (Kiel et al., 2005). Many HSCs were found associated with sinusoidal endothelial cells (ECs) and the endosteum (Kiel et al., 2005) as well as in contact with Cxcl12abundant reticular (CAR) cells surrounding ECs (Sugiyama et al., 2006). An extensive network of sinusoids in the central marrow allows the hematopoietic cells to transport in and out of the bone marrow. In addition to serving as a barrier, EC-specific deletion of Scf (Kitl), Cxcl12, Jag1, and Dll4 impacted HSC quiescence, differentiation, and localization (Ding et al., 2012;Ding and Morrison, 2013;Poulos et al., 2013). Stromal cells make up another major constituent and multiple sub-populations have been previously described, including Lepr + MSPCs, NG2 + (CSPG4, a pericyte marker) cells, and CAR cells. NG2 + cells unsheathe arterioles have been previously proposed to maintain HSC quiescence (Kunisaki et al., 2013). Osteoblasts derived from the bone marrow MSPCs were previously shown to support early lymphoid progenitors (Ding and Morrison, 2013).
Recent studies extended single-cell transcriptomic profiling to the non-hematopoietic compartment of the murine bone marrow to further understand the combination of signals and their cellular source responsible for HSC maintenance and differentiation. As the frequencies of bone marrow cell types vary by orders of magnitude (Gomariz et al., 2018), some groups utilized an approach of depleting abundant hematopoietic cell types to isolate the niche populations in an unbiased manner while others devised a positive selection approach to focus on previously identified subsets ( Table 1).
Our group took advantage of lineage-specific Cre-labeling to profile 9,622 cells representing major niche subsets previously shown to play a role in hematopoiesis using the 10× Genomics platform (Tikhonova et al., 2019). Clustering analysis detected two endothelial, four perivascular, and three osteo-lineage subpopulations, validated using orthogonal approaches such as immunofluorescence or flow cytometry-based on the identified biomarkers. The vascular endothelial subset marked by expression of Cdh5 (VE-Cadherin) contained cells representing Stab2 + and Flt4 + sinusoidal capillaries and a smaller (12% of ECs) fraction of Ly6a + and Cd34 + arterial cells. The Lepr + perivascular compartment consisted primarily of adipo-primed mesenchymal cells characterized by high expression of Cxcl12 and adipogenic markers such as Adipoq and Lpl as well as osteo-primed cells defined by gradually increasing expression of osteogenic genes such as Bglap and Spp1. This was a surprising finding, given that the bone of young animals contains few adipocytes. The osteo-lineage subset was comprised of three transcriptionally distinct clusters representing osteoblasts, chondrocytes, and fibroblasts. To infer the relationship between the mesenchymal lineage cells, the developmental trajectory was reconstructed using Monocle (Trapnell et al., 2014). Pseudotime ordering of the perivascular and terminally differentiated osteo-lineage cells revealed a continuum of cellular states with progressively increasing levels of expression of adipogenic markers toward one endpoint and osteogenic markers toward the opposite end. The adipo-primed sub-populations were enriched for previously identified human bone marrow mesenchymal stem cell gene signature (Ghazanfari et al., 2016), and functional studies showed higher activity of fibroblastic colony-forming units. In contrast to osteo-rimed progenitors that were tdTomato + but expressed lower levels of Lepr, adipoprimed progenitors were both high in Lepr expression and tdTomato + . Taken together these findings indicate that early mesenchymal progenitors reside within adipo-gene expressing Lepr + MSPC fraction. Wolock et al. (2019) examined the mesenchymal lineage cells to determine the relations and hierarchies of stromal progenitors as they differentiate toward the mature bone, fat, and cartilage cells. 2,847 sorted non-hematopoietic (CD45 − Ter119 − ) and non-endothelial (CD31 − ) cells were profiled by 3 droplet-based inDrop scRNA-seq with a median of 736 detected genes per Frontiers in Cell and Developmental Biology | www.frontiersin.org a cell. Spectral clustering identified 7 clusters labeled as MSC, adipocyte progenitors, pre-adipocytes, osteoblast/chondrocyte progenitors, pre-osteoblast/chondrocytes, pro-osteoblasts, and pro-chondrocytes. SPRING visualization exposed a transcriptional continuum across the adipogenic and osteogenic branches. Differentiation trajectories were inferred with Velocyto (La Manno et al., 2018), a method to quantify the relationship between the abundance of precursor and mature mRNA, predicting individual cells' future state. RNA velocity analysis suggested that the MSC subpopulation was atop the differentiation hierarchy with branches toward the pre-adipocyte, pro-osteoblast, and pro-chondrocyte endpoints. PBA was utilized as an independent methodology to resolve each of the three lineages' average transcriptional trajectory and order genes based on their expression pattern along those trajectories.
Utilizing a negative selection approach, Baryawno et al. (2019) profiled 20,896 bone marrow Lin-cells and identified 17 clusters originating from bone and bone marrow, including endothelial cells, MSCs, osteolineage cells, chondrocytes, fibroblasts, and pericytes. Multiple methods were utilized to resolve differentiation relationships between the populations, including the correlation of expression profiles between clusters and diffusion map analyses. The clusters' connectivity was evaluated with partition-based graph abstraction (PAGA) that generates a graph-like map of data maintaining both its continuous and disconnected structure. Three endothelial clusters were identified, spanning a transcriptional continuum from Flt4 + sinusoidal to Ly6a + arteriolar sub-populations. The latter additionally had a sub-population of cells from endosteal and bone marrow arteries, marked by expression of Vwf and Kitl. Mesenchymal cells, classified based on the expression of Lepr and Cxcl12, were further divided into four subsets with varying levels of expression of those markers. One was characterized by osteo-lineage markers Sp7 and Alpl, indicating differentiation toward that lineage. This relationship was further profiled by analyzing the two osteolineage clusters, one covering a range of differentiation states and a transcriptionally distinct subset of mostly cells from the bone fraction. An additional distinct subpopulation of perivascular mesenchymal stromal cells and pericytes was identified with low levels of Lepr but high levels of Nes and NG2 (Cspg4) as well as pericytes markers Acta2, Myh11, and Mcam. Baccin et al. (2020) combined single-cell with spatially resolved transcriptomics. To compensate for the highly variable abundance of cell types within the bone marrow, undigested bone marrow or enzymatically digested bones were depleted or enriched for populations of interest. The final dataset included 7,497 cells grouped into 32 clusters. The total bone marrow yielded hematopoietic populations. Those populations' depletion resulted in primarily erythroid progenitors with low expression of Cd45 and 2% non-hematopoietic cells. These rare populations were captured by depleting cells expressing CD45 or an erythroid marker CD71. The identified non-hematopoietic populations included Schwann cells, smooth muscle cells, myofibroblasts, Pdgfra + mesenchymal cells, and endothelial cells. The endothelial subset was split into arterial and sinusoidal cells. The mesenchymal subset comprised nine sub-populations, including chondrocytes, osteoblasts, three fibroblast clusters, Ng2 + Nestin + MSCs, and two populations similar to Cxcl12abundant reticular (CAR) cells. The two CAR sub-populations differed primarily in their expression of adipocyte (Adipoq) and osteo-lineage (Sp7, Bglap) genes. RNA velocity analysis of the mesenchymal lineage cells assigned the Ng2 + MSCs atop the differentiation hierarchy. Baccin et al. (2020) additionally developed a laser-capture microdissection protocol coupled with sequencing (LCM-seq) to overcome spatial transcriptomic approaches that have not been successful in the bone marrow due to dependence on high RNA quality or unfixed tissue material. The method yielded full-length transcriptomic data from bone marrow sections containing 200-300 cells in a layer. LCM-seq was applied to 76 micro-dissected regions from the diaphyseal bone marrow to localize the bone marrow populations to endosteal, sinusoidal, arteriolar, and non-vascular niches. The frequencies of scRNA-seq populations within the spatially resolved transcriptomic data were estimated with CIBERSORT to determine the bone marrow cell types' localization to candidate niches. Osteoblasts and chondrocytes localized to the endosteal niches. Arterial endothelial cells and smooth muscle cells were found in the arteriolar niches. Sinusoidal cells localized to the sinusoidal niches but were also detected in the (sub)-endosteal niches. Adipo-CAR cells were associated with high sinusoidal areas, but the Osteo-CAR cells were found in arteriolar or non-vascular niches. Additionally, to predict these cell types' spatial relationships exclusively based on scRNA-seq, the RNA-Magnet method was developed, taking advantage of known cellular adhesion receptors and cognate ligands. The inferred adhesiveness of the identified populations to distinct niches correlated with the localization as measured by spatial transcriptomics. RNA-Magnet was also applied to signaling mediators, such as cytokines and growth factors, expressed by all cells in the dataset. The analysis suggested that the HSPCs were more likely to receive signals from non-hematopoietic populations than the mature immune cells, indicating a shift from mesenchymal to immune signaling. Both Adipo-CAR and Osteo-CAR cells were the primary sources of signals recognized by progenitor populations. Zhong et al. (2020) performed additional profiling of mesenchymal lineage cells in young, adult, and aging mice. The study used a Col2-Cre mouse model as a strategy to capture the most primitive mesenchymal lineage cells. 13,759 cells with a median of 2,686 detected genes were captured from young (1-1.5-month-old) animals. Clustering analysis identified 22 populations, of which 9 were mesenchymal, 11 hematopoietic, 1 endothelial, and 1 mural. The 7,585 mesenchymal linage cells contained osteoblast, osteocyte, adipocyte, and chondrocyte populations as well as four progenitor subsets. Early mesenchymal progenitors (EMPs) were determined as the most primitive cells in the dataset based on the expression of stemness markers Sca1, Cd34, and Thy1. Intermediate mesenchymal progenitors and late mesenchymal progenitors expressed progressively higher levels osteo-lineage genes such as Sp7 and Col1a1. Differentiation trajectory inference with Monocle placed EMPs at the terminal end of one of the three branches, with terminally differentiated osteocytes and adipocytes at the other ends. Lineage committed progenitors were placed at the point of bifurcation. Pseudotemporal reconstruction was additionally performed with Slingshot, which builds minimum spanning trees from clusters as opposed to single cells with Monocle. UMAP cell embeddings were used for trajectory inference, resulting in an analogous branching structure.

INTEGRATED OVERVIEW OF THE BONE MARROW NICHE DATASETS
The single-cell studies have been able to reveal extensive heterogeneity in the populations thought to be homogeneous. One of the challenges associated with the rapid flow of this new information is an understanding of how the populations described in these studies are related. Furthermore, inconsistent nomenclature further complicates one's ability to compare and contrast the data. To address this gap, we reanalyzed the five discussed scRNA-seq datasets combined (Table 1). We calculated the cluster markers for the populations identified within each dataset and performed hierarchical clustering based on the scaled expression correlation distances of those genes ( Figure 1A). To better assess the heterogeneity of the described populations and their transcriptional relationships, we additionally performed a joint analysis of 32,743 cells from the three datasets encompassing the most diverse subsets. Using the Seurat anchor-based data integration pipeline (Stuart et al., 2019), we were able to overlay them and visualize the populations using a consistent representation ( Figure 1B). Despite distinct sample preparation strategies (enzyme cocktails, digestion time, sorting and sorting strategies), assayed fraction (bone marrow or a combination of bone marrow and bone), and computational analysis approaches (summarized in Table 2), the findings are largely in agreement among the experiments carried out by independent groups. Expectantly, the relative abundances of various subpopulations fluctuated across the datasets (Figure 1C). For example, Baryawno dataset contained a large fraction of chondrocytes and fibroblasts, whereas Baccin subset uniquely captured myofibroblasts and Schwann cells.
Based on the expression patterns of cluster markers and visualization of integrated datasets, we harmonized the various sub-populations identified in the original studies using a consistent nomenclature (Figure 1B). Bone marrow vascular ECs were represented by sinusoidal (blue) and arteriolar (orange). By sampling a higher number of cells, Baryawno et al. identified an additional smaller sub-population of cells from endosteal and bone marrow arteries (dark orange). The mesenchymal stem cells give rise to a diverse range of cells such as adipocytes, osteoblasts, and chondrocytes. Previous studies took advantage of either Cxcl12 or Lepr to isolate perivascular cells, but the extent of the overlap of the two populations was not clear. The scRNAseq work revealed a broad transcriptional heterogeneity within these subsets. Bone marrow mesenchymal progenitors split into adipo-MSPCs (red) and osteo-MS (brown). Adipogenesisprimed mesenchymal cells were characterized by high expression of both Lepr and Cxcl12. These key markers were gradually decreased as the cells progressed toward osteogenesis-primed subpopulations. Osteoblasts (olive), chondrocytes (cyan) and a variety of fibroblast populations (gray) were represented across different dataset. Single cell studies are limited by the efficiency of digestion protocols and rare populations may not be sufficiently sampled. A more granular deconvolution of smaller bone subsets will require focused functional studies. To resolve transcriptional relationships between distinctly labeled clusters, the readers can refer to Figure 1A. To confirm that the gene expression of the integrated dataset corresponded to the populations identified in the individual studies, we assessed most differentially expressed marker genes (Figures 2A,B) and mapped expression patters of key pro-hematopoietic factors in the bone (Figure 2C). Indeed, we found the expression of Cxcl12 and Kitl to be enriched in adipo-MSCPs progenitors and arteriolar ECs, while Il7 and Il34 were expressed primarily in adipo-MSPCs ( Figure 2C). Collectively, our brief survey highlights a concordance between five different bone marrow niche single-cell datasets and resolves inconsistent nomenclature offered by distinct scientific groups.

ZOOMING INTO THE VASCULAR NICHES
Vascular ECs comprise a highly organized network throughout the body, facilitating the distribution of oxygen and nutrients. In addition to their delivery function, ECs enable immune cell trafficking allowing for cell migration and immune surveillance. Vascular ECs display tissue-specific morphology and have been shown in some cases to carry organ-specific functions.

Operation
Tool/Algorithm Description Normalization, batch correction, clustering Seurat Stuart et al., 2019) An R package for quality filtering, normalization, dimensionality reduction, and visualization of scRNA-seq data. It additionally includes a method for integrated analysis of multiple datasets by identifying pairwise correspondences between single cells across those datasets.

Visualization t-SNE (van der Maaten and Hinton, 2008)
Non-linear dimensionality reduction technique based on Student-t distribution for converting data in a high-dimensional space to a low-dimensional one while avoiding overcrowding.
SPRING (Weinreb et al., 2018a) Force-directed graph layout for visualizing continuous topologies that generates a graph of nodes representing cells connected to their nearest neighbors in high-dimensional gene expression space.
UMAP (Becht et al., 2018;McInnes et al., 2018) Approximates a manifold and a constructs its fuzzy topological representation with the goal of preserving more of the global structure.
Pseudotime or trajectory inference DPT/Destiny (Haghverdi et al., 2015;Angerer et al., 2016) Uses diffusion maps to identify the low-dimensional structure, then identifies a pseudotime metric based on transition probabilities of differentiating toward various cell fates.
Monocle (Trapnell et al., 2014;Qiu et al., 2017) Constructs a minimum-spanning tree (MST) through the dimension-reduced space created by independent component analysis. Cells are ordered along the longest path through the MST. Monocle was later modified to use DDRTree for dimensionality reduction and ordering.
PAGA (Wolf et al., 2019) A partition-based graph abstraction tool that provides a coarse-grained representation by placing edges between cluster nodes with similar cells. Unlike many trajectory inference methods, it is able to account for disconnected topologies.
PBA (Tusi et al., 2018;Weinreb et al., 2018b) Uses nearest neighbor graph cell densities to predict fate probabilities and the direction of differentiation.
Slingshot (Street et al., 2018) Uses a cluster-based MST to identify the lineages and where they branch, then uses simultaneous principal curves to fit a smooth representations of each lineage.
STEMNET (Velten et al., 2017) Uses hierarchical clustering to define the most differentiated cell populations and then uses those populations as a training set for classifying priming in the less mature populations.

Velocity
Velocyto (La Manno et al., 2018) Predicts the future state of single cells based on the relative abundance of unspliced precursor and spliced mature mRNA.
Cell-cell interactions RNA-Magnet (Baccin et al., 2020) Estimates the likelihood of physical interactions between single cells based on the expression of cell-surface receptors and their binding partners.
Frontiers in Cell and Developmental Biology | www.frontiersin.org Blood-brain barrier and kidney glomeruli vasculature have been areas of active research, but the specifics of vascular cells in other tissues are less understood. High-throughput transcriptional profiling studies compared molecular architecture of ECs across various organs, revealing a marked transcriptional heterogeneity (Chi et al., 2003;Nolan et al., 2013). In addition to variability driven by tissues, each organ contains several vessel types, including arteries, veins, and capillaries. Recent single-cell RNAseq analysis of > 32,000 ECs isolated from 11 different mouse tissues at the steady-state showed that vascular heterogeneity is largely dictated by the anatomical location of ECs rather than the vessel type (Kalucka et al., 2020), suggesting tissues-mediated molding of the vasculature and tissues specific functions.
In the bone marrow, ECs serve multiple roles, including HSC maintenance and leukocyte trafficking. Highly branched sinusoidal capillaries (SECs) carry venous blood and make up the vast majority of vascular cells. Arteriole vessels (AECs) carry arterial blood into the bone, have few branches, and constitute around ∼10% of the total bone marrow vasculature (Xu et al., 2018;Tikhonova et al., 2019). Here, we specifically focused on the vascular subset from the combined single-cell bone marrow niche datasets (Figure 3A). We calculated module scores for ECs based on the sinusoidal and arterial gene expression signatures obtained from the bulk tissue analysis (Xu et al., 2018) and found that 67.3% of cells were enriched for sinusoidal and 23.5% for the arterial signatures ( Figure 3B). To further confirm the identity of the cells, we assessed the expression levels of population marker genes and found Flt4, Tfpi, Stab2, and Vcam1 to be highly specific to sinusoids ( Figure 3C) and Cxl12, Kitl, Gja4, Vegfc to arterioles ( Figure 3D). Previous studies indicated that sinusoids and arterioles possess different permeability properties . Notably, hematopoietic stem cell adhesion, rolling and trans-endothelial events were reported exclusively in SECs . Indeed, sinusoids are enriched for expression of adhesion molecules such as E-selectin (Sele), Vcam-1 (Vcam1), ICAM-1 (Icam1), while arterioles expressed ICAM-2 (Icam2) (Figure 3E), consistent with the idea of differential transendothelial migration. It is not clear if all sinusoids can support trans-endothelial migration and it is tempting to speculate that there might be specialized sinusoid subtypes facilitating leukocyte trafficking.
Within the bone marrow, transitional zone vessels connect bone marrow sinusoids and arterioles. Based on the cell-surface marker expression and phenotypic characteristics, Kasumbe et al. identified two subtypes of SECs: CD31 low EMCN low L-vessels form highly branched capillary networks throughout the bone marrow and represent the majority of the SECs. CD31 high EMCN high H-vessels were found in the metaphysis, organized as vessel columns, and in transitional zones bridging arterioles with the sinusoids . H-vessels are of a substantial clinical interest, as they were shown to regulate osteogenesis Ramasamy et al., 2014). We were not able to clearly identify transitional vessels with our own transcriptional dataset (Tikhonova et al., 2019). However, exploring the combined analysis based on the three studies allowed us to pinpoint a subset of ECs that were not enriched for either arterial or sinusoidal signatures. Further examination revealed that these transition cells expressed higher levels of Emcn and Cd31 compared to sinusoids, suggesting that they might represent transitional H vessels ( Figure 3F). We found this population to be enriched for arterial-associated genes Cd34 and Ly6c1 but also express unique markers such as Cotl1 and Sox4 ( Figure 3G). Collectively, this integrated analysis allows for deeper exploration of BM EC heterogeneity. Although transitional EC population was not clearly revealed by each individual analysis, combining the datasets allowed for their further exploration, underscoring the utility of publicly available data to the community.
A recent study reports even further specification among bone marrow SECs. Chen and Liu et al. demonstrated that Apelin + (Apln) bone marrow ECs control vascular regeneration following irradiation-induced vascular damage. These specialized SECs make up approximately 0.03% of the total bone marrow and are located in the metaphysis, endosteum, and the transition area between the metaphysis and diaphysis. Bulk RNA-seq of Apln + cells revealed an enrichment of pathways related to angiogenesis and vascular endothelial growth factor (VEGF) signaling. Importantly, irradiation triggers the expansion of Apln + ECs and vascular regeneration. Spatial organization of hematopoietic differentiation is another open question. Elegant work by Zhang et al. revealed that following acute systemic infection, monocyte-dendritic cell progenitors localize to a subset of blood vessels expressing colony-stimulating factor 1 (Csf1), known to be a key regulator of myelopoiesis. Indeed, EC-specific deletion of Csf1, results in a loss of non-classical monocytes and dendritic cells during steady-state and following infection (Zhang et al., 2021). It remains to be elucidated if specialized ECs also support progenitors of lymphoid and erythroid lineages. Further mining of single-cell bone marrow niche datasets, single-cell sequencing of ECs following stress scenarios, combined with receptor-ligand interaction analysis (Browaeys et al., 2019;Cabello-Aguilar et al., 2020;Efremova et al., 2020) and clever functional studies will reveal further specialization of vascular ECs.
The precise identity of the vascular cells that support HSCs continues to be debated. Combining computational modeling with whole-mount confocal immunofluorescence imaging techniques, Kunisaki et al. found quiescent HSCs to be closely associated with AECs (Kunisaki et al., 2013). Optically clearing studies of the bone marrow, that allowed to perform deep confocal imaging, demonstrate that independently of their cycling status, approximately 80% of HSCs are closely associated with SECs, with another 10% of HSCs being adjacent to AECs, and 10% being adjacent to transition zone vessels (Acar et al., 2015). In contrast, Kunisaki et al., using whole-mount confocal immunofluorescence imaging, showed that quiescent HSCs associate specifically with small arterioles that are preferentially found in endosteal bone marrow. Furthermore, Itkin et al. demonstrated that arterioles maintain hematopoietic stem cells in a low reactive species oxygen state , perpetuating the debate regarding the precise identity of vascular cells that support quiescent HSCs. Just as novel technologies have been driving forward our understanding of hematopoietic differentiation, further advances combining multi-modal singlecell approaches, such as pairing gene expression data with spatial tissue context, and in vivo imaging will shape our understanding of hematopoietic cell interactions with their niches and resolve standing questions.

CONCLUSION
In a span of less than two years, our understanding of the hematopoietic differentiation and bone marrow microenvironment has made a major leap forward. Despite these advances, key questions remain unanswered (Tikhonova et al., 2020). For example, cellular identity of the bone marrow mesenchymal stem cells remains to be elucidated. Despite highly comparable data, pseudo-temporal ordering revealed conflicting differential MSPC hierarchies and led to distinct conclusions. We believe that scRNA-seq analyses are able to guide our understanding of MSPC biology and their exploration can inspire creativity in our scientific inquiry. However, they cannot replace functional experiments required to answer complex biological questions.
Single-cell profiling of hematopoietic hierarchy highlighted a fluidic nature of differentiation. Although we understand some of the intrinsic signals implementing these decisions, we do not yet grasp which extrinsic factors mediate hematopoietic differentiation. We speculate that specific temporospatial combinations of bone marrow niche signals guide hematopoietic progenitors toward gradual commitment to a lineage choice.