The study subject of osteoblast included one 31-year-old male with osteoarthritis and osteopenia (BMD T-score: 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.

We isolated osteoblasts from the human femoral head according to a previously widely used protocol with only minor modifications (9). Briefly, bone tissue was divided into bone fragments with bone forceps and washed with phosphate-buffered saline (PBS). The fragments were placed into the digestive solution containing a highly purified, endotoxin-free type II collagenase (1 mg/ml, Cat: A004174-0001, Sangon Biotech, Shanghai, China) and digested for two rounds at 37℃, the first round for 30 min and the second round for 60 min. The solution was then collected and filtered through a 70 μm filter and incubated with red blood cell lysis buffer for five minutes. The collected cells were washed twice with PBS and incubated with human CD31/34/45-PE (Cat:303106, Cat:343606, Cat:304008, BioLegend, San Diego, CA, USA), and human ALPL-APC (Cat: FAB1448A, R&D Systems, Minneapolis, MN, USA) antibodies. After incubation, ALP⁺/CD45/34/31⁻ cells were collected as osteoblasts by fluorescence-activated cell sorting (FACS).

BMMSCs were isolated from human bone marrow cells extracted from the bone marrow cavity of the femoral diaphysis according to the widely used isolation protocol (10). Briefly, the bone marrow was diluted with twice the volume of PBS and mixed gently. The mixture was then evenly layered onto an equivalent volume of Ficoll (GE Healthcare, Chicago, Illinois, USA) and the buffy coat was separated by centrifugation (440 g, 35 min, 4℃). The separated buffy coat was transferred to a new 15 ml centrifuge tube and washed with PBS. Then, red blood cells were lysed with RBC Lysis Buffer (Thermo Fisher, Waltham, MA, USA). After washing twice with PBS, the remaining cells were further purified with CD271 magnetic MicroBeads for positive selection (11). BMMSCs were incubated for 10 min at 4℃∼8℃ with monoclonal antibody (mAb) against CD271. After washing, the cells were incubated with anti-IgG1 immunomagnetic beads for 15 min at 4℃. The cell suspension was placed on a column in a cell separator (Miltenyi Biotec), and the positive fraction was subjected to a second separation step. The cells were then counted and assessed for viability with a Countstar^Ⓡ Rigel S3 fluorescence cell analyzer (ALIT Life Science Co., Ltd, Shanghai, China).

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 amplification of the library, the two types of single-cell RNA-seq libraries were sequenced separately on the Illumina Novaseq6000 platform.

For separate preprocessing of BMMSCs and osteoblast data sets, Seurat (12) package V3 (https://satijalab.org/seurat/) was used for filtering, variable gene selection, dimensionality reduction analysis, and clustering. From the dataset of osteoblast, we obtained 7,506 ALPL/RUNX2/COL1A1^high/CD45^low cells. All osteoblasts were clustered by dimensionality reduction, and 9 clusters (O1-O9) were obtained. Among them, clusters O5 and O6 were nucleated erythrocyte groups that highly expressed hemoglobin coding genes (HBA1 and HBB); clusters O8 was smooth muscle cell group that expressed smooth muscle alpha (alpha)-2 actin coding gene (ACTA2) and smooth muscle myosin heavy chain (MYH11) (13); cluster O9 was neutrophil group that expressed the neutrophil related genes S100A8 and hematopoietic marker gene PTPRC (also known as CD45) (Supplementary Fig. S1C and S1D).

After the first pre-screening, we got 3,756 BMMSCs. For details, please refer to our previous research (7) and GSE147287. Among them, cluster B5 highly expressed hematopoietic marker gene PTPRC (CD45) and neutrophil marker gene S100A8 (Supplementary Fig. S2B), so the cluster B5 was eliminated in the subsequent analysis. The B6 cluster was separated from BMMSCs in the dimensionality reduction figure (Supplementary Fig. S2A). Based on the highly expressed IBSP and ENPP1 marker genes of cluster B6 (Supplementary Fig. S2B), cluster B6 might be in the late stage of osteoblast development, so we kept cluster B6 for subsequent 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 identifier counts >40,000, and a low number of detected genes (<200&>5,000) for analysis. Finally, we retained 8,841 cells (including 2,865 BMMSCs and 5,976 osteoblasts) for further analyses. First, the first 2,000 genes with the most obvious cell-to-cell differences were selected by the “FindVariableFeatures” function in the Seurat package. Subsequently, principal component analysis (PCA) analysis was performed on the submatrix consisting of the highly variable genes. The top 20 principal components were used for clustering analysis and visualization.

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).

To delineate the developmental progression of bone marrow mesenchymal cells to mature osteoblasts in pseudotime, we performed the diffusion map analysis (14) as implemen-ted 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 identified 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.

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

To infer the hypothetical intercellular communication, we used iTALK (15) R packages (https://github.com/Coolgenome/iTALK) to visualize the ligand-receptor networks. 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 classified 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.

We constructed the Protein-Protein Interaction (PPI) network based on STRING v10 (https://string-db.org/) database and Cytoscape (16) software. The first 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 (17) to identify hub genes of the PPI network.

The scRNA-seq data of human osteoblasts (8) and BMMSCs (7) 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 (12) and generated a two-dimensional (2D) embedding plot for visualization using the uniform manifold approximation and projection (UMAP) algorithm (18). 8 clusters (C1-C8) were identified based on highly variable genes (Fig. 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. The gene expression profiles between the two subjects were strongly correlated (R²=0.96, Fig. 1B), so we combined the data from the two subjects for subsequent analysis.

Next, we identified highly expressed genes in each cluster (Fig. 1C) and used Gene Ontology (GO) enrichment analysis to characterize the most significant gene expression signatures from each cell cluster (Fig. 1D). Most of clusters were enriched for GO terms related to “ossification”, “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) (Fig. 1C). However, we found cluster BMMSC2 was enriched for terms related to “lipid localization”, “lipid transport”, and “lipid homeostasis” (Fig. 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 (Fig. 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 “ossification”, “osteoblast differentiation” and other terms related to bone development (Fig. 1D), and did not express essential osteogenic genes, such as alkaline phosphatase (ALPL) and integrin subunit β1 (ITGB1) for osteoblast development (Fig. 1C). Hence, as reported in our previous study, cluster BMMSC3 was at the terminal stage of BMMSCs and lacked the capacity for cellular differentiation (7). Therefore, BMMSC1 was a BMMSCs subcluster with osteogenic potential and function. Cells from clusters BMMSC1, PreOB1, PreOB2, PreOB3, ImmatOB and MatOB were classified as osteogenesis cells and used in further analysis.

To assess if clustering was critically influenced 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 (12). We found different cell cycle stages were present at a roughly similar proportion in each subpopulation (Supplementary Fig. S3B and S3C), therefore, the cluster results were not affected by cell cycle heterogeneity.

Next, we inferred osteogenesis differentiation trajectory using unsupervised ordering of the scRNA-seq for osteogenesis cells by diffusion mapping (Fig. 2A) (14). By comparing the distribution of the pseudotime, we found that cluster BMMSC1was enriched in the early stage of pseudotime (Fig. 2B). The clusters ImmatOB and MatOB, which highly expressed specific genes of late osteoblasts stage like IBSP, BGLAP and SPP1 (Fig. 1C and 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 specific genes such as SPARC, COL1A1, BGLAP and SPP1 were upregulated from BMMSCs to mature osteoblasts (Fig. 2C). A subset of genes was consistently downregulated during the osteogenesis differentiation process, such as CXCL12 and LEPR (Fig. 2C). To support the trajectory inference, the expression of significant 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 (Fig. 2D). This result was consistent with the findings of transcriptional continuum in osteoblasts differentiation (19).

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 (Fig. 3A). Comparing the BMMSCs subclusters related to osteogenesis (BMMSC1) with adipogenesis (BMMSC2), we found that the expression of ADAM28 (ADAM metallopeptidase domain 28) was significantly upregulated in BMMSC1, but ADAM28 was almost absent in BMMSC2 (Fig. 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 (20). Overexpression of ADAM28 favored to the proliferation and lineage-specific differentiation of human dental papilla mesenchymal cells (hDPMCs) (21). 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 specific 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) (Fig. 3B). MGP promoted the bone formation by activating the Wnt/β-catenin pathway (22). IGFBP4 reduced the proliferation of BMMSCs, thus inhibiting excessive osteogenic differentiation in rats (23, 24).

Interestingly, we noticed that two preosteoblasts subclusters with opposite functions were novelly identified in our studies. PreOB2 highly expressed nuclear receptor subfamily 4 group A member 1 and 2 (NR4A1 and NR4A2) and PreOB3 highly expressed activating transcription factor 3 (ATF3), CC motif chemokine ligand 2 (CCL2), CXC motif chemokine 2 (CXCL2) and interferon regulatory factor 1 (IRF1) (Fig. 3B). NR4A1 and NR4A2 in cluster PreOB2 were critical negative regulators of osteoclastogenesis and bone resorption, suggesting their bone protective role (25, 26). In PreOB3, IRF1 (27) and ATF3 might act as a negative regulator of osteoblast differentiation. The overexpression of ATF3 inhibited transcription of pro-osteogenic differentiation-associated genes, such as bone morphogenetic protein 2 (BMP2) and ALPL (28). CCL2 played role in recruitment of monocytes and the formation of osteoclasts (29). Studies revealed that CXCL2 promoted osteoclastogenesis and prevented preosteoblast differentiation through inhibiting the ERK1/2 signaling pathway (30, 31). Interestingly, we found that immediate early response 3 (IER3) has the highest expression in the PreOB3 cluster (Fig. 3B). It has been reported that IER3 may be associated with multiple signaling pathways, such as nuclear factor kappa B (NF-κB) pathway and extracellular regulated protein kinases (ERK) pathway (32). We speculated that it may play a role in osteogenic differentiation.

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 (Fig. 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 (33) was highly expressed (Fig. 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 (34). Our results showed that CADM1 was also highly expressed (Fig. 3B). Therefore, CADM1 might play a similar role in cellular communication in human mature osteoblasts. Previous studies showed that lymphocyte-specific protein 1 (LSP1) was mainly expressed in white blood cells (35, 36). A recent study suggested endothelium also expressed LSP1 which played an essential role in transendothelial migration (36). We found that LSP1 was highly enriched in immature osteoblasts and mature osteoblasts (Fig. 3B).

To investigate the biological interaction between BMMSCs and osteoblasts, we performed intercellular communication analysis using an R Package “iTALK” (15), which can identify ligands and receptors in scRNA-seq data (Supplementary Table S1) and analyze mediated cross-signal of cell communication (Fig. 4A). The top 30 pairs of ligand-receptor related to intercellular communication in the osteogenic lineage were visualized (Fig. 4B). We observed strong interaction between mature osteoblasts and BMMSCs. The results showed integrin subunit β1 (ITGB1) receptors were widely present in various subclusters (Fig. 4C). However, the prerequisite for integrin as a receptor was that the integrin β chain (37) must be combined with the integrin α chain to form heterodimeric transmembrane proteins (38) 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 (39). COL1A1 regulated the process of osteogenic differentiation which was mainly mediated by integrin α1/β1 (38, 40, 41). The α chain integrin subunit alpha V (ITGAV) and integrin subunit alpha 1 (ITGA1) were mainly expressed in BMMSC1, PreOB1, PreOB2 and PreOB3 (Fig. 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 significant ligand-receptor pair. COL1A1 is a common marker of osteoblasts, and the expression of COL1A1 is up-regulated during the differentiation of osteoblasts from preosteoblasts to mature osteoblasts (42). Previous studies have shown that DDR2 is a potential cellular marker for osteogenic progenitors and principally binds to collagens type I, II, III and X (43, 44). 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 (44, 45). In our study, DDR2 was mostly expressed in BMMSC1, and COL1A1 was highly expressed in MatOB (Fig. 4C). Therefore, COL1A1 specifically expressed by mature osteoblasts could mediate the differentiation of BMMSCs towards osteogenic direction.

In order to identify significant 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 (16) 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 (Supplementary Fig. S4A and S4B). This included a ribosomal protein (RP) family network that played an important role in protein biosynthesis (Supplementary Fig. S4B). Two sub-networks were identified by using MCODE (Fig. 5A and 5B). SPP1 showed intensive interactions with 11 proteins including COL1A1 and BGLAP in one of the sub-networks (Fig. 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) (Fig. 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 specific 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 findings 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). SFRP4 expression has been reported in the osteogenic lineage, and previous studies using in vitro assays have shown that SFRP4 inhibits osteoblast proliferation and differentiation, at least in part, by attenuating classical Wnt signaling in vivo (50, 51). Functional studies showed that sFRP4 suppressed the osteoblastic differentiation of human BMSCs (52). WIF1 expression has been detected in osteoblasts in previous studies (53, 54), and WIF1, as a secretory regulator of Wnt signaling, inhibited Wnt signaling and negatively regulated osteogenic differentiation of mouse mesenchymal cells in vitro (55). In our study, SFRP4 was highly expressed in MatOB, WIF1 was highly expressed in ImmatOB and MatOB, and FZD1 was only expressed in BMMSC1 (Fig. 5C). We verified their expression in the latest sample data collected (Supplementary Fig. S4C), in which FZD1 was highly expressed in BMMSC, WIF1 and SFRP4 were highly expressed in MatOB. 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 (Fig. 5D).