Identification of phenotypically, functionally, and anatomically distinct stromal niche populations in human bone marrow based on single-cell RNA sequencing
Abstract
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.
Editor's evaluation
The manuscript by Li and coworkers is a landmark characterization of sorted human non-hematopoietic 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.
https://doi.org/10.7554/eLife.81656.sa0Introduction
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.
Our understanding of the BM hematopoietic microenvironment (HME) has evolved considerably over the past decade through a number of landmark studies that have identified the cellular identity, anatomy, and functions of different murine HME components (Baryawno et al., 2019; Tikhonova et al., 2019; Wolock et al., 2019; Baccin et al., 2020). Skeletal and stromal cell populations, endothelial cells and other non-hematopoietic cells but also differentiated hematopoietic cells have been described as important HME elements in mice (Méndez-Ferrer et al., 2008; Urbieta et al., 2010; Yamazaki et al., 2011; Ding et al., 2012; Casanova-Acebes et al., 2013; Zhao et al., 2014; Chan et al., 2018).
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
View detailed protocolHuman 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 (Li et al., 2014). 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)
View detailed protocolFor 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 fluorescence-minus-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)
Request a detailed protocolUnsorted and FACS-sorted primary BM-MNC were cultured at plating densities of 5×104 cells/cm2 cells and 10–50 cells/cm2, 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).
In vitro differentiation assays
Request a detailed protocolCultured mesenchymal stromal cells (cMSCs) were differentiated towards the adipogenic, osteoblastic, and chondrogenic lineage as described previously (Li et al., 2014). Briefly, cells were cultured for 14 days in AdipoDiff medium (Miltenyi) and cells were stained with Oil Red O (Sigma). For osteogenic differentiation, cells were cultured in osteogenesis induction medium for 21 days and calcium depositions were detected by Alizarin Red staining (Sigma). Osteogenesis induction medium: standard MSC medium supplemented with 0.05 mM L-ascorbic-acid-2-phosphate (Wako Chemicals, Neuss, Germany), 0.1 µM dexamethasone and 10 mM β-glycerophosphate (both from Sigma).
Chondrogenic differentiation was induced by culturing cell pellets (2.5x105 cells/pellet) for 28 days in chondrogenesis induction medium. Chondrogenesis induction medium: DMEM-high glucose supplemented with 0.1 µM dexamethasone, 1 mM sodium pyruvate, 0.35 mM L-proline (all from Sigma), 0.17 mM ascorbic acid, 1% ITS +culture supplements (BD Biosciences) and 0.01 µg/ml TGF-ß3 (R&D Systems). Pellets were paraformaldehyde (PFA)-fixed and frozen in O.C.T. Compound (Sakura, Zoeterwoude, Netherlands). Cryosections were stained with Alcian Blue, and nuclei were stained with Nuclear Fast Red (both from Sigma). Sections were analyzed with a Nikon Eclipse TS100 microscope equipped with a Nikon DS-L3 digital camera.
Antibodies
The following antibodies were used for FACS analysis and cell sorting: CD235a-PE-Cy5 [clone GA-R2 (HIR2)], CD71- PE-Cy5 (clone M-A712), CD45-APC-Cy7 (clone 2D1), CD52-Alexa Fluor (AF) 488 (clone 4C8), CD56 (NCAM1)- Brilliant Violet 605 (clone Leu-19, NKH1), CD81-PE (JS-81) (all BD Bioscience), LEPR-PE (clone 52263, R&D Systems) CD271-APC (clone ME20.4–1 .H4, Miltenyi).
Immunofluorescence staining
Request a detailed protocolHuman bone marrow paraffin sections (5 µm) were deparaffinized and rehydrated following standard protocols. Following antigen retrieval with antigen retrieval solution (DAKO, S1699) and blocking with goat serum (Invitrogen, #31873), slides were stained sequentially scanned, and restained following fluorescence bleaching of primary coupled antibodies in PBS with 50 mM Na2CO3 and 3% H2O2. An OlympusVS120 slide scanner was used for scanning. Antibodies used were CD45-AlexaFlour 647 (Bio-Rad, MCA87A647T), CD56 (NCAM1)-AF647 (Biolegend, 318313), mouse anti-human CD271 (R&D Systems, MAB367), and rabbit anti-human CD81 (Novus Biologicals, NBP1-77039) with secondary goat anti-mouse AF647 (Jackson ImmunoResearch, 115-605-166) and goat anti-rabbit AF488 (Jackson ImmunoResearch, 111-545-003) with DAPI (Sigma, D956410MG) as a nuclear counterstain, respectively.
For SPP1 staining, after deparaffinization, rehydration, and antigen retrieval, slides were permeabilized with 0.1% Triton X100 (Sigma Ultra, T9284) solution at 37 °C for 15 min, blocked using donkey serum (Jackson Immuno Research, 017-000-121), and stained overnight at 4 °C with goat anti-human SPP1 antibody (R&D Systems, AF1433-SP), followed by secondary staining with donkey anti-goat AF594 (Jackson ImmunoResearch, 706-585-147) and staining with CD56 (NCAM1)-AF647 plus DAPI. Slides were then scanned, bleached, and incubated with rabbit anti-human CD45 (Sigma-Aldrich, HPA000440) and mouse anti-human CD271 overnight, followed by secondary antibody staining with goat anti-mouse AF647 and goat anti-rabbit AF488 and scanned again.
For confocal laser scanning, slides went through the same standard deparaffinization, rehydration, and antigen retrieval procedure. For all primary antibody stainings, blocking was done with normal goat serum, followed by overnight incubation at 4 °C with rabbit and mouse primary antibodies. Thereafter, slides were incubated with secondary goat anti-mouse AF647 and goat anti-rabbit AF488 for 1 hr at room temperature (RT). Primary antibodies used were mouse anti-human CD271 (R&D Systems, MAB367), rabbit anti-human CD81 (Novus Biologicals, NBP1-77039), mouse anti-human CD81 (Abcam, ab79559), rabbit anti-human CD45 (Sigma-Aldrich, HPA000440), rabbit anti-human CD56 (Sigma-Aldrich, 156 R-96), and mouse anti-human CD56 (Fisher scientific, 10027562). For SPP1, NCAM1, CD271 co-stainings, slides were first blocked with goat serum and incubated with rabbit anti-human NCAM1 and mouse anti-human CD271 antibodies overnight at 4 °C. Thereafter, slides were incubated with secondary goat anti-mouse AF647 and goat anti-rabbit AF488 for 1 hr at RT followed by permeabilization with 0.1% Triton X100 (Sigma Ultra, T9284) solution at 37 °C for 15 min. After blocking with donkey serum donkey anti-human SPP1 antibody was incubated overnight at 4 °C followed by donkey anti-goat AF594 (Jackson ImmunoResearch, 706-585-147) incubation for 1 hr at RT.
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
Request a detailed protocolAll 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.
Co-culture of BM CD34+ cells with stromal cells
Request a detailed protocolFACS-sorted stromal cells (OC: CD45low/-CD235a-CD71-CD271+CD52-NCAM1+; MSSC: CD45low/-CD235a-CD71-CD271+CD52-NCAM1-CD81++) were culture expanded and plated as adherent feeder cells into 96-well plates at 10,000 cells per well and cultured in StemMACS MSC Expansion medium for 1 day. Then, medium was removed and 5000 BM CD34+ cells were added in serum-free expansion medium (SFEM, STEMCELL Technologies) supplemented with or without stem cell factor, thrombopoietin, and FLT-3 ligand (all at 25 ng/ml). Expanded cells were harvested, counted and the expression of CD34 and CD90 was analyzed by flow cytometry after 7 days of co-culture.
Single-cell RNA sequencing and data analysis
Single-cell RNA sequencing
Request a detailed protocolSingle-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/-CD235a- cells 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/-CD235a- and 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
Request a detailed protocolDimensionality 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
Request a detailed protocolGraph-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
Request a detailed protocolFor 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
Request a detailed protocolBone 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
Request a detailed protocolCell-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
Request a detailed protocolThe 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.
Results
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 non-hematopoietic cell-containing CD45low/-CD235a- population 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 CD45low/-CD235a- population (Figure 1—figure supplement 1A). In order to also capture and be able to analyze extremely rare stromal cell subsets, we therefore also sorted CD45low/-CD235a-CD271+ cells, which are highly enriched for CXCL12-expressing BM stromal stem and progenitor cells (Figure 1A, Figure 1—figure supplement 1B; Tormin et al., 2011; Li et al., 2014). As CD45low/-CD235a-CD271+ cells represent a subset of CD45low/-CD235a- cells, 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, Figure 1—figure supplement 1C) which formed 42 clusters corresponding to distinct cell types and differentiation stages, respectively (Figure 1B, Figure 1—figure supplement 1D–F, Supplementary file 1; Supplementary file 2; Supplementary file 3). As expected (Boulais et al., 2018; Baryawno et al., 2019), scRNAseq analysis of CD45low/-CD235a- cells captured several hematopoietic cell types including a cluster of hematopoietic stem/progenitor cells (HSPCs) enriched for CD34 expression and other HSPC markers (PROM1, CRHBP, AVP, MLLT3, FAM30A, GATA1) (Figure 1—figure supplement 1G, and Supplementary file 2 and Supplementary file 4). Importantly, non-hematopoietic cells including stromal cells as well as non-stromal cells, such as endothelial cells (cluster 28), KRT5-enriched basal cell-like cluster (cluster 0), and a cluster of cells enriched for neuronal markers expression (cluster 39) could be clearly identified (Figure 1C and D and Figure 1—figure supplement 1H). These clusters were detected in both CD45low/-CD235a- and CD45low/-CD235a-CD271+ cells (Figure 1—figure supplement 1D–E). 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 CD45low/-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 CD45low/-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).
The stromal clusters were then annotated based on their expression of previously reported BM stromal markers (Churchman et al., 2012; Li et al., 2014), which allowed us to discriminate six different cell types, that is multipotent stromal stem cells (MSSCs), highly adipocytic gene-expressing progenitors (HAGEPs), balanced progenitors, pre-osteoblasts, osteochondrogenic progenitors (OCs), and pre-fibroblasts (Figure 2A–B, Figure 2—figure supplement 1A–B).
Cluster 3 was annotated as MSSCs with multi-lineage differentiation capacity based on its high expression of stromal markers (CXCL12, LEPR, DCN, PTGDS) and genes indicating both adipogenic (CEBPD, LPL, PLIN1, ADIPOQ, CCL2, PPARG) and osteogenic (GAS6, FBN1, ALPL, RUNX1, SPARCL1, CDH11) differentiation capacities (Figure 2C and D). The annotation of cluster 3 as MSSCs was further confirmed by velocity analysis as shown in Figure 3A.
Cells in cluster 5 expressed adipogenic differentiation markers as well as a group of stress-related transcription factors such as FOS, FOSB, JUNB and EGR1, some of which have been shown to mark bone marrow adipogenic progenitors (Figure 2B–C and Figure 2—figure supplement 1C; 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 and Leś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 pre-osteoblast 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 (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 group A cells, which could furthermore be separated into different subpopulations based on the differential expression of CD81.
Based on this, we designed FACS strategies for the isolation of group A (CD45low/-CD235a-CD71-CD271+NCAM1-CD52-), group B (CD45low/-CD235a-CD71-CD271+NCAM1+CD52-) and group C (CD45low/-CD235a-CD71-CD271+NCAM1-CD52+) cells (Figure 4B). As CD71 (TFRC) expression was barely detected in stromal clusters, CD71 was included as an additional exclusion marker to remove residual CD71-expressing erythroid cells which were present in the CD45low/-CD235a- population (Figure 1—figure supplement 1G). Group A cells (CD45low/-CD235a-CD71-CD271+NCAM1-CD52-) were further separated into four subpopulations (A1-A4) based on CD81 expression, which corresponded to MSSCs (A1: CD81++), HAGEPs (A2: CD81+), balanced progenitors (A3: CD81+/-), and pre-osteoblasts (A4: CD81-) (Figure 4B). Quantitative real-time PCR of sorted cells demonstrated that the expression levels of CD81 in A1 and A2 were higher than A3 and A4, indicating these four phenotypically defined subsets reflected the transcriptionally identified clusters in group A (Figure 4—figure supplement 1B).
Functional analysis of sorted cells showed that the MSSC population (A1: CD45low/-CD235a-CD71-CD271+NCAM1-CD52-CD81++) had the highest potential to form fibroblastic colonies (colony-forming units, fibroblast: CFU-F) and multi-lineage differentiation capacity toward osteoblasts, adipocytes, and chondrocytes (Figure 4C–D, Figure 4—figure supplement 1C). Although the CFU-F frequency of HAGEPs (A2: CD45low/-CD235a-CD71-CD271+NCAM1-CD52-CD81+) was lower compared to MSSCs, they still gave rise to a significant fraction of colonies and showed full in vitro differentiation capacities (Figure 4C–D, Figure 4—figure supplement 1C).
Colony-forming capacities of balanced progenitors (A3: CD45low/-CD235a-CD71-CD271+NCAM1-CD52-CD81+-), pre-osteoblasts (A4: CD45low/-CD235a-CD71-CD271+NCAM1-CD52-CD81-) and OCs (B: CD45low/-CD235a-CD71-CD271+NCAM1+CD52-) were considerably lower in comparison with MSSCs and HAGEPs (Figure 4C and Figure 4—figure supplement 1C). Pre-fibroblasts (group C) showed the lowest potential to form fibroblastic colonies (Figure 4C and Figure 4—figure supplement 1C). Furthermore, the adipogenic potential of pre-fibroblasts was markedly reduced and the pre-fibroblast-derived chondrocyte pellets were smaller compared to the other groups while osteogenic differentiation was normal (Figure 4D), thus indicating a limited differentiation potential of this population. OCs (group B) demonstrated an enhanced capacity to differentiate towards osteoblasts and chondrocytes, however, at the expense of adipogenic potential (Figure 4D). These findings suggest a skewed differentiation potential of OCs, which is consistent with the lower expression of adipogenic markers and elevated expression of osteochondrogenic markers (Figure 2C and E). The in vitro differentiation capacities of balanced progenitors (A3) and pre-osteoblasts (A4) could not be assayed due to insufficient proliferation upon culture expansion.
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 pre-fibroblasts. 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 Figure 5—figure supplement 1A). 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 inter-stromal 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 CXCL12-mediated 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 5—figure 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 the potential to support hematopoietic cells via direct cell-to-cell contact as well as hematopoiesis-supporting 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 CXCL12- and 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 SPP1-CD44 (Figures 5B and 6D and Figure 5—figure supplement 1B). 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, Figure 5—figure supplement 1B, Figure 5—figure supplement 2G). Furthermore, Oncostatin M (OSM) expression was found in various hematopoietic cell types, indicating a possible regulation of stromal cells through its receptors, OSMR and LIFR (Figure 5B, Figure 5—figure supplement 1B, Figure 5—figure supplement 2H).
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).
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 SPP1-mediated 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 CD45low 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 pre-fibroblasts), 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 part – the 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., 2014; 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 cell-to-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 likelihood-based 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).
Data availability
The scRNA-seq matrix data generated in this study have been deposited in the GEO database (GSE190965).
-
NCBI Gene Expression OmnibusID GSE190965. Transcriptomic profiling of human bone marrow non-hematopoietic cells.
References
-
Generalizing RNA velocity to transient cell states through dynamical modelingNature Biotechnology 38:1408–1414.https://doi.org/10.1038/s41587-020-0591-3
-
Regulated expression of three C/EBP isoforms during adipose conversion of 3T3-L1 cellsGenes & Development 5:1538–1552.https://doi.org/10.1101/gad.5.9.1538
-
Development of humanized ossicles: bridging the hematopoietic gapTrends in Molecular Medicine 26:552–569.https://doi.org/10.1016/j.molmed.2020.01.016
-
Neutrophil-Derived microvesicles enter cartilage and protect the joint in inflammatory arthritisScience Translational Medicine 7:315ra190.https://doi.org/10.1126/scitranslmed.aac5608
-
Isolation and characterization of bone marrow multipotential mesenchymal progenitor cellsArthritis and Rheumatism 46:3349–3360.https://doi.org/10.1002/art.10696
-
Isolation and characterization of primary bone marrow mesenchymal stromal cellsAnnals of the New York Academy of Sciences 1370:109–118.https://doi.org/10.1111/nyas.13102
-
Expression of versican and ADAMTS1, 4, and 5 during bone development in the rat mandible and hind limbThe Journal of Histochemistry and Cytochemistry 53:1553–1562.https://doi.org/10.1369/jhc.5A6669.2005
-
Pdgfrα and CD51 mark human nestin+ sphere-forming mesenchymal stem cells capable of hematopoietic progenitor cell expansionThe Journal of Experimental Medicine 210:1351–1367.https://doi.org/10.1084/jem.20122252
-
Haematopoietic stem cell activity and interactions with the nicheNature Reviews. Molecular Cell Biology 20:303–320.https://doi.org/10.1038/s41580-019-0103-9
-
SoftwareMSC_BM_scripts, version swh:1:rev:22c09a6f37f4cd46f977b3a97b9e8139776030bfSoftware Heritage.
-
Single-Cell RNA sequencing deconvolutes the in vivo heterogeneity of human bone marrow-derived mesenchymal stem cellsInternational Journal of Biological Sciences 17:4192–4206.https://doi.org/10.7150/ijbs.61950
Article and author information
Author details
Funding
Swedish Cancer Foundation (20-1163PjF 01H)
- Stefan Scheding
Swedish Childhood Cancer Foundation (PR2018-0078)
- Stefan Scheding
Swedish Childhood Cancer Foundation (PR2021-0065)
- Stefan Scheding
Swedish Bloodcancer Association
- Stefan Scheding
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank Helene Larsson and Anna Jonasson for their help to collect bone marrow samples and the Lund Stem Cell Center FACS facility personnel for technical assistance. We would also like to thank Dr. Paul Bourgine for critically reading the manuscript. Funding: This work was supported by funds from the StemTherapy Program, the Swedish Cancer Foundation, the Swedish Childhood Cancer Foundation, the Swedish Bloodcancer Association (Blodcancerförbundet), Foundation Siv-Inger and Per-Erik Anderssons minnesfond, John Persson Foundation, ALF (Government Public Health Grant), and the Skåne County Council Research Foundation.
Ethics
Human bone marrow (BM) cells were collected at the Hematology Department, Skåne University Hospital Lund, Sweden, from consenting healthy donors. The use of human samples was approved by the Regional Ethics Review Board in Lund, Sweden.
Copyright
© 2023, Li et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 4,212
- views
-
- 685
- downloads
-
- 15
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Stem Cells and Regenerative Medicine
Adult mammals, unlike some lower organisms, lack the ability to regenerate damaged hearts through cardiomyocytes (CMs) dedifferentiation into cells with regenerative capacity. Developing conditions to induce such naturally unavailable cells with potential to proliferate and differentiate into CMs, that is, regenerative cardiac cells (RCCs), in mammals will provide new insights and tools for heart regeneration research. In this study, we demonstrate that a two-compound combination, CHIR99021 and A-485 (2C), effectively induces RCCs from human embryonic stem cell-derived TNNT2+ CMs in vitro, as evidenced by lineage tracing experiments. Functional analysis shows that these RCCs express a broad spectrum of cardiogenesis genes and have the potential to differentiate into functional CMs, endothelial cells, and smooth muscle cells. Importantly, similar results were observed in neonatal rat CMs both in vitro and in vivo. Remarkably, administering 2C in adult mouse hearts significantly enhances survival and improves heart function post-myocardial infarction. Mechanistically, CHIR99021 is crucial for the transcriptional and epigenetic activation of genes essential for RCC development, while A-485 primarily suppresses H3K27Ac and particularly H3K9Ac in CMs. Their synergistic effect enhances these modifications on RCC genes, facilitating the transition from CMs to RCCs. Therefore, our findings demonstrate the feasibility and reveal the mechanisms of pharmacological induction of RCCs from endogenous CMs, which could offer a promising regenerative strategy to repair injured hearts.
-
- Stem Cells and Regenerative Medicine
Bone regeneration is mediated by skeletal stem/progenitor cells (SSPCs) that are mainly recruited from the periosteum after bone injury. The composition of the periosteum and the steps of SSPC activation and differentiation remain poorly understood. Here, we generated a single-nucleus atlas of the periosteum at steady state and of the fracture site during the early stages of bone repair (https://fracture-repair-atlas.cells.ucsc.edu). We identified periosteal SSPCs expressing stemness markers (Pi16 and Ly6a/SCA1) and responding to fracture by adopting an injury-induced fibrogenic cell (IIFC) fate, prior to undergoing osteogenesis or chondrogenesis. We identified distinct gene cores associated with IIFCs and their engagement into osteogenesis and chondrogenesis involving Notch, Wnt, and the circadian clock signaling, respectively. Finally, we show that IIFCs are the main source of paracrine signals in the fracture environment, suggesting a crucial paracrine role of this transient IIFC population during fracture healing. Overall, our study provides a complete temporal topography of the early stages of fracture healing and the dynamic response of periosteal SSPCs to injury, redefining our knowledge of bone regeneration.