hBMSCs osteogenic, adipogenic, and chondrogenic differentiation 14-day transcriptome data cross-identify 31 DEGsA preliminary comparison of the six datasets (GSE36923, GSE44303, GSE109503, GSE140861, GSE28205, GSE37558) was performed before the differential gene analysis. Through the box-and-line plot of the distribution of the 6 datasets, we concluded that there was no batch effect within the datasets (Fig. 1a–c). According to the PCA plot of the distribution of the 6 datasets, the differences between different types of data were significant, and the analysis was feasible (Fig. 1d–i). Differential gene analysis of different types of data was performed separately, and DEGs coexpressed in similar datasets were selected to obtain 347 upregulated genes and 263 downregulated genes related to osteogenic differentiation, 658 upregulated genes and 704 downregulated genes related to adipogenic differentiation, and 106 upregulated genes and 48 downregulated genes related to cartilaginous differentiation (Fig. 2a–i). There was an antagonistic relationship between the three classifications. To identify the genes that were coexpressed in each category based on the Venn diagram, the genes coexpressed in osteogenic, adipogenic, and chondrogenic differentiation were compared two by two with opposite differentiation trends, and a total of 31 DEGs were obtained (Fig. 2j–o).Figure 1Quality of RNA-seq data of bone marrow-derived MSCs at the late stage of differentiation. Box plot: (a) osteogenic differentiation dataset (GSE28205, GSE37558); (b) adipogenic differentiation dataset (GSE36923, GSE44303); (c) chondrogenic differentiation dataset (GSE109503, GSE140861); PCA plot; (d) GSE37558; (e) GSE36923; (f) GSE140861; (g) GSE28205; (h) GSE44303; (i) GSE109503.Figure 2Target gene identification and enrichment analysis of RNA-seq data from bone marrow-derived MSCs at the late stage of differentiation (a,b) volcano plots of differential gene expression in osteogenic differentiation datasets (GSE28205 and GSE37558); (c,d) volcano plots of differential gene expression in adipogenic differentiation datasets (GSE36923 and GSE44303); (e,f) volcano plots of differential gene expression in cartilageogenic differentiation datasets (GSE109503 and GSE140861) volcano plots of differential gene expression; (g) heatmap of differential gene expression of osteogenic differentiation datasets (GSE28205 and GSE37558) after integration by the RRA algorithm; (h) heatmap of differential gene expression of adipogenic differentiation datasets (GSE36923 and GSE44303) after integration of by the RRA algorithm; (i) heatmap of differential gene expression of chondrogenic differentiation dataset (GSE109503 and GSE140861) heatmap of differential gene expression after integration by the RRA algorithm. Venn diagrams of DEGs with opposite differentiation trends in the three differentiation directions. A total of 31 DEGs were obtained. (j) Intersection of osteogenic upregulated differential genes and adipogenic downregulated differential genes; (k) intersection of osteogenic upregulated differential genes and chondrogenic downregulated differential genes; (l) intersection of chondrogenic upregulated differential genes and adipogenic downregulated differential genes; (m) intersection of adipogenic upregulated differential genes and chondrogenic downregulated differential genes; (n) intersection of adipogenic upregulated differential genes and chondrogenic downregulated differential genes; (o) intersection of adipogenic upregulated differential genes and osteogenic downregulated differential genes; (p) differential gene enrichment analysis (GO) bar graph; (q) differential gene enrichment analysis (KEGG) network graph. AD adipogenesis, OS osteogenesis, CH chondrogenesis.Enrichment analysis revealed that DEGs were concentrated in the extracellular matrix and related to ECM-receptor interactionsGene Ontology (GO) analysis revealed that the differentially expressed genes (DEGs) associated with an adjusted p-value (adj. P) < 0.05 were primarily enriched in biological processes and molecular functions related to collagen-containing extracellular matrix, hexosaminidase activity, binding to laminin, structural components of the extracellular matrix that confer tensile strength, and integrin binding (Fig. 2p). At KEGG pathway analysis revealed that DEGs with adj.P < 0.05 were associated with the ECM-receptor interaction, local adhesion, and glycolipid biosynthesis-ganglio series pathways (Fig. 2q).ceRNA network constructionUsing the String website, six target genes enriched in extracellular mesenchymal tissues were identified, namely, Hapln1, Col4a1, Lamc1, Itga10, Col10a14, and Tln2, and two genes enriched in the target genes involved in the biosynthesis of glycosphingolipid-ganglio series pathways were identified, namely, Hexb and St3gal5 (Fig. 3a,b). The gene symbols, abbreviations and functions are shown in Table 1. The two PPI networks of the target genes (60 nodes, 115 edges; 19 nodes, 20 edges) were obtained through the MCODE module of Cytoscape, with retention greater than or equal to 2, K-Core taken as 2, and Max. depth taken as 100 nodes. mRNAs are pink, miRNAs are purple, and lncRNAs are orange; the greater the connectivity of the nodes is, the greater the number of nodes. Lamc1 and Col4a1 are the core of the mRNA-miRNA-lncRNA interaction network and are closely related to noncoding RNAs that have been shown to play important roles in osteogenesis or adipogenesis, such as MALAT1, H19, XIST, and NEAT1.Figure 3Expression of target genes in the database (a) mRNA–miRNA–lncRNA interaction network map composed of Hapln1, Col4a1, Lamc1, Itga10, Col10a1, and Tln2; (b) mRNA–miRNA–lncRNA interaction network map composed of Hexb and St3gal5. Expression in the HPA database; (c) histogram of differential gene tissue expression levels; (d) histogram of differential gene gene expression levels; expression in the Bgee database; (e) histogram of differential gene tissue expression levels.Table 1 Eight hub genes and their functions.Target gene expression was higher at both the adipose gene level and tissue level than at bone tissueAnalysis of the Human Protein Atlas (HPA) database revealed that the expression levels of Col4a1 and Lamc1 were significantly higher in adipose tissue compared to bone marrow tissue, both at the gene and tissue levels. Similarly, Tln2 exhibited slightly elevated gene expression in adipose tissue relative to bone marrow (Fig. 3c,d). Interestingly, while Hexb and St3gal5 displayed moderately higher gene expression in adipose tissue than in bone marrow, their tissue-level expression patterns were reversed, with higher levels observed in bone marrow compared to adipose tissue. The gene expression scores of the hub genes in abdominal adipose tissue, omental fat pads, subcutaneous adipose tissue, synovial joints, cartilage tissue, trabecular bone tissue, tibia, bone marrow, and myeloid cells were compared using the data obtained from the BGEE database. Col4a1 and Lamc1 had the highest expression in abdominal adipose tissue, omental fat pads, and subcutaneous adipose tissue, and Hexb had the highest expression in synovial joints, trabecular bone tissues, and bone marrow cells. Hapln1 was most highly expressed in cartilage tissues, and only Tln2 and Itga10 were expressed in bone marrow (Fig. 3e).BMSCs are categorized into 11 clusters at single-cell resolutionWe used t-distributed stochastic neighbor embedding (tSNE) to investigate cellular heterogeneity within clusters, and after the integration of the two datasets (GSE128423, GSE156635) (n = 44053), we identified 21 cell clusters and 11 cell types-osteoblasts (n = 2411), chondrocytes (n = 8251), adipogenic cells (n = 5396), endothelial cells (n = 7682), fibroblasts (n= 6949), neutrophils (n = 912), lymphocytes (n = 3957), megakaryocytes (n= 682), erythroblasts (n= 735), stem cells (n = 4908) and blastocytes (n = 2170) -and the biomarker expression of each subgroup is given (Fig. 4a,b). Erythrocytes were not removed since our focus was on the directed differentiation of stem cells to adipocytes, osteoblasts, and chondrocytes. The distribution of the target genes is shown in Fig. 4c.Figure 4Subcluster identification of bone marrow-derived MSC single-cell sequencing data. (a) The TSNE algorithm was applied to the first 17 PCAs for dimensionality reduction, and 22 cell clusters were successfully classified; (b) cell clusters were manually annotated with CellMarker 2.0 according to the composition of marker genes and were successfully annotated into 11 subpopulations; (c) the expression levels of the target genes (Col4a1,Lamc1,Itga10,Hapln1,Hexb,Col10a1,Tln2,St3gal5)in the subpopulations; (d) the marker gene bubble map. Reliability of subcluster annotation of single-cell data.Scoring of scRNA-seq annotation reliability proves cluster annotation reliabilityWe chose three methods to evaluate the reliability of the single-cell annotations. First, we used the WNT pathway and the ADIPOGENESIS pathway to sort the cells (Fig. 5a,b). The WNT pathway was mainly distributed in osteoblasts, adipocytes, stromal cells, and chondroblasts, and the ADIPOGENESIS pathway was mainly distributed in adipocytes, stromal cells, and megakaryocytes, which was in line with our biological a priori experience. Then using the local similarity between cells, the stemness and differentiation potential of cells were assessed according to CytoTrace, and the degree of differentiation was in the following order: stem cells < chondrocytes < osteoblasts < adipoblasts. Col4a1 and Lamc1 were mainly distributed in adipocytes and osteoblasts; Tln2 and Hexb were expressed in chondrocytes, osteoblasts, and adipocytes; Col10a1 was distributed in a small number of chondrocytes, Itga10 and Hapln1 were distributed in a large number of chondrocytes, and St3gal5 was distributed in a large number of stem cells (Fig. 5c). Using Monocle2 to perform a time-matching analysis of the four types of cells, it was concluded that the BMSCs can be divided into two different types of sources, both of which can be differentiated into chondroblasts, and that the later stage is accompanied by the differentiation of the adipoblasts with the osteoblasts with a large number of adipoblasts (Fig. 5d). Finally, cell cycle maps for all the subpopulations were generated and we concluded that stem cells, lymphocytes, and neutrophils were in the early stage of differentiation, while the remaining subpopulations were in the late stage of differentiation (Fig. 5e).Figure 5Proof of subgroup annotation reliability (a) the WNT pathway and the adult ADIPOGENESIS pathway were selected to score the cell subclusters; (b) the differentiation of the four classes of cells was assessed according to Cytotrace; (c) the expression levels of the target genes (Col4a1,Lamc1,Itga10,Hapln1,Hexb,Col10a1,Tln2,St3gal5)according to Cytotrace; (d) the differentiation of the four classes of cells was assessed based on Monocle2; (e)the differentiation of the 11 classes of cells was assessed based on Tricycle.Cell cycle analysis of target genesBased on the scatter plots, Col4a1, Lamc1, Hapln1, Itga10, Tln2, Hexb, were expressed at all times in osteoblasts, Col10a1 and St3gal5 were expressed to a lesser extent, and Itga10 with Lamc1 became more abundantly expressed at mitotic phase (M). Scatter plots of Runx2 and Bglap expression cycles are given as reference (Fig. 6). In adipocytes, Col4a1 and Lamc1 were significantly more expressed than in osteoblasts, being abundantly expressed during mitosis (M), and Itga10, St3gal5, and Col10a1 were expressed during DNA synthesis (S), which is consistent with our previous prediction. Scatter plots of Ppar\(\gamma \) and Fabp4 are given as reference (Fig. 7). In chondrocytes, Hapln1 and Itga10 expression is reduced and Tln2 expression is elevated in mitosis (M). Scatter plots of Sox9 and Col2a1 are given as reference (Fig. 8).Figure 6Expression period of target genes of Lamc1 interaction in osteoblasts (a) Col4a1; (b) Lamc1; (c) Itga10; (d) Hapln1; (e) Hexb; (f) Col10a1; (g) Tln2; (h) St3gal5;(i) Bglap; (j) Runx2.Figure 7Expression period of target genes of Lamc1 interaction in in adipocytes (a) Col4a1; (b) Lamc1; (c) Itga10; (d) Hapln1; (e) Hexb; (f) Col10a1; (g) Tln2; (h) St3gal5;(i) Adipoq; (j) Fabp4.Figure 8Expression period of target genes of Lamc1 interaction in chondrocy (a) Col4a1; (b) Lamc1; (c) Itga10; (d) Hapln1; (e) Hexb; (f) Col10a1; (g) Tln2; (h) St3gal5; (i) Col2a1; (j) Sox9.Intercellular interactions in adipogenic, osteogenic differentiation target two pairs of ligands that interact with Lamc1In BMSCs, osteoblasts communicate with adipogenic cells when adipogenic cells are used as receptors, and all three types of cells communicate with stem cells when stem cells are used as receptors (Fig. 9a). Considering that the pathway enriched for target genes is the ECM-receptor interaction pathway, to determine the mechanism of interactions in the stroma, we investigated the most significant ligand–receptor interactions in the LAMININ pathway (Fig. 9b). Lamc1-CD44 and Lamc1-Dag1 exhibited extreme differentiation, followed by interactions at the level of individual ligands, and the ligands that interact with Lamc1 were identified. LAMC1-CD44 interaction may affect stem cell differentiation to bone or cartilage and adipose, and LAMC1-DAG1 interaction may affect osteogenic, adipogenic, and chondrogenic differentiation relationships (Fig. 9c–e).Figure 9The molecular mechanism of Lamc1 interaction (a) hierarchical plot of cellular communication; (b) the most important ligand-receptor interactions in the laminin pathway (top 20); (c) LAMC1-Dag1 interaction chord diagram; (d) LAMC1-CD44 interaction chord diagram; (e) Demonstrate the specific ligand–receptor interaction mechanisms of stem cells, osteoblasts, and adipogenic and chondrogenic cells in the Laminin pathway.Target genes promote osteogenic differentiation and inhibit adipogenic differentiationThe target genes Col10a1 and Hapln are predicted to play important roles in chondrogenic differentiation, and we have confirmed this in past reports. Similarly, Col4a1 and Itga10 have been shown to play important roles in osteogenic differentiation, and therefore, no subsequent validation will be performed. During osteogenic differentiation, the alizarin red staining results revealed gradual reddening with increasing days of induction, with MC3T3-E1 cells reaching the reddest at 7 days (Fig. 10a,b). We detected the expression of target genes and marker genes at 3 days, 5 days, and 7 days. Lamc1, Tln2, Hexb, and St3gal5 tended to be upregulated, and the expression of the marker gene Runx2 gradually increased to a maximum at 7 days (Fig. 10c–e). During adipogenic differentiation, as the duration of induction increased, the red intensity of the 3T3-L1 cells gradually decreased, and on day 8, the red intensity of the 3T3-L1 cells decreased (Fig. 10f,g). We detected the expression of target genes and marker genes after 4 days, 6 days, and 8 days. Lamc1, Tln2, Hexb, and St3gal5 tended to be downregulated, and the expression of the marker gene Ppar\(\gamma \) gradually increased to a maximum on 8 days (Fig. 10h–j). Laminin can increase the expression of calcium and increase the concentration of osteocalcin. After Lamc1 was knocked down in MC3T3-E1 and 3T3-L1 cells, the gene expression levels of Runx2 and Ocn were significantly downregulated, and those of Ppar\(\gamma \) and Fabp4 were also downregulated (Fig. 11c–e,h–j).Similarly, the red color of the knocked-down cells was significantly lighter than that of the controls (Fig. 11a,b,f,g). Similarly, the previous experiment was repeated with the housekeeping gene GAPDH and yielded similar results (Supplementary Fig. S1).Figure 10Results of target gene expression experiments. (a) Alizarin red staining of osteoblasts induced for 3, 5 and 7 days; (b) capture images with camera,100\(\times \) microscope,scale bar \(=500\;\upmu \hbox {m}\); qPCR results of target gene expression in osteoblasts: (c) 3 days; (d) 5 days; (e) 7 days. (f) Oil Red O staining of adipocytes induced for 4, 6 and 8 days; (g) capture images with camera, 100\(\times \) microscope, scale bar \(=\) \(500\;\upmu \hbox {m}\); qPCR results of target gene expression in adipocyte: (h) 4 days; (i) 6 days; (j) 8 days.Figure 11Results of target gene expression experiments after knockdown of Lamc1. (a) Alizarin red staining of osteoblasts induced for 3, 5 and 7 daysafter knockdown of Lamc1; (b) capture images with camera, 100\(\times \) microscope,scale bar \(=\) \(500\;\upmu \hbox {m}\); qPCR results of target gene expression after knockdown of Lamc1 in osteoblasts: (c) 3 days; (d) 5 days; (e) 7 days. (f) Oil Red O staining of adipocytes induced for 4, 6 and 8 days after knockdown of Lamc1; (g) capture images with camera, 100 \(\times \) microscope, scale bar \(=\) \(500\;\upmu \hbox {m}\); qPCR results of target gene expression after knockdown of Lamc1 in adipocyte: (h) 4 days; (i) 6 days; (j) 8 days.
