Integration of single-cell transcriptome sequencing data from osteoblasts and BMMSCs
The scRNA-seq data of human osteoblasts and BMMSCs 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 and generated a two-dimensional (2D) embedding plot for visualization using the uniform manifold approximation and projection (UMAP) algorithm. 8 clusters (C1-C8) were identified 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 profile between the two subjects was strongly correlated (R2 = 0.96, Figure 1B). Hence the batch effect between the two subjects was relatively small and did not affect subsequent analysis.
Next, we identified highly expressed genes in each cluster (Figure 1C) and used Gene Ontology (GO) enrichment analysis to characterize the most significant gene expression signatures from each cell cluster (Figure 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) (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 “ossification”, “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. 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. 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). 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 specific 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 specific 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 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 (Figure 2D). This result was consistent with the findings of transcriptional continuum in osteoblasts differentiation.
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 significantly 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. Overexpression of ADAM28 favored to the proliferation and lineage-specific differentiation of human dental papilla mesenchymal cells (hDPMCs). 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) (Figure 3B). MGP promoted the bone formation by activating the Wnt/β-catenin pathway. IGFBP4 reduced the proliferation of BMMSCs, thus inhibiting excessive osteogenic differentiation in rats[19, 20].
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) (Figure 3B). NR4A1 and NR4A2 in cluster PreOB2 were critical negative regulators of osteoclastogenesis and bone resorption, suggesting their bone protective role[21, 22]. In PreOB3, IRF1 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 [24, 25]. CCL2 played role in recruitment of monocytes and the formation of osteoclasts[26, 27]. 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 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. 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-specific 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. 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", 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 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. 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 significant 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 specifically expressed by mature osteoblasts could mediate the differentiation of BMMSCs towards osteogenic direction.
Protein-protein interaction network
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 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 identified 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. 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. Similarly, E2F transcription factor 1 (E2F1) promoted differentiation and mineralization of osteoblasts by upregulating the expression of FZD1. 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]. Functional studies showed that sFRP4 suppressed the osteoblastic differentiation of human BMSCs, and WIF1 played as a negative regulator of osteoblastic differentiation in mouse mesenchymal cells. 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).