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

doi:10.21203/rs.3.rs-1209239/v1

PDF

Research Article

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

https://doi.org/10.21203/rs.3.rs-1209239/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 04 Jan, 2022

You are reading this latest preprint version

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 (scRNA-seq) transcriptomes data from BMMSCs and osteoblasts. We successfully identified 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.

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.

Single-cell RNA sequencing (scRNA-seq)

Bone marrow mesenchymal stem cells (BMMSCs)

Osteoblasts

Cellular heterogeneity

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–5]. Osteoblastogenic development is a complex process that requires the tight regulation of specific 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 co-receptor 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 identified. 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 inflammation[8]. Our previous studies performed single-cell sequencing of human BMMSCs and osteoblasts. In BMMSCs population, we identified osteoblast precursors, adipocyte precursors, chondrocyte precursors and terminal-stage cells subpopulations in LEPR^highCD45^low BMMSCs[9]. In osteoblasts, we identified 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.

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 two-dimensional (2D) embedding plot for visualization using the uniform manifold approximation and projection (UMAP) algorithm[12]. 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 (R² = 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[13]. 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[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 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[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 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[16]. Overexpression of ADAM28 favored to the proliferation and lineage-specific 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 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[18]. 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[23] 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[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-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[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 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[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 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[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]. 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).

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 findings 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 fluctuation of RUNX2 expression in our results confirmed 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 defined 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 specific 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 calcified 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 significant 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 first, 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 verified the known BMMSCs and osteoblasts subtypes (BMMSC1, BMMSC2, PreOB1, ImmatOB and MatOB), and identified 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 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.

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.

Preprocessing of the scRNA-seq data

For separate preprocessing of BMMSCs and osteoblast data sets, Seurat[11] 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^lowcells. 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)[59]; cluster O9 was neutrophil group that expressed the neutrophil related genes S100A8 and hematopoietic marker gene PTPRC(also known as CD45) (Figure S1C-S1D).

After the first 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 figure (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 identifier counts >40000, and a low number of detected genes (<200&>5000) for analysis. Finally, we retained 8,841 cells (including 2,865 BMMSCs and 5,976 osteoblasts) for further analyses. First, the first 2000 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.

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

Pathway enrichment analysis

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.

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

Protein-protein interaction network analysis

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

Acknowledgement

This research was benefited by grants from the National Institutes of Health (R01AR069055, U19AG055373, P20GM109036, R01AG061917), National Natural Science Foundation of China(NSFC) (Grant No. 81902277and 81570807), National Key R&D Program of China (Grant No. 2017YFC1001100), Natural Science Foundation of Hunan Province (S2019JJQNJJ2093), Changsha Science and technology project (kq1907153), Central South University (Grant Nos. 164990007, 2018zzts886), and Xiangya Clinical Big Data Project (xyyydsj9).

Author contributions

HWD and LJT conceived the study and designed the experiments. YG and ZW provided technical support. XHL, HXZ, CZ, YC, YL and XQ performed most of the experiments. YL analyzed the data, wrote the manuscript and designed the figures. YC revised the manuscript. LJT, HS, HMX and HWD corrected the language, managed editing of the manuscript, supervised writing and approved the final draft. All authors have read and agreed to the published version of the manuscript.

Conflict of interest

All authors have no conflicts 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

The study was approved by the Medical Ethics Committee of Xiangya Hospital of Central South University, and the IRB approval number is No. 201912315.

Martin TJ, Seeman E. Bone remodelling: its local regulation and the emergence of bone fragility. Best Pract Res Clin Endocrinol Metab. 2008;22:701–22. doi:10.1016/j.beem.2008.07.006.
Hadjidakis DJ, Androulakis II. Bone remodeling. Ann N Y Acad Sci. 2006;1092:385–96.
Muruganandan S, Roman AA, Sinal CJ. Adipocyte differentiation of bone marrow-derived mesenchymal stem cells: Cross talk with the osteoblastogenic program. Cell Mol Life Sci. 2009;66:236–53.
Muruganandan S, Sinal CJ. The impact of bone marrow adipocytes on osteoblast and osteoclast differentiation. IUBMB Life. 2014;66:147–55.
Chen Q, Shou P, Zheng C, Jiang M, Cao G, Yang Q, et al. Fate decision of mesenchymal stem cells: Adipocytes or osteoblasts? Cell Death Differ. 2016;23:1128–39.
Zieba JT, Chen YT, Lee BH, Bae Y. Notch signaling in skeletal development, homeostasis and pathogenesis. Biomolecules. 2020;10:1–18.
Maeda K, Kobayashi Y, Koide M, Uehara S, Okamoto M, Ishihara A, et al. The regulation of bone metabolism and disordersby wnt signaling. Int J Mol Sci. 2019;20.
Liu S, Stroncek DF, Zhao Y, Chen V, Shi R, Chen J, et al. Single cell sequencing reveals gene expression signatures associated with bone marrow stromal cell subpopulations and time in culture. J Transl Med. 2019;17:1–13. doi:10.1186/s12967-018-1766-2.
Wang Z, Li X, Yang J, Gong Y, Zhang H, Qiu X, et al. Single-cell RNA sequencing deconvolutes the in vivo heterogeneity of human bone marrow-derived mesenchymal stem cells. Int J Biol Sci. 2021;17:4192–206.
Gong Y, Yang J, Li X, Zhou C, Chen Y, Wang Z, et al. A systematic dissection of human primary osteoblasts in vivo at single-cell resolution. Aging (Albany NY). 2021;13:20629–50.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36:411–20.
Becht E, McInnes L, Healy J, Dutertre C-A, Kwok IWH, Ng LG, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2018.
Wang Z, Li X, Yang J, Gong Y, Zhang H, Qiu X, et al. Single-cell RNA sequencing deconvolutes the in vivo heterogeneity of human bone marrow-derived mesenchymal stem cells. bioRxiv. 2020;:
2020.04.06.027904
. doi:10.1101/2020.04.06.027904.
Angerer P, Haghverdi L, Büttner M, Theis FJ, Marr C, Buettner F. destiny: diffusion maps for large-scale single-cell data in R. Bioinformatics. 2016;32:1241–3.
Rutkovskiy A, Stensløkken K-O, Vaage IJ. Osteoblast Differentiation at a Glance. Med Sci Monit Basic Res. 2016;22:95–106.
Zhao Z, Wen LY, Jin M, Deng ZH, Jin Y. ADAM28 participates in the regulation of tooth development. Arch Oral Biol. 2006;51:996–1005.
Zhao Z, Tang L, Deng Z, Wen L, Jin Y. Essential role of ADAM28 in regulating the proliferation and differentiation of human dental papilla mesenchymal cells (hDPMCs). Histochem Cell Biol. 2008;130:1015–25.
Zhang J, Ma Z, Yan K, Wang Y, Yang Y, Wu X. Matrix Gla Protein Promotes the Bone Formation by Up-Regulating Wnt/β-Catenin Signaling Pathway. Front Endocrinol (Lausanne). 2019;10 December:1–11.
Li H, Yu S, Hao F, Sun X, Zhao J, Xu Q, et al. Insulin-like growth factor binding protein 4 inhibits proliferation of bone marrow mesenchymal stem cells and enhances growth of neurospheres derived from the stem cells. Cell Biochem Funct. 2018;36:331–41.
Wu J, Wang C, Miao X, Wu Y, Yuan J, Ding M, et al. Age-Related Insulin-Like Growth Factor Binding Protein-4 Overexpression Inhibits Osteogenic Differentiation of Rat Mesenchymal Stem Cells. Cell Physiol Biochem. 2017;42:640–50.
Martínez-González J, Badimon L. The NR4A subfamily of nuclear receptors: new early genes regulated by growth factors in vascular cells. Cardiovasc Res. 2005;65:609–18.
Li X, Wei W, Huynh H, Zuo H, Wang X, Wan Y. Nur77 prevents excessive osteoclastogenesis by inducing ubiquitin ligase Cbl-b to mediate NFATc1 self-limitation. Elife. 2015;4 JULY 2015:1–17.
Salem S, Gao C, Li A, Wang H, Nguyen-Yamamoto L, Goltzman D, et al. A novel role for interferon regulatory factor 1 (IRF1) in regulation of bone metabolism. J Cell Mol Med. 2014;18:1588–98.
Han Y, Zhang K, Hong Y, Wang J, Liu Q, Zhang Z, et al. miR-342-3p promotes osteogenic differentiation via targeting ATF3. FEBS Lett. 2018;592:4051–65.
Park J kyung, Jang H, Hwang SS, Kim EJ, Kim DE, Oh KB, et al. ER stress-inducible ATF3 suppresses BMP2-induced ALP expression and activation in MC3T3-E1 cells. Biochem Biophys Res Commun. 2014;443:333–8. doi:10.1016/j.bbrc.2013.11.121.
Mulholland BS, Forwood MR, Morrison NA. Monocyte Chemoattractant Protein-1 (MCP-1/CCL2) Drives Activation of Bone Remodelling and Skeletal Metastasis. Curr Osteoporos Rep. 2019;17:538–47.
Siddiqui JA, Partridge NC. CCL2/monocyte chemoattractant protein 1 and parathyroid hormone action on bone. Front Endocrinol (Lausanne). 2017;8 MAR:1–7.
Yang Y, Zhou X, Li Y, Chen A, Liang W, Liang G, et al. CXCL2 attenuates osteoblast differentiation by inhibiting the ERK1/2 signaling pathway. J Cell Sci. 2019;132.
Ha J, Lee Y, Kim HH. CXCL2 mediates lipopolysaccharide-induced osteoclastogenesis in RANKL-primed precursors. Cytokine. 2011;55:48–55. doi:10.1016/j.cyto.2011.03.026.
Semler O, Garbes L, Keupp K, Swan D, Zimmermann K, Becker J, et al. A mutation in the 5′-UTR of IFITM5 creates an in-frame start codon and causes autosomal-dominant osteogenesis imperfecta type v with hyperplastic callus. Am J Hum Genet. 2012;91:349–57.
Inoue T, Hagiyama M, Enoki E, Sakurai MA, Tan A, Wakayama T, et al. Cell adhesion molecule 1 is a new osteoblastic cell adhesion molecule and a diagnostic marker for osteosarcoma. Life Sci. 2013;92:91–9. doi:10.1016/j.lfs.2012.10.021.
Pulford K, Jones M, Banham AH, Haralambieva E, Mason DY. Lymphocyte-specific protein 1: A specific marker of human leucocytes. Immunology. 1999;96:262–71.
Liu L, Cara DC, Kaur J, Raharjo E, Mullaly SC, Jongstra-Bilen J, et al. LSP1 is an endothelial gatekeeper of leukocyte transendothelial migration. J Exp Med. 2005;201:409–18.
Wang Y, Wang R, Zhang S, Song S, Jiang C, Han G, et al. iTALK: an R Package to Characterize and Illustrate Intercellular Communication. bioRxiv. 2019;:507871.
Plow EF, Haas TA, Zhang L, Loftus J, Smith JW. Ligand binding to integrins. J Biol Chem. 2000;275:21785–8.
Campbell ID, Humphries MJ. Integrin structure, activation, and interactions. Cold Spring Harb Perspect Biol. 2011;3:1–14.
Marie PJ. Targeting integrins to promote bone formation and repair. Nat Rev Endocrinol. 2013;9:288–95. doi:10.1038/nrendo.2013.4.
Chen Q, Shou P, Zhang L, Xu C, Zheng C, Han Y, et al. An osteopontin-integrin interaction plays a critical role in directing adipogenesis and osteogenesis by mesenchymal stem cells. Stem Cells. 2014;32:327–37.
Reyes CD, García AJ. Alpha2beta1 integrin-specific collagen-mimetic surfaces supporting osteoblastic differentiation. J Biomed Mater Res A. 2004;69:591–600.
Gronthos S, Simmons PJ, Graves SE, G. Robey P. Integrin-mediated interactions between human bone marrow stromal precursor cells and the extracellular matrix. Bone. 2001;28:174–81.
Fu HL, Valiathan RR, Arkwright R, Sohail A, Mihai C, Kumarasiri M, et al. Discoidin domain receptors: Unique receptor tyrosine kinases in collagen-mediated signaling. J Biol Chem. 2013;288:7430–7.
Leitinger B, Kwan APL. The discoidin domain receptor DDR2 is a receptor for type X collagen. Matrix Biol. 2006;25:355–64.
Ruiz PA, Jarai G. Collagen I induces discoidin domain receptor (DDR) 1 expression through DDR2 and a JAK2-ERK1/2-mediated mechanism in primary human lung fibroblasts. J Biol Chem. 2011;286:12912–23.
Yang H, Sun L, Cai W, Gu J, Xu D, Deb A, et al. DDR2, a discoidin domain receptor, is a marker of periosteal osteoblast and osteoblast progenitors. J Bone Miner Metab. 2020;38:670–7. doi:10.1007/s00774-020-01108-y.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Yu S, Yerges-Armstrong LM, Chu Y, Zmuda JM, Zhang Y. AP2 suppresses osteoblast differentiation and mineralization through down-regulation of Frizzled-1. Biochem J. 2015;465:395–404.
Yu S, Yerges-Armstrong LM, Chu Y, Zmuda JM, Zhang Y. E2F1 effects on osteoblast differentiation and mineralization are mediated through up-regulation of frizzled-1. Bone. 2013;56:234–41. doi:10.1016/j.bone.2013.06.019.
Katoh Y, Katoh M. Comparative genomics on Dkk1 orthologs. Int J Oncol. 2005;27:275–9.
Pez F, Lopez A, Kim M, Wands JR, Caron De Fromentel C, Merle P. Wnt signaling and hepatocarcinogenesis: Molecular targets for the development of innovative anticancer drugs. J Hepatol. 2013;59:1107–17. doi:10.1016/j.jhep.2013.07.001.
Yamada A, Iwata T, Yamato M, Okano T, Izumi Y. Diverse functions of secreted frizzled-related proteins in the osteoblastogenesis of human multipotent mesenchymal stromal cells. Biomaterials. 2013;34:3270–8. doi:10.1016/j.biomaterials.2013.01.066.
Cho SW, Yang JY, Sun HJ, Jung JY, Her SJ, Cho HY, et al. Wnt inhibitory factor (WIF)-1 inhibits osteoblastic differentiation in mouse embryonic mesenchymal cells. Bone. 2009;44:1069–77. doi:10.1016/j.bone.2009.02.012.
Yue R, Zhou BO, Shimada IS, Zhao Z, Morrison SJ. Leptin Receptor Promotes Adipogenesis and Reduces Osteogenesis by Regulating Mesenchymal Stromal Cells in Adult Bone Marrow. Cell Stem Cell. 2016;18:782–96.
Aubin JE. Regulation of Osteoblast Formation and Function. Reviews in Endocrine and Metabolic Disorders. 2001;2:81–94.
Komori T. Regulation of osteoblast differentiation by transcription factors. J Cell Biochem. 2006;99:1233–9.
Insua A, Monje A, Wang H-L, Miron RJ. Basis of bone metabolism around dental implants during osseointegration and peri-implant bone loss. J Biomed Mater Res A. 2017;105:2075–89.
Matic I, Matthews BG, Wang X, Dyment NA, Worthley DL, Rowe DW, et al. Quiescent Bone Lining Cells Are a Major Source of Osteoblasts During Adulthood. Stem Cells. 2016;34:2930–42.
Mikami Y, Yamamoto K, Akiyama Y, Kobayashi M, Watanabe E, Watanabe N, et al. Osteogenic gene transcription is regulated via gap junction-mediated cell-cell communication. Stem Cells Dev. 2015;24:214–27.
Yang X, Yang J, Lei P, Wen T. LncRNA MALAT1 shuttled by bone marrow-derived mesenchymal stem cells-secreted exosomes alleviates osteoporosis through mediating microRNA-34c/SATB2 axis. Aging (Albany NY). 2019;11:8777–91.
Brun J, Lutz KA, Neumayer KMH, Klein G, Seeger T, Uynuk-Ool T, et al. Smooth muscle-like cells generated from human mesenchymal stromal cells display marker gene expression and electrophysiological competence comparable to bladder smooth muscle cells. PLoS One. 2015;10:1–21.
Chin C-H, Chen S-H, Wu H-H, Ho C-W, Ko M-T, Lin C-Y. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8 Suppl 4 Suppl 4:S11.

No competing interests reported.

PDF

Version 1

posted 04 Jan, 2022

You are reading this latest preprint version

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

Status:

Version 1

Abstract

Figures

Background

Results

Discussion

Methods

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1