Dissection of cellular communication between human primary osteoblasts and bone marrow mesenchymal stem cells in vivo at single-cell resolution

Osteoblasts are derived from bone marrow mesenchymal stem cells (BMMSCs) and play important role in bone remodeling. While our previous studies have investigated the cell subtypes and heterogeneity in osteoblasts and BMMSCs separately, cell-to-cell communications between osteoblasts and BMMSCs in vivo in humans have not been characterized. In this study, we performed a systematic integration analysis with our single-cell RNA sequencing (scRNA-seq) transcriptomes data from BMMSCs and osteoblasts. We successfully identied a novel preosteoblasts subtype which highly expressed ATF3, CCL2, CXCL2 and IRF1. Biological functional annotations of the transcriptomes suggested that the novel preosteoblasts subtype may inhibit osteoblasts differentiation, maintain cells to a less differentiated status and recruit osteoclasts. Ligand-receptor interaction analysis showed strong interaction between mature osteoblasts and BMMSCs. Meanwhile, we found FZD1 was highly expressed in BMMSCs of osteogenic differentiation direction. WIF1 and SFRP4, which were highly expressed in mature osteoblasts were reported to inhibit osteogenic differentiation. We speculated that WIF1 and sFRP4 expressed in mature osteoblasts inhibited the binding of FZD1 to Wnt ligand in BMMSCs, thereby further inhibiting osteogenic differentiation of BMMSCs.


Abstract Background
Osteoblasts are derived from bone marrow mesenchymal stem cells (BMMSCs) and play important role in bone remodeling. While our previous studies have investigated the cell subtypes and heterogeneity in osteoblasts and BMMSCs separately, cell-to-cell communications between osteoblasts and BMMSCs in vivo in humans have not been characterized.

Results
In this study, we performed a systematic integration analysis with our single-cell RNA sequencing (scRNAseq) transcriptomes data from BMMSCs and osteoblasts. We successfully identi ed a novel preosteoblasts subtype which highly expressed ATF3, CCL2, CXCL2 and IRF1. Biological functional annotations of the transcriptomes suggested that the novel preosteoblasts subtype may inhibit osteoblasts differentiation, maintain cells to a less differentiated status and recruit osteoclasts. Ligandreceptor interaction analysis showed strong interaction between mature osteoblasts and BMMSCs.
Meanwhile, we found FZD1 was highly expressed in BMMSCs of osteogenic differentiation direction. WIF1 and SFRP4, which were highly expressed in mature osteoblasts were reported to inhibit osteogenic differentiation. We speculated that WIF1 and sFRP4 expressed in mature osteoblasts inhibited the binding of FZD1 to Wnt ligand in BMMSCs, thereby further inhibiting osteogenic differentiation of BMMSCs.

Conclusions
At the single cell level, this study provided insights into the cell-to-cell communications between BMMSCs and osteoblasts and mature osteoblasts may mediate negative feedback regulation of osteogenesis process.

Background
Bone is a metabolically active organ that undergoes continuous remodeling throughout life [1,2]. Bone remodeling involves the bone resorption by osteoclasts followed by the formation of bone matrix through the osteoblasts. Osteoblasts are derived from bone marrow mesenchymal stem cells (BMMSCs). In addition to osteoblasts, BMMSCs can differentiate into adipocytes and chondrocytes. The adipogenic and osteoblastogenic differentiation are competitively balanced [3][4][5]. Osteoblastogenic development is a complex process that requires the tight regulation of speci c transcription factors activation or suppression in response to local signaling pathways [6,7]. For instance, Wnt signaling is implicated in cell fate determination, proliferation, differentiation and apoptosis of osteoblasts [6]. The combination of the Wnt and the seven-pass transmembrane receptor Frizzled (FZD) or single-pass transmembrane coreceptor complex LRP5/6 activates the canonical Wnt signaling pathway to promote downstream osteoblast target gene expression [7].
Cell heterogeneity is an important contributor to biological function. During the entire osteoblast development process from BMMSCs to mature osteoblasts, several cell subtypes have been identi ed. Two subpopulations of BMMSCs were found: a skeletal stem cell subpopulation that differentiated into bone, cartilage and fat, and another subpopulation that modulated immune function and in ammation [8]. Our previous studies performed single-cell sequencing of human BMMSCs and osteoblasts. In BMMSCs population, we identi ed osteoblast precursors, adipocyte precursors, chondrocyte precursors and terminal-stage cells subpopulations in LEPR high CD45 low BMMSCs [9]. In osteoblasts, we identi ed two known cell subtypes and a novel undetermined osteoblasts subtype. Meanwhile, our studies revealed the functional characteristics of each osteoblast subtype in bone metabolism, angiogenesis modulation and hematopoietic stem cells (HSC) niche regulation [10]. The osteoblast development from BMMSCs to mature osteoblast was a continuous process. But in our previous studies, BMMSCs and osteoblasts were analyzed separately, which neglected the continuity of differentiation and the cell to cell communication of different cell subtypes. Therefore, due to the lack of comprehensive combined analysis of BMMSCs and osteoblasts, the regulation mechanisms involved in osteoblast development process and cell to cell communication between BMMSCs and osteoblasts at single cell level remained limited.
In this study, we conducted an integrated analysis of BMMSCs and osteoblasts at the single-cell level, reconstructed pseudotime trajectories and analyzed cell to cell communication between BMMSCs and osteoblasts. We suggested that SFRP4 and WIF1 which highly expressed in mature osteoblasts, inhibited the binding of FZD1 to Wnt ligand in BMMSCs, thereby further inhibiting the regulation of osteogenic differentiation.

Results
Integration of single-cell transcriptome sequencing data from osteoblasts and BMMSCs The scRNA-seq data of human osteoblasts [10] and BMMSCs [9] from two subjects used in this study were detailed in our previous publication. After quality control and excluding contaminant cell types, 8,841 cells (5,976 osteoblasts and 2,865 BMMSCs) were used in the integration analysis. We performed unsupervised detection using principal component analysis in Seurat [11] and generated a twodimensional (2D) embedding plot for visualization using the uniform manifold approximation and projection (UMAP) algorithm [12]. 8 clusters (C1-C8) were identi ed based on highly variable genes ( Figure   1A). C1-C3 were BMMSCs (BMMSC1, BMMSC2, BMMSC3). C4-C6 were preosteoblasts1-3 clusters (PreOB1, PreOB2, PreOB3). C7 and C8 were immature osteoblasts cluster (ImmatOB) and mature osteoblasts cluster (MatOB) respectively. We observed that the two BMMSCs subjects did not blend well with each other after integration analysis ( Figure 1A), but the gene expression pro le between the two subjects was strongly correlated (R 2 = 0.96, Figure 1B). Hence the batch effect between the two subjects was relatively small and did not affect subsequent analysis.
Next, we identi ed highly expressed genes in each cluster ( Figure 1C) and used Gene Ontology (GO) enrichment analysis to characterize the most signi cant gene expression signatures from each cell cluster ( Figure 1D). Most of clusters were enriched for GO terms related to "ossi cation", "bone development", "osteoblast differentiation", "bone mineralization". Most of clusters expressed osteoblasts related marker genes (eg, ALPL, COL1A1, ITGB1) and mature osteoblasts marker genes (eg, IBSP, BGLAP, SPP1) ( Figure 1C). However, we found cluster BMMSC2 was enriched for terms related to "lipid localization", "lipid transport", and "lipid homeostasis" ( Figure 1D). Also, we found that LPL which was a marker gene expressed earlier in adipogenesis showed increased expression in both the directions of osteogenic (cluster BMMSC1) and adipogenic (cluster BMMSC2) differentiation, but adiponectin (ADIPOQ) which was a representative adipokine only expressed in cluster BMMSC2 ( Figure 1C). Therefore, we suggested that cluster BMMSC2 was a BMMSCs subcluster with adipogenic potential and function. Similarly, the cluster BMMSC3 from BMMSCs did not enrich for "ossi cation", "osteoblast differentiation" and other terms related to bone development ( Figure 1D), and did not express essential osteogenic genes, such as alkaline phosphatase (ALPL) and integrin subunit β1(ITGB1) for osteoblast development ( Figure 1C). Hence, as reported in our previous study, cluster BMMSC3 was at the terminal stage of BMMSCs and lacked the capacity for cellular differentiation [13]. Therefore, BMMSC1 was a BMMSCs subcluster with osteogenic potential and function. Cells from clusters BMMSC1, PreOB1, PreOB2, PreOB3, ImmatOB and MatOB were classi ed as osteogenesis cells and used in further analysis.
To assess if clustering was critically in uenced by the cell cycle, we calculated the cell cycle score to determine the stages of the cells, in which cell cycle scores were based on its expression of G2/M and S phase markers [11]. We found different cell cycle stages were present at a roughly similar proportion in each subpopulation ( Figure S3B-S3C), therefore, the cluster results were not affected by cell cycle heterogeneity.

Developmental trajectory of osteogenesis cells
Next, we inferred osteogenesis differentiation trajectory using unsupervised ordering of the scRNA-seq for osteogenesis cells by diffusion mapping (Figure 2A) [14]. By comparing the distribution of the pseudotime, we found that cluster BMMSC1was enriched in the early stage of pseudotime ( Figure 2B). The clusters ImmatOB and MatOB, which highly expressed speci c genes of late osteoblasts stage like IBSP, BGLAP and SPP1( Figure 1C,2B), were mostly presented at the later stage. Therefore, clusters ImmatOB and MatOB were at the terminal stages throughout the pseudotime. Next, we checked the transcriptional continuous expression of osteoblast characteristic genes during the process of osteogenesis differentiation. The osteogenesis cells speci c genes such as SPARC, COL1A1, BGLAP and SPP1 were upregulated from BMMSCs to mature osteoblasts ( Figure 2C). A subset of genes was consistently downregulated during the osteogenesis differentiation process, such as CXCL12 and LEPR ( Figure 2C). To support the trajectory inference, the expression of signi cant osteogenesis regulatory factors and transcription factors were mapped to osteogenic lineage. We observed that IBSP was expressed in mature osteoblasts as previously reported, and genes expressed in early osteogenic development such as LEPR and CXCL12 were not expressed in mature osteoblasts ( Figure 2D). This result was consistent with the ndings of transcriptional continuum in osteoblasts differentiation [15].

Transcriptional characteristics of human osteogenesis cells
Based on the functional enrichment analysis of the highly expressed genes, we found that the GO terms related to the Wnt/β-catenin pathway, such as "regulation of Wnt signaling pathway", "positive regulation of Wnt signaling pathway", "positive regulation of canonical Wnt signaling pathway", "non-canonical Wnt signaling pathway" were mainly enriched in BMMSC1 ( Figure 3A). Comparing the BMMSCs subclusters related to osteogenesis (BMMSC1) with adipogenesis (BMMSC2), we found that the expression of ADAM28 (ADAM metallopeptidase domain 28) was signi cantly upregulated in BMMSC1, but ADAM28 was almost absent in BMMSC2 ( Figure 3B). ADAM28 belongs to the ADAM family of multifunctional cell surface and secreted glycoproteins. ADAM28 was expressed in multiple human tooth germ and tooth mesenchyme [16]. Overexpression of ADAM28 favored to the proliferation and lineage-speci c differentiation of human dental papilla mesenchymal cells (hDPMCs) [17]. We speculated that ADAM28 played a similar role in BMMSCs and could induce BMMSCs to enter osteogenic differentiation instead of adipogenic differentiation.
Preosteoblasts, which were in the transitional stage from skeletal progenitors to osteoblasts, not only expressed BMMSCs speci c genes, such as leptin receptor (LEPR) and thymus cell antigen 1 (THY1, an important regulator of osteogenic differentiation and adipogenic differentiation), but also highly expressed matrix gla protein (MGP) and insulin-like growth factor-binding protein 4 (IGFBP4) ( Figure 3B). MGP promoted the bone formation by activating the Wnt/β-catenin pathway [18]. IGFBP4 reduced the proliferation of BMMSCs, thus inhibiting excessive osteogenic differentiation in rats [19,20].
Studies revealed that CXCL2 promoted osteoclastogenesis and prevented preosteoblast differentiation through inhibiting the ERK1/2 signaling pathway [28,29]. Interestingly, we found that immediate early response 3 (IER3) has the highest expression in the PreOB3 cluster ( Figure 3B), but there were few studies on IER3 in bone. So, we suggested that IER3 might be a crucial gene.
GO enrichment analysis showed that terms which were related to bone development in PreOB2 included "positive regulation of osteoblast differentiation", "regulation of bone mineralization", "extracellular structure organization", "extracellular matrix organization" etc. On the contrary, the GO terms related to "osteoclast differentiation" and "regulation of osteoclast differentiation" were enriched in PreOB3 ( Figure  3A). Combining pseudotime trajectory analysis and transcription factor annotation, we found that PreOB1 would transition to PreOB2 and PreOB3 after entering the end of development, but only PreOB2 could develop into mature osteoblasts. Cluster PreOB3 and PreOB2 had different effects on osteogenesis. PreOB3 could recruit osteoclasts and inhibit the proliferation and differentiation of preosteoblasts through expressing ATF3, CCL2 and CXCL2 genes.
Next, we checked the expression of mature osteoblasts signature genes in the integrated data. Except the expression of BGLAP, COL1A1 and SPP1, we also found that interferon-induced transmembrane protein 5 (IFITM5) with bone formation and osteoblast maturation functions [30] was highly expressed ( Figure  3B). In a study, Takao and colleagues demonstrated that cell adhesion molecule 1 (CADM1) expressed transiently during osteoblastic maturation in mouse. CADM1, located on the "lateral" membrane of osteoblasts, was mainly responsible for interstitial communication [31]. Our results showed that CADM1 was also highly expressed ( Figure 3B). Therefore, CADM1 might play a similar role in cellular communication in human mature osteoblasts. Previous studies showed that lymphocyte-speci c protein 1 (LSP1) was mainly expressed in white blood cells [32,33]. A recent study suggested endothelium also expressed LSP1 which played an essential role in transendothelial migration [33]. We found that LSP1 was highly enriched in immature osteoblasts and mature osteoblasts ( Figure 3B).

Ligand-receptor interactions in various stages of osteogenesis cells
To investigate the biological interaction between BMMSCs and osteoblasts, we performed intercellular communication analysis using an R Package "iTALK" [34], which can identify ligands and receptors in scRNA-seq data and analyze mediated cross-signal of cell communication ( Figure 4A). The top 30 pairs of ligand-receptor related to intercellular communication in the osteogenic lineage were visualized ( Figure  4B). We observed strong interaction between mature osteoblasts and BMMSCs. The results showed integrin subunit β1 (ITGB1) receptors were widely present in various subclusters ( Figure 4C). However, the prerequisite for integrin as a receptor was that the integrin β chain [35,36] must be combined with the integrin α chain to form heterodimeric transmembrane proteins [37] in order to bind the ligand. In our results, secreted phosphoprotein 1 (SPP1) and type I collagen (COL1A1) in mature osteoblast clusters could interact with ITGB1 receptor. SPP1 could inhibit adipogenesis and promote osteogenic differentiation of BMMSCs through integrin αv/β1 [38]. COL1A1 regulated the process of osteogenic differentiation which was mainly mediated by integrin α1/β1 [37,39,40]. The α chain integrin subunit alpha V (ITGAV) and integrin subunit alpha 1(ITGA1) were mainly expressed in BMMSC1, PreOB1, PreOB2 and PreOB3 ( Figure 4C). Therefore, although ITGB1 was expressed in all subclusters of osteogenic lineage, COL1A1 and SPP1 expressed in the mature osteoblasts may promote osteogenic differentiation of BMMSCs and preosteoblasts.
Discoidin domain receptor 2 (DDR2) and COL1A1 were another signi cant ligand-receptor pair. DDR2 was mostly expressed in BMMSC1 ( Figure 4C), and principally bonded to collagens type I, II, III and X [41,42]. After being activated by collagen ligands, the downstream signals propagated by DDR2 stimulated the ERK1/2 and p38 MAP kinase pathways to increase osteoblast differentiation [43,44]. Therefore, COL1A1 speci cally expressed by mature osteoblasts could mediate the differentiation of BMMSCs towards osteogenic direction.

Protein-protein interaction network
In order to identify signi cant genes that play a critical important role in mature osteoblasts and BMMSCs interactions. We mapped top 30 genes of the two clusters to a protein-protein interaction network based on the STRING dataset, and then used Cytoscape [45] to identify biologically meaningful genes and display their relationships. We found some interesting protein-protein interactions in mature osteoblasts (MatOB) and BMMSCs (BMMSC1). A total of 41 nodes and 120 edges were constructed in the PPI network ( Figure S4). This included a ribosomal protein (RP) family network that played an important role in protein biosynthesis ( Figure S4B). Two sub-networks were identi ed by using MCODE (Figure 5A-5B). SPP1 showed intensive interactions with 11 proteins including COL1A1 and BGLAP in one of the sub-networks ( Figure 5A). We were interested in another sub-network including frizzled class receptor 1(FZD1), WNT inhibitory factor 1(WIF1) and secreted frizzled related protein 4(SFRP4) ( Figure  5B).
FZD1 was a transmembrane receptor that could bind to Wnt ligand and further mediate Wnt canonical and non-canonical pathways [46]. As a member of FZD protein family, the speci c mechanism of how FZD1 regulated cell differentiation had not yet been elucidated. Some studies showed that transcription factor activating protein 2 (AP2) was a negative regulator of osteoblast differentiation and mineralization, and its inhibitory effect may be mediated through the down-regulation of FZD1 expression [46]. Similarly, E2F transcription factor 1 (E2F1) promoted differentiation and mineralization of osteoblasts by upregulating the expression of FZD1 [47]. Therefore, the above ndings suggested that the expression of FZD1 had a positive effect on the differentiation of osteoblasts. When WIF1 and SFRP family members (sFRP1, 2, 3, 4, 5) competitively bound to soluble Wnt ligands, the binding of the Wnt/FZD complex was inhibited [48,49]. Functional studies showed that sFRP4 suppressed the osteoblastic differentiation of human BMSCs [50], and WIF1 played as a negative regulator of osteoblastic differentiation in mouse mesenchymal cells [51]. In our study, SFRP4 was highly expressed in MatOB, WIF1 was highly expressed in ImmatOB and MatOB, and FZD1 was only expressed in BMMSC1 ( Figure 5C). We speculated that mature osteoblasts inhibited the binding of FZD1 and Wnt ligand in BMMSC through expressing sFRP4 and WIF1, hence the inhibited Wnt pathway weakened the osteogenic differentiation ability of BMMSCs ( Figure 5D).

Discussion
BMMSCs are the major source of osteoblasts and adipocytes in bone [52]. Various factors participate in the process of osteogenesis via complex biochemical reactions and cellular communication. Our study integrated scRNA-seq transcriptomes of BMMSCs and osteoblasts in humans, and provided important insights into the differentiation fate of BMMSCs and osteoblasts. In addition, we also analyzed the interaction between BMMSCs and osteoblasts in terms of the ligand receptor pairs.
Previous studies have shown that during osteogenic development at least three stages were recognized: preosteoblasts, immature osteoblasts and mature osteoblasts [31,53]. Some in vivo ndings combined with the previous in vitro data indicated that RUNX2 triggered the expression of downstream genes at an early stage of osteoblast differentiation, but maintained osteoblasts in an immature stage. RUNX2 inhibited osteoblast maturation and the transition to osteocytes, when the expression of RUNX2 began to decrease, the immature osteoblasts progressed to a mature state [54]. The obvious uctuation of RUNX2 expression in our results con rmed the role of RUNX2 in determining the development stage of the osteogenic subgroups. In our previous study, the cluster B6 ( Figure S2A) from BMMSCs were de ned as chondrocytes [13]. After integration analysis, we found that B6 appeared in immature osteoblasts cluster ( Figure 1A). Cluster B6 highly expressed NCAM1, OMD and WIF1 genes as chondrocytes do. However, cluster B6 also highly expressed late-stage markers for osteoblast differentiation, such as IBSP and ENPP1 ( Figure S2B). In contrast, we considered that cluster B6 was closer to late-stage osteoblasts than chondrocytes.
We described two subclusters of osteoblasts (PreOB2 and PreOB3), both of which expressed the marker genes of preosteoblast, such as LEPR, THY1, MGP and IGFBP4. Differently, the cluster PreOB2 highly expressed NR4A1 and NR4A2, which were mainly involved in the regulation of osteoblast differentiation and inhibiting the recruitment of osteoclasts. The cluster PreOB3, which was never reported before, highly expressed ATF3, CCL2, CXCL2, IRF1 and IER3, and had the function of inhibiting cell differentiation, maintaining cells less differentiated status and recruiting osteoclasts. We speculated that high expression of NR4A1 and NR4A2 in preosteoblasts is essential to osteoblast maturation. On the contrary, when preosteoblasts highly expressed ATF3, IRF1 and CXCL2, preosteoblasts probably maintained a state of arrested differentiation and played the function of recruiting osteoclasts. It was also possible that PreOB3 was bone lining cells since PreOB3 was in a less active state. However, the speci c genes expressed in bone lining cells were currently unclear [55,56]. Most studies indicated that the expression of genes in bone lining cells is similar to that in preosteoblasts. Further morphological studies are needed to verify if PreOB3 was bone lining cells.
Each cell in a multicellular organism was not independent, but connected and communicated with other cells in a variety of ways. BMMSCs and osteoblasts play a crucial role in bone metabolism, and they communicate with others through some of cell signal molecules. Under the condition of co-culture with osteoblasts, it was found that the expression levels of ALP and BSP in BMMSCs were higher than that in BMMSCs mono-cultured. And the promotion of ALP and BSP transcription was mediated by gap junctions in direct cell contact [57]. By adding human primary BMMSCs-derived exosomes to the cultured osteoblasts, it was found that the exosomes were endocytosed into osteoblasts to increase calci ed nodules and enhance the intensity of ALP staining [58]. Another rat experiment showed that miR-395 in BMMSCs-derived exosomes can target STAT1 and inhibit STAT1 levels to promote the proliferation and differentiation of osteoblasts. In our analysis, we investigated the cells interactions between BMMSCs and osteoblasts. We observed strong interaction between mature osteoblasts and BMMSCs and COL1A1-DDR2, COL1A1-ITGB1 and SPP1-ITGB1 were possible signi cant ligand-receptor pairs involved in the interaction. Furthermore, we suggested that mature osteoblasts might inhibit the osteogenic ability of BMMSCs through expressing SFRP4 and WIF1 to inhibit the competitive binding of FZD1 to Wnt ligands. Therefore, mature osteoblasts may mediate negative feedback regulation of osteogenesis process.
There were two limitations in the present study. At rst, all cells were collected from diseased individuals. Comparing with healthy individuals, the types or proportions of osteoblast subclusters we got from the data might have some deviations. On the other hand, age was the main factor affecting the transcription pattern of BMMSCs [13]. In this study, the subjects of BMMSCs consisted of a 67-year-old postmenopausal female with osteoporosis and an 85-year-old male with osteoarthritis. The osteoblasts were collected from the femur head of one 31-year-old male patient with osteoarthritis and osteopenia. Although there were certain limitations, we attempted to reveal a continuum of the dynamic gene expression pattern of BMMSCs to osteoblasts, and our dataset and analysis will serve as a resource for future studies investigating osteogenesis differentiation.
In summary, we have sorted out the development trajectory from BMMSCs to mature osteoblasts in vivo and delineated key transcription factors. At the single-cell level, we veri ed the known BMMSCs and osteoblasts subtypes (BMMSC1, BMMSC2, PreOB1, ImmatOB and MatOB), and identi ed a novel osteoblast subtype PreOB3 which highly expressed ATF3, CCL2, CXCL2 and IRF1. We suggested that WIF1 and sFRP4 expressed in osteoblasts inhibited the binding of FZD1 to Wnt ligand in BMMSCs, thereby further inhibiting osteogenic differentiation.

Study population and single cell sequencing
The study subject of osteoblast included one 31-year-old male with osteoarthritis and osteopenia (BMD Tscore: 0.6 at lumbar spine, -1.1 at the total hip). And study subjects of bone marrow mesenchymal stem cells consisted of two subjects: one 67-year-old female with osteoporosis (T-score: -3.3 at lumbar spine, -3.7 at the total hip), the other 84-year-old male with osteoarthritis (T-score: -0.9 at lumbar spine, 2.7 at the total hip).
All subjects were screened with a detailed questionnaire, medical history, physical examination, and bone mineral density (BMD) testing before surgery. Subjects were excluded from this study if they had conditions that affect bone metabolism, including diseases of the kidney, liver, parathyroid, and thyroid, or any of the following conditions: diabetes mellitus, hyperprolactinemia, oophorectomy, ankylosing spondylitis, malabsorption syndromes, malignant tumors, hematologic diseases, or previous pathologic fractures.
The RNA library was constructed using Single Cell 3'Library Gel Bead Kit V3 on the 10x platform. The sequencing principle was based on the droplet method. The single cell suspension was added to the 10x Genomics Chromium Controller for single cell capture and sequencing tag addition. After PCR ampli cation of the library, the two types of single-cell RNA-seq libraries were sequenced separately on the Illumina Novaseq6000 platform.
After the rst pre-screening, we got 3756 BMMSCs. For details, please refer to our previous research [13] and GSE147287. Among them, cluster B5 highly expressed hematopoietic marker gene PTPRC (CD45) and neutrophil marker gene S100A8 ( Figure S2C), so the cluster B5 was eliminated in the subsequent analysis. The B6 cluster was separated from BMMSCs in the dimensionality reduction gure ( Figure S2A). Based on the highly expressed RUNX2, IBSP, and ENPP1 marker genes of cluster B6 ( Figure  S2C), cluster B6 might be in the late stage of osteoblast development, so we kept cluster B6 for subsequent analysis.

Integrated analysis
To analyze the datasets from different cell types, we also used Seurat V3 for integrated analysis. To unify quality control criteria, we removed cells with a high fraction of mitochondrial transcripts (20% for BMMSCs and osteoblasts, respectively), a low unique molecular identi er counts >40000, and a low

Cell cycle scores
To assess cell cycle effects on the dataset, we used the "CellCycleScoring" function in the Seurat package to score each cell. Each cell was divided into G1, G2M, and S phases based on Seurat's built-in cell cycle marker gene dataset (https://satijalab.org/seurat/v3.2/cell_cycle_vignette.html).

Diffusion pseudotime analysis
To delineate the developmental progression of bone marrow mesenchymal cells to mature osteoblasts in pseudotime, we performed the diffusion map analysis [14] as implemented in the R package destiny (bioconductor.org/packages/destiny). The method inferred the low-dimensional manifolds by estimating the eigenvalues and eigenvectors for the diffusion operator related to the data. We used the variable genes identi ed by Seurat as the basis to construct a single cell differentiation trajectory. Then, the diffusion map was constructed using the "DiffusionMap" function of R package destiny.

Pathway enrichment analysis
Gene Ontology (GO) term enrichment analysis was performed using the standalone clusterPro ler package (https://github.com/YuLab-SMU/clusterPro ler). When the adjusted P value was less than 0.05, the term we obtained was considered to be signi cantly enriched.

Cell to cell ligand-receptor interaction analysis
To infer the hypothetical intercellular communication, we used iTALK [34] R packages (https://github.com/Coolgenome/iTALK) to visualize the ligand-receptor networks ( Figure 4A). We used the built-in database of ligand-receptor interaction for iTALK, which collected a total of 2,648 nonredundant and known interacting ligand-receptor pairs. The ligand receptor pairs were classi ed into 4 categories based on the primary function: cytokines/chemokines, immune checkpoint genes, growth factors, and others. iTALK incorporates commonly used algorithms for batch effects corrections and differential gene expression analysis.

Protein-protein interaction network analysis
We constructed the Protein-Protein Interaction (PPI) network based on STRING v10 (https://stringdb.org/) database and Cytoscape [45] software. The rst 30 genes of each subcluster were selected, and then the genes of the two subclusters were input into the STRING v10 database to obtain the basic protein-protein interaction network. We imported the network data into Cytoscape software, and used the built-in program cytoHubba[60] to identify hub genes of the PPI network.
Declarations language, managed editing of the manuscript, supervised writing and approved the nal draft. All authors have read and agreed to the published version of the manuscript.

Con ict of interest
All authors have no con icts of interest to declare.

Availability of data and materials
The scRNA-seq datasets of human BMMSCs and human osteoblasts can be accessed from Gene Expression Omnibus (GEO) database under the accession numbers of GSE147287 and GSE147390, respectively.
Ethics approval and consent to participate  Correlation was tested by Pearson correlation coe cient (R 2 = 0.96, p < 0.01).
(C) Violin plots showed the log-transformed normalized expression levels of the most signi cant marker genes in each cluster.
(D) Dot plot showed selected biological processes (BPs) in the GO analysis of clusters. X-axis, gene ratio; Y-axis, enriched BP terms in clusters; color (red, high; blue, low), -log10(P-adjusted) of each term.