Abstract

   Bone marrow-derived mesenchymal stem cells (BMSCs) exhibit
   multi-lineage differentiation potential and robust proliferative
   capacity. The late stage of differentiation signifies the functional
   maturation and characterization of specific cell lineages, which is
   crucial for studying lineage-specific differentiation mechanisms.
   However, the molecular processes governing late-stage BMSC
   differentiation remain poorly understood. This study aimed to elucidate
   the key biological processes involved in late-stage BMSC
   differentiation. Publicly available transcriptomic data from human
   BMSCs were analyzed after approximately 14 days of osteogenic,
   adipogenic, and chondrogenic differentiation. Thirty-one differentially
   expressed genes (DEGs) associated with differentiation were identified.
   Pathway enrichment analysis indicated that the DEGs were involved in
   extracellular matrix (ECM)-receptor interactions, focal adhesion, and
   glycolipid biosynthesis, a ganglion series process. Subsequently, the
   target genes were validated using publicly available single-cell
   RNA-seq data from mouse BMSCs. Lamc1 exhibited predominant distribution
   in adipocytes and osteoblasts, primarily during the G2/M phase. Tln2
   and Hexb were expressed in chondroblasts, osteoblasts, and adipocytes,
   while St3gal5 was abundantly distributed in stem cells. Cell
   communication analysis identified two receptors that interact with
   LAMCI. q-PCR results confirmed the upregulation of Lamc1, Tln2, Hexb,
   and St3gal5 during osteogenic differentiation and their downregulation
   during adipogenic differentiation. Knockdown of Lamc1 inhibited
   adipogenic and osteogenic differentiation. In conclusion, this study
   identified four genes, Lamc1, Tln2, Hexb, and St3gal5, that may play
   important roles in the late-stage differentiation of BMSCs. It
   elucidated their interactions and the pathways they influence,
   providing a foundation for further research on BMSC differentiation.

   Keywords: Mesenchymal stem cell, Osteogenesis, Adipogenesis, Laminin,
   Single-cell RNA-seq

   Subject terms: Cell biology, Computational biology and bioinformatics,
   Developmental biology, Molecular biology, Stem cells

Introduction

   Mesenchymal stem cells (MSCs) are pluripotent that can be obtained from
   various tissues, umbilical cord, amniotic fluid, and
   placenta^[30]1–[31]4. It is capable of expanding and differentiating in
   vitro into a variety of mesodermal-type lineages, including bone,
   adipose, cartilage, muscle, tendon, and stroma, supporting
   hematopoiesis and the vascular system, and undergoes migration and
   paracrine signaling, which is not only crucial in embryonic development
   but also plays an important role in tissue homeostasis and repair
   throughout the life of an organism^[32]5,[33]6. MSCs from different
   sources exhibit different differentiation tendencies, and studies have
   shown that BMSCs exhibit greater APL activity, calcium deposition, and
   transformation capacity at an earlier age,compared with adipose-derived
   MSCs^[34]7. Therefore, BMSCs were chosen for studying the molecular
   mechanisms of osteogenic differentiation.

   Key transcription factors for osteoblast differentiation are required
   for the expression of osteoblast-specific genes, including
   Runt-associated transcription factor 2 (Runx2) and osteocalcin (Ocn),
   alkaline phosphatase (ALP), type I collagen, bone bridging proteins,
   and bone sialic acid protein^[35]8–[36]13. Bone morphogenetic proteins
   (BMPs), especially BMP-2, BMP-4, and BMP-7, activate the Smad signaling
   pathway, which in turn leads to the expression of Runx2 and Ocn and
   inhibits adipogenic differentiation^[37]14–[38]16. In addition to BMP,
   the Wnt/
   [MATH: <mi>β</mi> :MATH]
   -catenin signaling pathway plays a key role in osteoblast
   differentiation^[39]17. Activation of this pathway leads to the
   accumulation of
   [MATH: <mi>β</mi> :MATH]
   -catenin in the cytoplasm, which then translocates to the nucleus,
   activates the transcription of osteoblast-specific genes, and inhibits
   adipogenesis. On the other hand, adipocyte differentiation is regulated
   by peroxisome proliferator-activated receptor
   [MATH: <mi>γ</mi> :MATH]
   (PPAR
   [MATH: <mi>γ</mi> :MATH]
   )^[40]18. PPAR
   [MATH: <mi>γ</mi> :MATH]
   is considered a major regulator of adipogenesis because it controls the
   expression of many adipocyte-specific genes. PPAR
   [MATH: <mi>γ</mi> :MATH]
   activation not only promotes adipogenesis but also inhibits
   osteogenesis^[41]19. The balance between osteoblast and adipocyte
   differentiation is a tightly regulated process^[42]20. Therefore,
   understanding the mechanisms that regulate this balance is essential
   for developing therapeutic approaches for diseases such as osteoporosis
   and obesity, which are characterized by an imbalance between bone and
   fat formation.

   In mesenchymal stem cells, cellular matrix components play a key
   regulatory role in the directed differentiation of mesoderm into
   osteoblasts or adipogenic cells^[43]21. The extracellular matrix is a
   complex network of many molecules, including collagens, glycoproteins
   (GPs), and ECM-associated proteins^[44]22. Specific components of the
   ECM can influence cell fate decisions; for instance, collagen promotes
   osteoblast differentiation, whereas fibronectin may promote adipogenic
   cell differentiation^[45]23,[46]24. Cell adhesion molecules and cell
   signaling molecules in the cell matrix are important regulators of the
   differentiation process, transmitting signals and modulating cell
   behavior by interacting with receptors on the cell
   surface^[47]25,[48]26. For example, integrins increase the ability of
   cells to differentiate osteoblastically, maintain bone homeostasis, and
   regulate bone mass while inhibiting the adipogenic differentiation of
   BMSCs^[49]27–[50]29. Taken together, the study of cellular matrix
   components contributes to our understanding of the mechanisms of germ
   layer-directed differentiation.

   BMSCs exhibit the most pronounced difference between days 14 and 17 of
   differentiation, and the difference between days 17 and 21 of
   differentiation is not significant^[51]7. Therefore, our study combined
   tri-lineage differentiation with bulk RNA-seq analysis and single-cell
   RNA-seq (scRNA-seq) analysis to determine which factors simultaneously
   play key roles in osteogenic, adipogenic, and chondrogenic terminal
   phase differentiation, and cell staining and RT-qPCR were used to
   verify the changes in target genes during differentiation, aiming to
   provide a new basis for MSC therapies.

Results

hBMSCs osteogenic, adipogenic, and chondrogenic differentiation 14-day
transcriptome data cross-identify 31 DEGs

   A preliminary comparison of the six datasets ([52]GSE36923,
   [53]GSE44303, [54]GSE109503, [55]GSE140861, [56]GSE28205, [57]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. [58]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. [59]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. [60]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.
   [61]2j–o).

Figure 1.

   [62]Figure 1
   [63]Open in a new tab

   Quality of RNA-seq data of bone marrow-derived MSCs at the late stage
   of differentiation. Box plot: (a) osteogenic differentiation dataset
   ([64]GSE28205, [65]GSE37558); (b) adipogenic differentiation dataset
   ([66]GSE36923, [67]GSE44303); (c) chondrogenic differentiation dataset
   ([68]GSE109503, [69]GSE140861); PCA plot; (d) [70]GSE37558; (e)
   [71]GSE36923; (f) [72]GSE140861; (g) [73]GSE28205; (h) [74]GSE44303;
   (i) [75]GSE109503.

Figure 2.

   [76]Figure 2
   [77]Open in a new tab

   Target 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 ([78]GSE28205 and [79]GSE37558); (c,d) volcano
   plots of differential gene expression in adipogenic differentiation
   datasets ([80]GSE36923 and [81]GSE44303); (e,f) volcano plots of
   differential gene expression in cartilageogenic differentiation
   datasets ([82]GSE109503 and [83]GSE140861) volcano plots of
   differential gene expression; (g) heatmap of differential gene
   expression of osteogenic differentiation datasets ([84]GSE28205 and
   [85]GSE37558) after integration by the RRA algorithm; (h) heatmap of
   differential gene expression of adipogenic differentiation datasets
   ([86]GSE36923 and [87]GSE44303) after integration of by the RRA
   algorithm; (i) heatmap of differential gene expression of chondrogenic
   differentiation dataset ([88]GSE109503 and [89]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 interactions

   Gene 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. [90]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. [91]2q).

ceRNA network construction

   Using 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. [92]3a,b). The
   gene symbols, abbreviations and functions are shown in Table [93]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 3.

   [94]Figure 3
   [95]Open in a new tab

   Expression 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.
   Gene symbol Describe Function
   Col4a1 Collagen Type IV Alpha 1 Chain Type IV collagen proteins are
   integral components of basement membranes
   Lamc1 Laminin Subunit Gamma 1 Laminins, a family of extracellular
   matrix glycoproteins, are the major noncollagenous constituent of
   basement membranes. They have been implicated in a wide variety of
   biological processes including cell adhesion, differentiation,
   migration, signaling, neurite outgrowth, and metastasis
   Itga10 Integrin Alpha-10 They participate in cell adhesion as well as
   cell-surface-mediated signaling
   Hapln1 Hyaluronan And Proteoglycan Link Protein 1 It is predicted to
   enable hyaluronic acid binding activity, to be an extracellular matrix
   structural constituent conferring compression resistance. And it is
   predicted to be involved in skeletal system development
   Col10a1 Collagen Type X Alpha 1 Chain This gene is a short chain
   collagen expressed by hypertrophic chondrocytes during endochondral
   ossification
   Tln2 Talin 2 It is thought to associate with unique transmembrane
   receptors to form novel linkages between extracellular matrices and the
   actin cytoskeleton
   Hexb Hexosaminidase Subunit Beta Related to hydrolase activity
   St3gal5 ST3 Beta-Galactoside AlphA-2,3-Sialyltransferase 5 The protein
   encoded by this gene is a type II membrane protein
   [96]Open in a new tab

Target gene expression was higher at both the adipose gene level and tissue
level than at bone tissue

   Analysis 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. [97]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. [98]3e).

BMSCs are categorized into 11 clusters at single-cell resolution

   We used t-distributed stochastic neighbor embedding (tSNE) to
   investigate cellular heterogeneity within clusters, and after the
   integration of the two datasets ([99]GSE128423, [100]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. [101]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.
   [102]4c.

Figure 4.

   [103]Figure 4
   [104]Open in a new tab

   Subcluster 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
reliability

   We 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. [105]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. [106]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. [107]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. [108]5e).

Figure 5.

   [109]Figure 5
   [110]Open in a new tab

   Proof 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 genes

   Based 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. [111]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
   [MATH: <mi>γ</mi> :MATH]
   and Fabp4 are given as reference (Fig. [112]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. [113]8).

Figure 6.

   [114]Figure 6
   [115]Open in a new tab

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

   [116]Figure 7
   [117]Open in a new tab

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

   [118]Figure 8
   [119]Open in a new tab

   Expression 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 Lamc1

   In 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. [120]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. [121]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. [122]9c–e).

Figure 9.

   [123]Figure 9
   [124]Open in a new tab

   The 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
differentiation

   The 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. [125]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. [126]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. [127]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
   [MATH: <mi>γ</mi> :MATH]
   gradually increased to a maximum on 8 days (Fig. [128]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
   [MATH: <mi>γ</mi> :MATH]
   and Fabp4 were also downregulated (Fig. [129]11c–e,h–j).Similarly, the
   red color of the knocked-down cells was significantly lighter than that
   of the controls (Fig. [130]11a,b,f,g). Similarly, the previous
   experiment was repeated with the housekeeping gene GAPDH and yielded
   similar results (Supplementary Fig. [131]S1).

Figure 10.

   [132]Figure 10
   [133]Open in a new tab

   Results of target gene expression experiments. (a) Alizarin red
   staining of osteoblasts induced for 3, 5 and 7 days; (b) capture images
   with camera,100
   [MATH: <mo>×</mo> :MATH]
   microscope,scale bar
   [MATH: <mrow><mo>=</mo><mn>500</mn><mspace
   width="0.277778em"></mspace><mi
   mathvariant="normal">μ</mi><mtext>m</mtext></mrow> :MATH]
   ; 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
   [MATH: <mo>×</mo> :MATH]
   microscope, scale bar
   [MATH: <mo>=</mo> :MATH]
   [MATH: <mrow><mn>500</mn><mspace width="0.277778em"></mspace><mi
   mathvariant="normal">μ</mi><mtext>m</mtext></mrow> :MATH]
   ; qPCR results of target gene expression in adipocyte: (h) 4 days; (i)
   6 days; (j) 8 days.

Figure 11.

   [134]Figure 11
   [135]Open in a new tab

   Results 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
   [MATH: <mo>×</mo> :MATH]
   microscope,scale bar
   [MATH: <mo>=</mo> :MATH]
   [MATH: <mrow><mn>500</mn><mspace width="0.277778em"></mspace><mi
   mathvariant="normal">μ</mi><mtext>m</mtext></mrow> :MATH]
   ; 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
   [MATH: <mo>×</mo> :MATH]
   microscope, scale bar
   [MATH: <mo>=</mo> :MATH]
   [MATH: <mrow><mn>500</mn><mspace width="0.277778em"></mspace><mi
   mathvariant="normal">μ</mi><mtext>m</mtext></mrow> :MATH]
   ; qPCR results of target gene expression after knockdown of Lamc1 in
   adipocyte: (h) 4 days; (i) 6 days; (j) 8 days.

Methods

Data collection and processing

   Gene expression profiling data ([136]GSE36923, [137]GSE44303,
   [138]GSE109503, [139]GSE140861, [140]GSE28205 and [141]GSE37558) were
   obtained from the Gene Expression Omnibus (GEO). [142]GSE109503
   contained RNA-seq data, and the remaining five gene sets were obtained
   from bulk RNA-seq data. [143]GSE109503 and [144]GSE140861 are
   chondrogenic differentiation datasets; [145]GSE109503 has four subsets
   (two controls and two experimental groups); [146]GSE140861 has six
   subsets (three controls and three experimental groups); [147]GSE28205
   and [148]GSE37558 are osteogenic differentiation datasets;
   [149]GSE28205 has six subsets (three control groups and three
   experimental groups); [150]GSE37558 has six subsets (three control
   groups and three experimental groups); [151]GSE36923 and [152]GSE44303
   are adipogenic differentiation datasets; [153]GSE36923 has seven
   subsets (four control groups and three experimental groups); and
   [154]GSE44303 has seven subsets (four controls and three experimental
   groups).

   Single-cell RNA-seq data ([155]GSE128423 and [156]GSE156635) were
   obtained from the Gene Expression Omnibus (GEO) based on a 10
   [MATH: <mo>×</mo> :MATH]
   Genomics assay with eight normal BMSC tissue samples. The summary
   results of the gene expression profiling dataset information are
   provided in Table [157]2.

Table 2.

   Summary information on gene expression profiling datasets.
   DataSet Samples Organism Type of cell Differentiation time Group DOI
   [158]GSE36923

   [159]GSM906367,

   [160]GSM906368,

   [161]GSM906369.
   Homo sapiens BMSCs 0 Day Control 10.1016/j.biocel.2013.11.010

   [162]GSM906370,

   [163]GSM906371,

   [164]GSM906372.
   15 Days Adipogenic
   [165]GSE44303

   [166]GSM1082817,

   [167]GSM1082818,

   [168]GSM1082819.
   Homo sapiens BMSCs 0 Day Control 10.1093/nar/gku567

   [169]GSM1082844,

   [170]GSM1082845,

   [171]GSM1082846.
   14 Days Adipogenic
   [172]GSE109503

   [173]GSM2944707,

   [174]GSM2944708,

   [175]GSM2944709.
   Homo sapiens BMSCs 0 Day Control 10.1096/fj.201800534R

   [176]GSM1082844,

   [177]GSM1082845,

   [178]GSM1082846.
   14 Days Chondrocyte
   [179]GSE140861

   [180]GSM4188963,

   [181]GSM4188970.
   Homo sapiens BMSCs 0 Day Control 10.1089/ten.TEA.2016.0559

   [182]GSM4188967,

   [183]GSM4188974.
   14 Days Chondrocyte
   [184]GSE28205

   [185]GSM698427,

   [186]GSM698433,

   [187]GSM698439.
   Homo sapiens BMSCs 0 Day Control 10.1002/jbmr.1578

   [188]GSM698429,

   [189]GSM698431,

   [190]GSM698436.
   14 Days Osteoblast
   [191]GSE37558

   [192]GSM921574,

   [193]GSM921575,

   [194]GSM921576,

   [195]GSM921578.
   Homo sapiens BMSCs 0 Day Control 10.1186/1471-2164-15-965

   [196]GSM921584,

   [197]GSM921585,

   [198]GSM921586.
   12 Days Osteoblast
   [199]GSE128423

   [200]GSM3674224,

   [201]GSM3674225,

   [202]GSM3674226,

   [203]GSM3674227,

   [204]GSM3674228,

   [205]GSM3674229.
   Mus musculus BMSCs Undifferentiation – 10.1016/j.cell.2019.04.040
   [206]GSE156635

   [207]GSM4735397,

   [208]GSM4735398.
   Mus musculus BMSCs Undifferentiation – 10.1016/j.celrep.2021.109352
   [209]Open in a new tab

   Microarray data were normalized ,and variance was analyzed using the
   Limma package in the R/Bioconductor software for matrix data from each
   GEO dataset, and RNA-seq data were normalized ,as was variance analyzed
   using the edgeR package. DEGs of the same cell type identified from
   each of the six datasets were integrated based on different
   differentiation categories using the RobustRankAggreg package. |log2FC|
   [MATH: <mo>≥</mo> :MATH]
   1 and P
   [MATH: <mo>≤</mo> :MATH]
   0.05 were considered to indicate statistical significance.

   PCA plots for each dataset were plotted using the FactoMineR package
   and the factoextra package in the R/Bioconductor software. pheatmap
   package was used to plot heat maps of the genes after data merging.
   ggplot2 package was used to plot volcano plots after data merging. Venn
   plots were generated using TBtools software (version x64 v1.09876). GO
   function enrichment and KEGG pathway enrichment were performed for 31
   differential genes (DEGs) using the clusterProfiler package in R
   software. adj. P
   [MATH: <mo>≤</mo> :MATH]
   0.05 was considered statistically significant. Bar and bubble plots
   were drawn using the ggplot2 package, and network plots were drawn
   using the enrichplot package.

   scRNA-seq data were preprocessed using the Seurat package in R software
   for matrix data from each GEO dataset, the PercentageFeatureset
   function was used to determine the proportion of mitochondrial genes,
   and correlation analyses were used to investigate the relationships
   between sequencing depth and mitochondrial gene sequences and total
   intracellular sequences. Each gene was expressed in at least 3 cells.
   The gene expression in each cell was greater than 500 and less than
   6000, the mitochondrial content was less than 30%, and the UMI in each
   cell was at least greater than 200. The scRNA-seq data were normalized
   by the NormalizeData method after data filtering, and highly variable
   genes were identified by the FindVariableFeatures method to compute the
   mean expression and dispersion and converted to logarithmic format.
   Next, the data from the two scRNA-seq datasets were integrated using
   the Harmony package, and then we performed principal component analysis
   (PCA) and reduced the data to the first 17 PCA components (the number
   of components selected based on the standard deviation of the principal
   components—in the platform area of the Elbow Graph). To visualize the
   data, we used the Seurat package for t-distribution random neighbor
   embedding (tSNE) to project the cells in two-dimensional space. We
   visualized clusters on 2D maps generated using t-distributed random
   neighbor embedding (t-SNE). Clusters were annotated according to
   CellMarker 2.0 ([210]http://yikedaxue.slwshop.cn/), and the annotated
   clusters were scored using the AUCell package, selecting the WNT
   pathway-a key pathway for osteogenesis-and the ADIPOGENESIS pathway-a
   key pathway for adipogenesis. The cell cycle of the annotated clusters
   was inferred using the tricycle package, and the cell cycle in which
   the target genes are located was determined. The Cytotrace package was
   used to determine the stemness and differentiation potential of stem
   cells, chondrocytes, osteoblasts, and adipocytes, and based on the
   results of Cytotrace, Monocle2 was used to determine the direction of
   differentiation of the four types of cells and to determine the site of
   expression of the target gene.

   The CellChat package was used to query the role of intercellular
   communication, and the ECM-Receptor database was selected for
   subsequent analysis based on the KEGG enrichment results. The identify
   OverExpressedGenes function of the CellChat package was used to
   identify highly expressed genes, the identify OverExpressedInteractions
   function was used to identify highly expressed pathways, the project
   Data function was used to project the relationship between the highly
   expressed genes and the pathways to the PPI network, the
   computeCommunProb function was used to calculate the probability of
   cell-to-cell communication, data with the minimum number of cells less
   than 10 were removed, the computeCommunProbPathway function was used to
   calculate the cell-to-cell communication at the level of signaling
   pathways, the aggregateNet function was used to calculate the
   aggregated cell-to-cell communication network, and the “LAMININ”
   signaling pathway was extracted.

   The mechanism map was drawn via Figdraw
   ([211]https://www.figdraw.com/static/index.html).

PPI network construction

   The String website ([212]https://cn.string-db.org/) analyzed protein
   functional interactions and selected interactions with a composite
   score>0.4. The RNAInter database ([213]http://www.rnainter.org/) is a
   complete source of RNA interaction data for studying the relationships
   between mRNAs, lncRNAs, and microRNAs. The RNAInter database was used
   to identify miRNAs that interact with mRNAs, as well as lncRNAs that
   interact with miRNAs, and interactions with a score > 0.5 was selected.
   Cytoscape (version v3.9.1) is an open-source bioinformatics software
   platform for visualizing molecular interaction networks. The
   relationship between target genes and non-coding RNAs was plotted using
   Cytoscape.

Tissue distribution and gene expression identification

   Protein and gene expression levels in bone marrow and adipose tissue
   were compared with those in the HPA database using The Human Protein
   Atlas ([214]https://www.proteinatlas.org/) as a validation tool to
   determine whether there are differences in the expression of these
   genes in osteogenic and adipogenic differentiation. Expression scores
   in abdominal adipose tissue, omental fat pads, subcutaneous adipose
   tissue, tibia, synovial tissue, synovial joints, cartilage tissue,
   trabecular bone tissue, and bone marrow tissue were used to validate
   whether there are differences in the expression of these genes in
   osteogenic and adipogenic as well as cartilaginous differentiation
   using the BGEE database ([215]https://bgee.org/). Bar graphs were drawn
   using GraphPad Prism 9.

Cell culture and differentiation

   MC3T3-E1 cells (Procell Life, Science & Technology, China) were
   cultured in MEM
   [MATH: <mi>α</mi> :MATH]
   medium supplemented with 10% fetal bovine serum (FBS), 100 units/ml
   penicillin and streptomycin at
   [MATH: <mrow><mn>37</mn><mmultiscripts><mspace
   width="0.277778em"></mspace><mrow></mrow><mo>∘</mo></mmultiscripts><mi>
   C</mi></mrow> :MATH]
   in a 5%
   [MATH: <msub><mtext>CO</mtext><mn>2</mn></msub> :MATH]
   incubator. The medium was changed every 2–3 days and the cells were
   used for experiments in which the cells were grown to the logarithmic
   growth phase. For osteogenic differentiation, MC3T3-E1 cells were
   cultured in osteogenic medium supplemented with 10 mM
   [MATH: <mi>β</mi> :MATH]
   -glycerophosphate and 50
   [MATH: <mi mathvariant="normal">μ</mi> :MATH]
   g/mL L ascorbic acid osteogenic inducer. The cells were induced for 3,
   5, or 7 days. 3T3-L1 cells (Procell Life, Science & Technology, China)
   were cultured in DMEM medium supplemented with 10% fetal bovine serum
   (FBS), 100 units/ml penicillin and streptomycin at
   [MATH: <mrow><mn>37</mn><mmultiscripts><mspace
   width="0.277778em"></mspace><mrow></mrow><mo>∘</mo></mmultiscripts><mi>
   C</mi></mrow> :MATH]
   in a 5%
   [MATH: <msub><mtext>CO</mtext><mn>2</mn></msub> :MATH]
   incubator.The medium was changed every 2-3 days, and the cells were
   used for experiments in which the cells were grown to the logarithmic
   growth phase. For adipogenic differentiation, 3T3-L1 cells were
   cultured in adipogenic induction medium supplemented with 0.5
   [MATH: <mi mathvariant="normal">μ</mi> :MATH]
   M 3-isobutyl-1-methyl xanthine , 1
   [MATH: <mi mathvariant="normal">μ</mi> :MATH]
   M dexamethasone and 10
   [MATH: <mi mathvariant="normal">μ</mi> :MATH]
   g/mL insulin adipogenic inducer for 2 days, and the medium was replaced
   with adipogenic maintenance medium supplemented with 10
   [MATH: <mi mathvariant="normal">μ</mi> :MATH]
   g/mL insulin for 2 days. The medium was finally replaced with DMEM, and
   the culture was continued for 2–4 days.

Alizarin red staining method and Oil Red O staining

   After fixation with 80% ethanol for 30 min, cells were incubated with
   fresh alizarin red solution (Cyagen) for 10 min. After fixation with 4%
   neutral formaldehyde for 30 min, cells were incubated with saturated
   Oil Red O staining solution (Beyotime Biotechnology; Shanghai; China)
   for 10 min.

Quantitative reverse transcription PCR

   Total RNA was extracted using triazole lysate lysis (Sangon Biotech;
   Shanghai; China), isopropanol, chloroform, and 75% ethanol solutions.
   Total RNA reverse transcription was induced using the HiScirpt Reverse
   Transcriptase kit (Novozymes Bio; Nanjing; China) The genes and primers
   used in the osteogenesis and adipogenesis process are listed in Table
   [216]3. Performed with a QuantStudio 5 real-time fluorescent
   quantitative PCR system (Applied Biosystems; Shanghai; China). Each
   sample was performed in triplicate. The sequences of the mRNAs were
   calculated according to the
   [MATH: <mrow><mi mathvariant="normal">Δ</mi><mi
   mathvariant="normal">Δ</mi><msub><mtext>C</mtext><mtext>T</mtext></msub
   ></mrow> :MATH]
   relative stacking method. mRNAs were sequenced as Table [217]3.

Table 3.

   qPCR primer sequences.
        Gene                     Forward primers         Reverse primers
        Lamc1                    TGCCGGAGTTTGTTAATGCC    TGGTTGTTGTAGTCGGTCAGG
        Tln2                     AGGAAAGGGATTTGGCTAGAAGC CCAACATTCGGATTTTCTGAGGC
        Hexb                     CTGGTGTCGCTAGTGTCGC     CAGGGCCATGATGTCTCTTGT
        St3gal5                  AAAAGAATGCACTATGTGGACCC ATAGCCGTCTTCGCGTACCT
        Runx2                    CGGTCTCCTTCCAGGATGGT    GCTTCCGTCAGCGTCAACA
        Ocn                      GAACAGACAAGTCCCACACAGC  TCAGCAGAGTGAGCAGAAAGAT
   Ppar
   [MATH: <mi>γ</mi> :MATH]
        GGAAGACCACTCGCATTCCTT    GTAATCAGCAACCATTGGGTCA
        Fabp4                    ATCAGCGTAAATGGGGATTTGG  GTCTGCGGTGATTTCATCGAA
   [MATH: <mi>β</mi> :MATH]
 -actin GTGACGTTGACATCCGTAAAGA   GCCGGACTCATCGTACTCC
   [218]Open in a new tab

Lamc1 siRNA interference

   Lamc1 siRNA was obtained from GenePharma Corporation. siRNA was
   transfected into MC3T3-El cells using the gp transfection mate
   (GenePharma Corporation) according to the manufacturer’s instructions.
   24 h later, induction was initiated, and cells were collected 3–7 days
   later for alizarin red staining, and real-time quantitative RT PCR
   analysis. The siRNA was transfected into 3T3-L1 cells using gp
   transfection mate (GenePharma Corporation) according to the
   manufacturer’s instructions. After 90% confluence, induction was
   initiated, and the cells were collected 4–8 days later for oil red o
   staining, and real-time quantitative qPCR analysis experiments. The
   sequence of the siRNA was as follows. NC:5
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   -UUCUCCGAACGUGUCACGUTT-3
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   and 5
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   -ACGUGACACGUUCGGAGAATT-3
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   ; Lamc1-siRNA: 5
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   -GCCGUAAUCUCAGACAGUUTT-3
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   and 5
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   - AACUGUCUGAGAUUACGGCTT-3
   [MATH:
   <mmultiscripts><mrow></mrow><mrow></mrow><mo>′</mo></mmultiscripts>
   :MATH]
   .

Molecular docking analysis

   The Hdock website ([219]http://hdock.phys.hust.edu.cn/) was used to
   perform simulated docking studies of Laminin-111 (Entry ID:5mc9)with
   DAG1 (Entry ID: 5llk) or CD44 (Entry ID:4pz3), the PymoL software
   (version v3.0.0) was used to perform pre-docking treatments (water
   molecules from the structure and co-crystallized ligands were removed,
   H atoms were added) and Pymol was used for the visual presentation of
   the docking complexes, with the results with the highest docking scores
   being selected.

Discussion

   To investigate the molecular mechanisms of the multidirectional
   differentiation of BMSCs, especially osteogenesis, adipogenesis and
   chondrogenesis, BMSCs exhibited significant cellular heterogeneity, and
   an integrated analysis of the multiscale transcriptome was performed in
   the present study. These DEGs were analyzed by bulk RNA-seq and
   scRNA-seq, and it was found that the late stage of differentiation was
   the most significant in terms of interactions with the ECM receptor. It
   is well known that the ECM enhances cell recruitment through cell
   surface receptors, which determine cellular interactions with the ECM
   and trigger specific cellular functions such as adhesion, migration,
   proliferation, and differentiation^[220]30. Several studies have
   reported the critical role of the matrix in regulating the
   differentiation of MSCs toward the bone or adipose
   lineage^[221]30,[222]31.Recent studies have shown that stem cells can
   mechanically sense the hardness of their microenvironment, and the
   direction of differentiation of MSCs may change with changes in the
   matrix composition^[223]32,[224]33. For example, the extracellular
   matrix supplemented with glucan sulfate promotes the osteogenic
   differentiation of mesenchymal stem cells^[225]34, and type I collagen
   inhibits adipogenesis and differentiation by activating yes-associated
   proteins^[226]35. The latest research by Yangzi Jiang et al. showed
   that ECM produced by adult stem cells from different tissues has
   different cartilage induction abilities^[227]36. Osteoblasts are the
   main bone-forming cells, and they produce extracellular proteins that
   constitute the major components of bone. During bone formation, BMSCs
   differentiate into osteoblasts, which produce bone by synthesizing an
   extracellular matrix composed of various proteins. The initial
   deposition of the extracellular matrix is called the osteoid, which is
   then mineralized into bone tissue by the accumulation of hydroxyapatite
   (
   [MATH:
   <mrow><msub><mtext>Ca</mtext><mn>1</mn></msub><mn>0</mn><msub><mrow><mo
   stretchy="false">(</mo><msub><mtext>PO</mtext><mn>4</mn></msub><mo
   stretchy="false">)</mo></mrow><mn>6</mn></msub><msub><mrow><mo
   stretchy="false">(</mo><mtext>OH</mtext><mo
   stretchy="false">)</mo></mrow><mn>2</mn></msub></mrow> :MATH]
   ) as calcium phosphate^[228]37. As osteoblasts mature, they transform
   into osteocytes, and revealing the detailed mechanisms that regulate
   the later stages of osteoblast differentiation and function is
   important for clinical applications^[229]38. However, studies on the
   mechanisms underlying the differentiation of bone marrow-derived MSCs
   in the late stages of differentiation still need to be performed. We
   therefore selected an mRNA microarray dataset at the end of
   differentiation to identify genes that showed significant heterogeneity
   in the differentiation of BMSCs into adipocytes, osteoblasts, and
   chondrocytes. lamc1 tended to increase osteoblastic differentiation and
   decrease adipocytes, and Col4a1 had the opposite effect on
   differentiation. Highly expressed Lamc1 and Col4a1 play key roles in
   the protein interaction network, and their interactions with the key
   molecules involved in osteogenic, adipogenic, and chondrogenic
   differentiation, H19, XIST, NEAT1, and MALAT1, have been identified.
   According to our results, Hui-jian Chen et al. reported that Lama4 and
   Col4a1 are expressed in white adipose tissue at different sites and are
   upregulated during adipocyte differentiation^[230]39. Interestingly,
   the mRNA expression of Col4a1 was significantly increased after
   adipogenic induction of MSCs in a recent study, which is consistent
   with our expectations^[231]40. More importantly, laminin in the cell
   matrix contains three subunits:
   [MATH: <mi>α</mi> :MATH]
   ,
   [MATH: <mi>β</mi> :MATH]
   , and
   [MATH: <mi>γ</mi> :MATH]
   . Lamc1 is the most commonly used
   [MATH: <mi>γ</mi> :MATH]
   chain and is present in most laminin molecules^[232]41. During the
   process of osteogenesis, Laminin activates downstream signaling
   pathways such as Rho GTPase and MAPK by interacting with Integrin,
   which in turn regulates the processes of cellular bone morphogenesis
   and bone matrix deposition, and participates in the directional
   differentiation of stem cells to osteoblasts^[233]42,[234]43.Lamc1 is
   present in LN-111 and LN-511, it has been shown to promote osteogenic
   differentiation by binding to integrins^[235]41,which are key factors
   for the development of osteoblasts. The HPA database, Bgee database,
   and single-cell transcriptome data analysis revealed that Lamc1
   expression was significantly greater in adipocytes than in osteoblasts,
   but Lamc1 tended to promote bone formation and inhibit adipogenic, and
   this tendency became more obvious with increasing induction time.
   Further studies showed that after Lamc1 knockdown, the expression of
   osteogenic and adipogenic markers was significantly decreased, and
   differentiation was inhibited. It is therefore reasonable to assume
   that Lamc1 plays an important role in regulating the balance of bone
   adipogenesis^[236]44.

   The distribution of target genes in bone marrow-derived MSCs, the cell
   cycle, and ligand–receptor interactions were determined by single-cell
   transcriptome analysis. Among these genes, Lamc1, which is expressed
   mainly during the G2/M phase in adipogenic cells and may affect
   intracellular signaling through interactions with Dag1, was used as a
   focus for determining the balance between osteogenic and adipogenic
   differentiation. Dag1 (also known as
   [MATH: <mi>α</mi> :MATH]
   -glycoprotein) is involved in the adhesion between the extracellular
   matrix and cells, helps to maintain cell aggregation and structural
   stability, and participates in the assembly process of laminin and the
   basement membrane^[237]45. Recent studies have shown that Dag1 plays a
   role in cartilage formation and osteoblast differentiation in
   vivo^[238]46,[239]47.Another finding is that the interaction of Dag1
   with laminin has been reported^[240]48.To prove our prediction,
   molecular docking experiments were performed on laminin and Dag1. The
   experimental results showed that laminin and Dag1 formed a total of 11
   hydrogen bonds (Fig. [241]12a). Our data support the conclusion that
   the interaction of Dag1 and Lamc1 may influence the balance between the
   adipogenic and osteogenic differentiation of MSCs through cell
   adhesion, signaling, and gene regulation (Fig. [242]12c). Similarly, a
   total of four hydrogen bonds were formed after molecular docking
   experiments were performed on laminin and CD44, initially validating
   the interaction between Lamc1 and CD44 (Fig. [243]12b).

Figure 12.

   [244]Figure 12
   [245]Open in a new tab

   Molecular docking and mechanism maps (a) Results of molecular docking
   of laminin to DAG1; (b) Results of molecular docking of laminin to
   CD44; (c) mechanism maps.

   In conclusion, this study aimed to identify DEGs that may be involved
   in osteogenic differentiation, adipogenic differentiation, and the
   onset and progression of chondrogenic differentiation of bone marrow
   mesenchymal stem cells (BMSCs). A total of eight hub genes were
   identified, four of which were discovered and proposed for the first
   time in this study. Although these findings are preliminary, they
   suggest the existence of previously unknown regulatory relationships
   governing the differentiation process of BMSCs. Unfortunately, our
   experiments did not reveal new regulatory relationships between
   chondrogenesis and adipogenesis or osteogenesis, which may be related
   to our choice of a more conservative data merging approach.However,
   this conservative strategy ensured the accuracy of bioinformatic
   predictions. Despite these limitations, the identified hub genes and
   their potential interactions provide valuable insights and lay the
   foundation for further investigation into the molecular mechanisms
   underlying the multi-lineage differentiation of BMSCs.

Conclusion

   This study aimed to identify the key biological processes involved in
   the late stages of differentiation of BMSCs. The following conclusions
   were demonstrated: eight target genes were identified in the context of
   trilineage differentiation antagonism, which is concentrated in the
   extracellular matrix and is associated with ECM-receptor
   interactions.Notably, Lamc1 expression was higher in adipose tissue
   than in bone tissue. Furthermore, Lamc1 exhibited maximal expression
   during the mitotic phase of the cell cycle. Bioinformatic analysis
   suggested that Lamc1 interacts with CD44 and may influence the
   differentiation of stem cells towards osteoblastic, chondrogenic, or
   adipogenic lineages. Additionally, the interaction between Lamc1 and
   Dag1 may affect the balance of trilineage differentiation in BMSCs.
   Experimental validation confirmed that Lamc1 was upregulated during
   osteogenic differentiation and downregulated during adipogenic
   differentiation, indicating its essential role in regulating the
   differentiation process.

Supplementary Information

   [246]Supplementary Information.^ (1MB, pdf)

Author contributions

   L.Z. and S.L. wrote the main manuscript text and Y.P. typesetted and
   checked the manuscripts. J.Z. led and supervised the planning and
   execution of the experiments. All authors reviewed the manuscript.

Funding

   This research was supported by the National Natural Science Foundation
   of China (Grant numbers 81960171 and 82360167) and Guizhou Provincial
   Science and Technology Program (Guizhou Science and Technology
   Foundation-ZK[2022] General 620).

Data availibility

   The datasets covered in this reasearch are all publicly available
   transcriptomic data from the Gene Expression Omnibus (GEO). Ullah M,
   Sittinger M, Ringe J. Transdifferentiation of adipogenically
   differentiated cells into osteogenically or chondrogenically
   differentiated cells: phenotype switching via dedifferentiation. Int J
   Biochem Cell Biol. 2014 Jan;46:124-37. doi:
   10.1016/j.biocel.2013.11.010. Epub 2013 Nov 22. [247]GSE36923. Tomaru
   Y, Hasegawa R, Suzuki T, Sato T, Kubosaki A, Suzuki M, Kawaji H,
   Forrest AR, Hayashizaki Y; FANTOM Consortium; Shin JW, Suzuki H. A
   transient disruption of fibroblastic transcriptional regulatory network
   facilitates trans-differentiation. Nucleic Acids Res. 2014
   Aug;42(14):8905-13. doi: 10.1093/nar/gku567. Epub 2014 Jul 10.
   [248]GSE44303. Huynh NPT, Zhang B, Guilak F. High-depth transcriptomic
   profiling reveals the temporal gene signature of human mesenchymal stem
   cells during chondrogenesis. FASEB J. 2019 Jan;33(1):358-372. doi:
   10.1096/fj.201800534R. Epub 2018 Jul 9. [249]GSE109503. Somoza RA,
   Correa D, Labat I, Sternberg H, Forrest ME, Khalil AM, West MD, Tesar
   P, Caplan AI. Transcriptome-Wide Analyses of Human Neonatal Articular
   Cartilage and Human Mesenchymal Stem Cell-Derived Cartilage Provide a
   New Molecular Target for Evaluating Engineered Cartilage. Tissue Eng
   Part A. 2018 Feb;24(3-4):335-350. doi: 10.1089/ten.TEA.2016.0559. Epub
   2017 Jul 28.[250]GSE140861. Dani N, Olivero M, Mareschi K, van Duist
   MM, Miretti S, Cuvertino S, Patané S, Calogero R, Ferracini R,
   Scotlandi K, Fagioli F, Di Renzo MF. The MET oncogene transforms human
   primary bone-derived cells into osteosarcomas by targeting committed
   osteo-progenitors. J Bone Miner Res. 2012 Jun;27(6):1322-34. doi:
   10.1002/jbmr.1578. [251]GSE28205. Alves RD, Eijken M, van de Peppel J,
   van Leeuwen JP. Calcifying vascular smooth muscle cells and
   osteoblasts: independent cell types exhibiting extracellular matrix and
   biomineralization-related mimicries. BMC Genomics. 2014 Nov
   7;15(1):965. doi: 10.1186/1471-2164-15-965. [252]GSE37558. Baryawno N,
   Przybylski D, Kowalczyk MS, Kfoury Y, Severe N, Gustafsson K,
   Kokkaliaris KD, Mercier F, Tabaka M, Hofree M, Dionne D, Papazian A,
   Lee D, Ashenberg O, Subramanian A, Vaishnav ED, Rozenblatt-Rosen O,
   Regev A, Scadden DT. A Cellular Taxonomy of the Bone Marrow Stroma in
   Homeostasis and Leukemia. Cell. 2019 Jun 13;177(7):1915-1932.e16. doi:
   10.1016/j.cell.2019.04.040. Epub 2019 May 23. [253]GSE128423. Sivaraj
   KK, Jeong HW, Dharmalingam B, Zeuschner D, Adams S, Potente M, Adams
   RH. Regional specialization and fate specification of bone stromal
   cells in skeletal development. Cell Rep. 2021 Jul 13;36(2):109352. doi:
   10.1016/j.celrep.2021.109352.[254]GSE156635.

Competing interests

   The authors declare no competing interests.

Footnotes

   Publisher's note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

   These authors contributed equally: Lixia Zhao and Shuai Liu.

Supplementary Information

   The online version contains supplementary material available at
   10.1038/s41598-024-69629-4.

References