Transcriptome Sequencing Reveals Key Genes in Three Early Phases of Osteogenic, Adipogenic, and Chondrogenic Differentiation of Bone Marrow Mesenchymal Stem Cells in Rats

Introduction

The microenvironment within the bone marrow cavity consists of both stromal lineage and hematopoietic lineage. Mesenchymal stem cells (MSCs) are discovered by Friedenstein et al. from stromal tissues of the bone marrow (BM) more than 50 years ago and have the potential to give rise to multiple mesodermal progenies such as osteocytes, adipocytes, and chondrocytes (Arthur and Gronthos, 2020; Jiang et al., 2021; Johnson et al., 2021; Saeedi et al., 2021). Meanwhile, bone MSCs (BMSCs) possess immune-modulatory and anti-inflammatory abilities (Li et al., 2012). Therefore, cell therapy with BMSCs has been proposed as an eligible paradigm for treating trauma, bone defects, and bone metabolic disorders (e.g., osteonecrosis and osteonecrosis) (Galderisi et al., 2021; Mabuchi et al., 2021). However, to date, clinical trials of MSC transplantation for osteoporosis predominantly focused on applying autologous cells, and no results have been reported (Jiang et al., 2021). Many unsolved problems and challenges still exist (Galderisi et al., 2021; Mabuchi et al., 2021; Saeedi et al., 2021). The first one is the patient selection and source of stem cells. A small number of MSCs can be readily obtained for autologous cellular transplantation, and therefore, extensive expansion may be required for therapeutic utility (Tamama et al., 2006). Meanwhile, there is a lack of homogenous criteria for isolating and cultivating MSCs and standardized procedures for transplantation. Last but not least, the differentiation fate of transplanted BMSCs is not controllable in vivo.

Adipogenesis, osteogenesis, and chondrogenesis of BMSCs are closely related. The specific differentiation direction is highly regulated by biological, physical, and chemical incentives (Hu et al., 2018). As the most studied terminally differentiated cell of BMSCs, the osteocyte is a mechanosensory commander maintaining homeostasis of the bone by osteoclasts and osteoblasts (Kim and Adachi, 2021). An imbalance resulting BMSCs toward a higher rate of adipogenesis is detrimental to bone health and associates with loss of bone mass and development of musculoskeletal diseases. Bone marrow adipocytes (BMAs) account for 50–70% of the marrow volume in adults (Piotrowska and Tarnowski, 2021). Besides acting as an energy reservoir, this kind of adipose tissue has been proved to possess paracrine, endocrine, and immunomodulatory functions (Sulston and Cawthorn, 2016). It, therefore, is critical for the homeostasis of the bone environment. Chondrocytes can be found in the articular surface of adults, which has poor regenerative capacity. Autologous chondrocyte transplantation is appealing for treating osteoarthritis, but the requirement for biopsies of chondrocytes from a healthy area of the cartilage cap is an evident disadvantage (Boeuf and Richter, 2010). Multiple key genes, non-coding RNAs (miRNA, lncRNA, and circRNA), and signaling pathways deciding the fate of BMSCs have been proposed to act at different time points during differentiation. For instance, the early stage of adipogenic differentiation is marked by augmentation of AMP cyclic production and phosphorylation of CREB (Chen et al., 2016). While during early osteogenesis, there is an upregulation of hedgehog proteins, Wnt/β-catenin signaling, BMPs, endocrine hormones, growth factors, and transcription factors (TFs) represented by RUNX2 (Pierce et al., 2019). SOX9 is pivotal for chondrogenesis by securing chondrocyte lineage commitment, promoting cell survival, and transcriptionally activating the genes for many cartilage-specific structural components (Lefebvre and Dvir Ginzberg, 2017).

Lineage commitment refers to the process of differentiation in which BMSCs lose their multipotency. During this phase, changes consistent with their respective lineage occur in the transcriptional profile, cellular metabolism, and morphology of BMSCs (van de Peppel et al., 2017). A primary question about the differentiation process is at what moment during the differentiation course MSCs become committed to a specific phenotype. Robert and his colleagues suggested that key genes related to osteogenesis were not significantly modulated in the first 24 h of osteogenesis of MSCs, and conversely, in MSCs treated with the adipogenic medium for 24 h, TFs related to adipogenesis and genes associated with lipid metabolism were already differentially expressed (Robert et al., 2020). van de Pepple et al. (2017) divided the early stage (4 d) of osteogenic and adipogenic lineage commitment of human BMSCs into three distinct phases: initiation of differentiation (0–3 h), lineage acquisition (6–24 h), and early lineage progression (48–96 h). The first phase is characterized by expression changes of many TFs and during the second phase, more than 50% of regulated genes are direct targets of the crucial TFs identified in the first phase. In the third phase, the number of differentially expressed genes (DEGs) is stable, indicating that differentiating cells have reached a stable phenotype. However, the author only applied BMSCs of human sources and key genes during chondrogenesis which are not investigated.

Many studies suggested that BMSCs from four species (human, pig, rat, and guinea pig) functioned across species barriers and act as promising sources for xenograft (Li et al., 2012). Currently, a change at the transcriptional level of rat BMSCs during the lineage-commitment phase of differentiation, especially chondrogenesis, remains unclear. In this study, using high-throughput sequencing data and multiple bioinformatics methods, we attempted to reveal key genes deciding osteogenic, adipogenic, and chondrogenic differentiation of rat BMSCs in early phases and analyze the enriched functions and their correlations. Meanwhile, potentially related microRNAs (miRNAs) and signaling pathways were predicted to elucidate underlying regulatory networks. These findings may help to characterize regulatory programs controlling the early stages of lineage commitment of rat BMSCs, to develop novel therapies for age- and pathology-related musculoskeletal diseases and to protect transplanted xenogeneic MSCs from immune detection and enable their long-term survival.

Materials and Methods

This study was approved by the Animal Ethics Committee of Shandong Provincial Hospital affiliated to Shandong First Medical University. The experimental design is depicted in Figure 1A.

FIGURE 1

FIGURE 1. Flow chart of osteogenic, adipogenic, and chondrogenic differentiation of rat BMSCs, data gathering and analysis (A), and histochemical staining of calcium (B) with Alizarin red or adipocyte (C) with oil red O or cartilage (D) with Alcian blue staining after 21 days of osteogenic or adipogenic or chondrogenic differentiation. Scale bars, 50 μm.

BMSC Isolation and Induced Differentiation

BMSCs were extracted from femurs and tibia of three female Sprague–Dawley (SD) rats about 3 weeks old by washing the marrow cavity using the high glucose Dulbecco’s Modified Eagle medium (DMEM) (Gibco, Rockville, MD, United States). The cells were cultured in the DMEM supplemented with 10% fetal bovine serum (FBS) and placed in an incubator with a temperature of 37°C, 5% CO₂ concentration, and 95% relative humidity. BMSCs at passage 2 were seeded in 12-well cell culture plates. After 2 days of incubation (query), these cells were induced to give rise to osteoblasts by using the osteogenic differentiation medium (Promocell, Germany, C-28013), or to adipocytes using the adipogenic differentiation medium (Promocell, Germany, C-28016), or to chondrocytes using the chondrogenic differentiation medium (Promocell, C-28012). Differentiation media were replaced every 3–4 days.

BMSCs Staining

For confirming the result of induced differentiation, Alizarin red staining, Oil Red O (ORO) staining, and Alcian blue staining assays were conducted, as previously described (Bruedigam et al., 2011; Gari et al., 2016).

RNA Isolation

To obtain enough RNA for the gene expression profiling analyses, we pooled three individual cultures in TRIzol (Life Technologies). A total of three experimental samples per time point (3 time points, namely, 3, 12, and 72 h) per lineage (osteogenic, adipogenic, or chondrogenic) were used for gene expression profiling analyses. In addition to the control samples, there were a total of 30 samples.

Next-Generation High-Throughput Sequencing and Identification of DEGs

The RNA quality assessment, mRNA library construction, and mRNA sequencing were performed using methods described in the previous study (Zhang et al., 2021). Based on the annotation information in the platform, the probe sets were transformed into the corresponding gene symbol. The mean value was calculated if multiple probes corresponded to the same gene symbol. The data were normalized using quantile normalization with Illumina in R software (Version 3.6.2). After data normalization, the DEGs were identified using the limma package (Version 3.42.2) in R software. Individual p-values were converted to adjusted p-values (adj. p. val) by false discovery rate correction of the Benjamini and Hochberg test. The threshold of |log₂fold change (FC)| ≥ 1 with an adj. p. val <0.05 was used as the threshold for selecting DEGs. The volcano plot of the DEGs was drawn using the R ggplot2 package.

Functional and Pathway Enrichment Analyses of DEGs

The DEGs were uploaded to an online biological information database, the Database for Annotation, Visualization, and Integrated Discovery (DAVID) version 6.8 Beta (https://david-d.ncifcrf.gov/) for gene ontology (GO) enrichment (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (https://www.genome.jp/kegg/pathway) analysis, which was visualized in the R ggplot2 package. The significance criterion was set at p. val <0.05.

Protein–Protein Interaction Network Construction and Module Analysis

The protein–protein interaction (PPI) network of identified genes was constructed using the multiple protein online tool in the STRING database (Liu et al., 2021), a Search Tool for the Retrieval of Interacting Genes database (version 11.0b, http://string-db.org), and visualized by Cytoscape software (Version 3.6.2). The most significant modules in the PPI network were visualized using the molecular complex detection (MCODE) tool in Cytoscape, and the hub genes were selected by the cytoHubba plugin. The selection criteria were as follows: MCODE scores >5, node score cutoff = 0.2, degree cutoff = 2, max depth = 100, and k-score = 2.

Key Gene Selection

The hub genes (top 10) in the PPI network were identified by using three methods, including maximal clique centrality (MCC), degree, and maximum neighborhood component (MNC) from the cytoHubba plugin in Cytoscape. The biological process analysis of the top 10 genes was performed and visualized using the Biological Networks Gene Ontology tool (BiNGO) (version 3.0.3) plugin in Cytoscape.

DEG–miRNA Network Construction

After hub genes were selected, the DEG–miRNA pairs were identified using the miRDB database (http://mirdb.org/) and visualized in Cytoscape.

Results

Confirmation of the Three-Directional Differentiation of Rat BMSCs

On day 21 after induced differentiation, Alizarin red staining for calcium in the osteogenic differentiation group showed mineralization of the extracellular matrix (Figure 1B); meanwhile, intracellular lipid vesicles accumulated in the adipogenic group (Figure 1C), and cartilage formed in the chondrocyte differentiation group (Figure 1D). Taken together, these results confirmed that rat BMSCs differentiated in accordance with their expected multi-lineage potential.

Data Normalization

The box plot visualization of the gene expression of all samples (Figure 2A) was used to evaluate the data normalization and cross-comparability, which demonstrated that the black lines were almost in the same position after data normalization. Principal component analysis (PCA) confirmed the biological variability between different samples, indicating globally distinct expression profiles (Figure 2B).

FIGURE 2

FIGURE 2. (A) Distribution of expression of 27 samples before and after normalization. (B) Distribution of expression of 27 samples involving principal component analysis (PCA) for confirming biological variability between different samples. (C) Dynamic transcriptional changes upon osteogenic, adipogenic, and chondrogenic differentiation of rat BMSCs with the number of significant differential genes relative to the previous time point (based on three independent experiments) and the number of differentially expressed genes per hour compared to the previous time point (based on three independent experiments). AD, adipogenesis; OS, osteogenesis; CH, chondrogenesis.

Identification of DEGs at Three Time Points

Analysis of gene expression dynamics during rat BMSC differentiation revealed that transcript levels changed significantly at 3 h upon induction of differentiation (Figure 2C). During induction of osteogenic differentiation, 792 genes (398 upregulated and 394 downregulated ones) were significantly different at 3 h, and this increased further to 1,042 (480 upregulated and 562 downregulated ones) at 12 h. The number of significantly modulated genes decreased to 638 (314 upregulated and 324 downregulated ones) at 72 h. During adipogenic differentiation, the number of DEGs was 969 at 3 h, 577 at 12 h, and 680 at 72 h. During chondrogenic differentiation, the number of DEGs was 819 at 3h, 678 at 12 h, and 832 at 72 h. Remarkably, at 3 h after differentiation, the number of upregulated DEGs was higher than the downregulated ones. Next, we calculated the number of differentially expressed DEGs at each time point compared with the preceding time point, and in three directions of differentiation, the number of DEGs that changed per hour decreased as time extended (Figure 2C).

The volcano plot of DEGs at 3, 12, and 72 h of osteogenesis and the comparison results between every two lineages at 3 h after inductive differentiation are presented in Figure 3. Supplementary Figure S1 describes the visualized results of DEGs at three time points of adipogenic and chondrogenic differentiation.

FIGURE 3

FIGURE 3. Volcano plot of differentially regulated genes during osteogenesis at 3 h, (A) 12 h, and (B) 72 h (C) Compared with undifferentiated cells (t = 0 h); Volcano plot of differentially regulated genes between every two lineages at 3 h after inductive differentiation (t = 0 h). (D) Comparison between osteogenesis and adipogenesis, (E) Comparison between osteogenesis and chondrogenesis, and (F) Comparison between chondrogenesis and adipogenesis. AD, adipogenesis; OS, osteogenesis; CH, chondrogenesis.

Identification of Specific DEGs at Three Time Points

In osteogenesis of BMSCs, 203, 381, and 287 DEGs were identified to be differentially expressed only at 3, 12, and 72 h after inductive differentiation (Figure 4A). In adipogenesis of BMSCs, the numbers were 548, 130, and 295 (Figure 4B), and for chondrogenesis of BMSCs, the numbers were 320, 138, and 371 (Figure 4C). Notably, 195 (e.g., Mex3b and Sertad1), 168(e.g., Dtx4 and Il6), and 229 (e.g., Glra1 and Fdps) genes showed differentially expressed levels at all three time points in osteogenesis, adipogenesis, and chondrogenesis, respectively.

FIGURE 4

FIGURE 4. Venn Diagrams of differentially expressed genes during osteogenesis (A), adipogenesis (B), and chondrogenesis (C) at three time points compared with undifferentiated cells (t = 0 h). Venn Diagrams of differentially expressed genes at 3 h (D), 12 h (E), and 72 h (F) of three directions of differentiation compared with undifferentiated cells (t = 0 h). AD, adipogenesis; OS, osteogenesis; CH, chondrogenesis.

Comparison of DEGs Between Different Lineages

At 3 h after induction of differentiation, it was noticed that there were 178, 350, and 220 DEGs only differentially expressed in osteogenesis, adipogenesis, and chondrogenesis, respectively (Figure 4D). At 12 h after inductive differentiation, the numbers were 475, 136, and 252 (Figure 4E), and at 72 h, the numbers were 198, 244, and 448 (Figure 4F). Notably, 382 (e.g., Hmg20a and Mmp24), 184 (e.g., Galnt2 and Cxcl6), and 174 (e.g., Ccnd2 and Sgcd) genes showed differential expression at all three directions of differentiation at 3, 12, and 72 h after inductive differentiation, respectively.

Functional Enrichment Analysis of DEGs

GO and KEGG pathway enrichment analyses of differentially expressed DEGs at different time points of three directions of differentiation were performed to identify the most relevant biological processes (BPs), molecular functions (MFs), cellular components (CCs), and pathways. The top enriched terms of DEGs in BP, CC, MF, and KEGG at 3 h after osteogenesis are presented in Figure 5A and Supplementary Figure S2A. The DEGs at 3 h of osteogenesis were significantly enriched in the “PI3K-Akt signaling pathway” and BP terms such as “negative regulation of transcription from RNA polymerase II promoter” and “response to drug.” The enriched terms at 12 and 72 h after osteogenesis are presented in Figures 5B,C and Supplementary Figures S2B,C, while results of the other two directions of differentiation are presented in Supplementary Figures S3, S4.

FIGURE 5

FIGURE 5. GO including the biological process and KEGG enrichment analysis of differentially expressed genes at 3 h, (A) 12 h, and (B) 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h).

The DEGs between osteogenesis and adipogenesis at 3 h after inductive differentiation were also significantly enriched in BP terms such as “positive regulation of cell migration” and “negative regulation of transcription from RNA polymerase II promoter.” Results of functional enrichment analysis of the DEGs between every two lineage commitment at 3 h of inductive differentiation are presented in Figure 6 and Supplementary Figure S5.

FIGURE 6

FIGURE 6. GO including the biological process and KEGG enrichment analysis of differentially expressed genes between every two lineages at 3 h after inductive differentiation (t = 0 h). (A) Comparison between osteogenesis and adipogenesis, (B) Comparison between osteogenesis and chondrogenesis, and (C) Comparison between chondrogenesis and adipogenesis.

PPI Network Construction and Module Analysis

The most significant module of DEGs at 3 h of osteogenesis consisted of 65 nodes and 907 edges, with the Rrp15 gene being identified to have a relatively high connectivity degree (Figure 7A). The most significant modules at 12 and 72 h after osteogenesis are presented in Figures 7B,C, respectively. Meanwhile, the most significant modules of the other two-directional differentiation at three time points are presented in Supplementary Figure S6. The most significant module about interactions of DEGs between osteogenesis and adipogenesis at 3 h after inductive differentiation consisted of 7 nodes and 14 edges, with the Vcl gene being identified to have a relatively high connectivity degree (Figure 7D). The most significant modules regarding interactions of DEGs between osteogenesis and chondrogenesis and between adipogenesis and chondrogenesis at 3 h after inductive differentiation are presented in Figures 7E,F. Hub genes were selected by CytoHubba. The top 10 hub genes, which were selected based on the three most commonly used classification methods in cytoHubba, are presented in Figure 8 and Supplementary Figure S7. The top 10 hub genes selected by degree in different comparisons are presented in Table 1. The fold change (FC) and adj. p. val of the top 10 hub genes at 3, 12, and 72 h of osteogenesis compared with undifferentiated cells are presented in Table 2. The biological process analysis of the top 10 hub genes was performed and visualized using BiNGO in Cytoscape, as shown in Figure 9 and Supplementary Figures S8, S9.

FIGURE 7

FIGURE 7. Most significant module in protein–protein interaction (PPI) network analysis of differentially expressed genes at 3 h, (A) 12 h, and (B), and 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h); The most significant module of differentially regulated genes between osteogenesis and adipogenesis, (D) between osteogenesis and chondrogenesis, and (E) between chondrogenesis and adipogenesis (F) at 3 h after inductive differentiation (t = 0 h). Red indicates upregulated genes, and green indicates the downregulated genes.

FIGURE 8

FIGURE 8. Top 10 key genes at 3 h, (A) 12 h, and (B) 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h); The top 10 key genes identified from comparisons between osteogenesis and adipogenesis, (D) between osteogenesis and chondrogenesis, and (E) between chondrogenesis and adipogenesis (F) at 3 h after inductive differentiation (t = 0 h), which were selected based on the three most commonly used classification methods (MCC, MNC, and Degree) in cytoHubba.

TABLE 1

TABLE 1. List of the top 10 hub genes selected by degree methods in cytoHubba.

TABLE 2

TABLE 2. The fold change (FC) and adj. p. val of the top 10 hub genes at 3, 12, and 72 h of osteogenesis compared with undifferentiated cells.

FIGURE 9

FIGURE 9. Biological process analysis of the top 10 differentially expressed genes was constructed using BiNGO at 3 h, (A) 12 h, and (B) 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h). The color depth of nodes refers to the corrected p-value of ontologies. The size of nodes refers to the number of genes that are involved in the ontologies. p < 0.001 was considered statistically significant.

Construction of DEG-miRNA

The miRNA-DEG pairs were identified through network analysis of 10 DEGs using the miRDB databases. At 3 h of osteogenesis, a total of 32 associations between 22 miRNAs and 7 DEGs were identified, and then, the network was visualized in Cytoscape (Figure 10A). As is shown, Rrm1 interacts with four miRNAs (rno-miR-320-3p, rno-miR-145-3p, rno-miR-330-3p, and rno-miR-291a-5p) and four genes (Top2a, Ccna2, Mcm4, and Kif11). The DEG-miRNA networks regarding osteogenesis at 12 and 72 h are presented in Figures 10B,C. The DEG-miRNA networks about comparisons between every two lineage commitment at 3 h after inductive differentiation are presented in Figures 10D–F. The DEG-miRNA networks regarding adipogenesis and chondrogenesis at 3, 12, and 72 h after differentiation are presented in Supplementary Figure S10.

FIGURE 10

FIGURE 10. Regulatory network of DEG-miRNA was obtained from the miRDB database using the DEGs at 3 h, (A) 12 h, and (B) 72 h (C) in osteogenesis compared with undifferentiated cells (t = 0 h), and the DEGs identified from comparisons between osteogenesis and adipogenesis, (D) between osteogenesis and chondrogenesis, and (E) between chondrogenesis and adipogenesis (F) at 3 h after inductive differentiation (t = 0 h). Blue indicates genes, and green indicates targeted miRNA. DEG, differentially expressed gene; miRNA, microRNA.

Discussion

There were no previous studies paying special attention to the lineage commitment of rat BMSCs (Arthur and Gronthos, 2020; Zhu et al., 2020). In the current investigation, three time points (3, 12, and 72 h) after inductive differentiation were chosen to represent the three phases of early lineage commitment of rat BMSCs. Next-generation high-throughput sequencing was adapted to generate potential key genes. At the first time point (3 h), gene expression analyses revealed many DEGs in all three directions of differentiation. The first phase of early lineage commitment represents the initiation stage of the differentiation program. It is characterized by expression changes of many transcription factors that regulate the lineage commitment and set the multipotential stem cells toward a stable phenotype (van de Peppel et al., 2017). As for the second phase, differentiating lineages begin to deviate. In this study, most genes related to osteogenesis were modulated later than those related to adipogenesis, which confirmed the finding of the study conducted by Robert and his colleagues (Robert et al., 2020). In the third stage, the differentiating cells reach a stable phenotype (van de Peppel et al., 2017).

Integrated analysis of differential gene expression at 3 h of osteogenesis indicated that the top 10 DEGs (Ccna2, Cdc20, Il6, Mad2l1, Top2a, Kif11, Mcm4, Asf1b, Rrm1, and Pa2g4) in the PPI network were the hub genes; only one gene (Il6) showed a decreasing trend compared to the control group. Similarly, at 12 h of osteogenesis, nine (Mad2l1, Top2a, Ccna2, Cdc20, Mcm6, Mcm3, Kif11, Mcm4, and Cdc6) of the identified top 10 hub genes in the PPI network were upregulated except Il6. This trend was reversed at 72 h of osteogenesis as it has been shown that at this moment, nine (Il6, Cd44, Il1b, Ccl2, Kif11, Asf1b, Cdc20, Mcm3, and Pbk) of the top 10 DEGs were downregulated except Notch1. For adipogenesis, eight (Ccna2, Cdc20, Kif11, Top2a, Mad2l1, Mcm6, Asf1b, and Mcm3) hub genes were upregulated and two (Il6 and Gsk3b) were downregulated at 3 h of differentiation, while for chondrogenesis, only two (Il6 and Sspo) hub genes were downregulated at this time point. In fact, there is rather a large overlap regarding hub genes identified at 3 h of osteogenesis, chondrogenesis, and adipogenesis. Three genes (Ccna2, Cdc20, and Il6) were upregulated in all three groups, which suggested that most identified hub genes at this phase may not play independent roles in determining the lineage commitment of BMSCs but act as common participants in initiating differentiation. Ccna2 is ubiquitously expressed and acts as an essential regulator of the cell division cycle at the G1/S transition and mitotic entry (Bendris et al., 2012; Loukil et al., 2014). Cdc20 is a key regulator of ubiquitylation of FRET and could colocalize with Ccna2, LC3-B, p62 or lysosomes (Loukil et al., 2014). Il6 is a well-known pleiotropic cytokine and shows multiple functions such as the regulation of inflammation and immunity (Hunter and Jones, 2015). Deng et al. (2020) have demonstrated that the expression of Il6 increased during adipogenesis of MSCs and was positively correlated with the Oil Red O quantification result.

Unique DEGs that occurred in only one lineage may play pivotal roles in determining the fate of BMSCs. Five genes (LOC689899, Mex3b, Sertad1, LOC100362684, and Hopx) were dysregulated in three time points of osteogenesis but showed non-significant changes during two other lineages (Supplementary Table S1). Mex3B, which is reportedly to induce apoptosis after cell stimulations such as heat, radiation, and reactive oxygen species, is downregulated in osteogenesis (Umehara et al., 2020). Hopx has been demonstrated to be a promoter of BMSC proliferation and an inhibitor of adipogenesis by suppressing adipogenic pathway–related genes such as PLIN1, FABP4, and PLIN5 (Hng et al., 2020). Dtx4 was exclusively upregulated in three time points during adipogenesis and was involved in cell development and cell differentiation (Wang et al., 2017). Wang et al. (2017) indicated that knockdown of this gene could inhibit adipogenesis through restraining C/EBPα and PPARγ and therefore arresting mitotic clonal expansion. Sixteen genes were upregulated through the early lineage of chondrogenesis. Among them, the expression of Ibsp is driven by RUNX3 and MEF2C, which appears more critical than RUNX2 for chondrogenesis (Dreher et al., 2020). It could be noticed that a lot of genes are exclusively dysregulated at the first phase of differentiation, and their effect for determining the fate of rat BMSCs needs to be further confirmed (Supplementary Table S1). A total of 207 unique DEGs were only identified at 3 h of adipogenesis, which suggested that key genes related to adipogenesis were significantly modulated earlier than osteogenesis.

Then, various bioinformatics analysis methods were used to reveal the enriched functional terms regarding identified DEGs. As a result, GO terms of 792 DEGs identified at 3 h of osteogenesis for biological process categories included “negative regulation of transcription from RNA polymerase II promoter” and “response to drug.” Furthermore, the KEGG pathway enrichment analysis revealed that these DEGs were predominantly associated with the “PI3K-Akt signaling pathway.” The DEGs in the “PI3K-Akt signaling pathway” were Chrm2, Pdgfra, Ngfr, Col24a1, Lama1, Lama4, Flt4, Pik3r5, Rbl2, Il6, Ppp2r1b, Fgf9, Col5a1, Col5a2, Col4a3, Eif4ebp1, Tek, Ywhah, Ifnar1, Pck2, and Epha2. As to 12 h of osteogenesis, the most involved BP terms are “oxidation-reduction process” and “negative regulation of cell proliferation,” which suggest that BMSCs start from proliferation to differentiation at this phase. While in this phase, the most enriched KEGG pathway is “metabolic pathways.” Similar changes occurred early in adipogenesis in comparison with osteogenesis. It can be noticed that at 3 h of adipogenesis, the most enriched BP terms are “oxidation-reduction process,” “positive regulation of cell proliferation,” “negative regulation of apoptotic process,” and “response to drug,” while the most involved KEGG pathways are “metabolic pathways.”

miRNAs have been demonstrated to be important regulatory factors in the differentiation of BMSCs by binding the 3′-untranslated region of targeted mRNAs (Zhou et al., 2020; Liu et al., 2021). Thus, we constructed a DEG-miRNA regulatory network to show the potential interactions among identified DEGs and RNAs during osteogenesis, adipogenesis, and chondrogenesis. In the networks of DEGs about comparisons between adipogenesis and osteogenesis at 3h, rno-miR-449-5p may play an important role because it connects with three key genes. This miRNA, along with rno-miR-153-5p, also occurred in the network about comparisons between osteogenesis and chondrogenesis. Rno-miR449a-5p are essential for both mitochondria metabolic dysfunction and phenotype transformation of pulmonary arterial smooth muscle cells (Zhang et al., 2019), but currently, there is no direct evidence to prove the role of this miRNA for the lineage commitment of rat BMSCs.

For the first time, the gene transcription of rat BMSCs during osteogenesis, adipogenesis, and chondrogenesis was systematically detected and analyzed. The current study provided directions and molecular targets for further investigations on the lineage commitment of rat BMSCs by high-throughput sequencing and rigorous bioinformatics analysis. However, it should be admitted that there were shortcomings that merit consideration. First, this is only a very fundamental research study using two-dimensional culture. The differentiation in the internal environment may possess different characteristics in comparison with the in vitro experiment. These results cannot completely represent the normal physiological processes. The role of predicted genes and miRNAs in the three-directional differentiation of BMSCs merits further investigation. Second, only three time points distributed so-called “three phases of lineage commitment” were selected in this study as representatives. More detailed results may be acquired if more time points during differentiation were selected. Third, the DEG-miRNA network was constructed by using the miRDB database. In the future, we will detect the expression of non-coding RNA using our samples and then construct a competing endogenous RNA network. Last but not least, because there is a large number of differentially expressed genes identified in the diverse comparisons, validation by Western blot assay or protein array was not available. In the future, specific molecules could be selected for in-depth investigation.

Conclusion

These hub genes, functional terms, and predicted targeted genes identified in the current study may provide information for further clarification of the molecular mechanism underlying the lineage commitment of rat bone mesenchymal stem cells. Hub genes including Ccna2, Cdc20, and Il6 may act as typical participants in initiating osteogenesis, adipogenesis, and chondrogenesis. Mex3b, Sertad1, and Hopx showed enhanced expression throughout three early phases during the osteogenic differentiation but no significant change in other two-directional differentiation. A similar pattern of Dtx4 and Ibsp expressions occurred in adipogenesis and chondrogenesis, respectively. Additionally, rno-miR-449-5p is a potential regulator for the early lineage commitment of osteogenesis. Our findings will help understand the underlying mechanism determining the differentiation fate of BMSCs and provide theoretical support for the clinical treatment of osteoporosis, osteoarthritis, and other age-related bone diseases.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, accession ID: GSE185140.

Ethics Statement

This study was approved by the Animal Ethics Committee of Shandong Provincial Hospital affiliated to Shandong First Medical University

Author Contributions

FL and QZ conceived and designed the study. FL, QZ, JD, and DZ analyzed RNA sequencing data. FL, QZ, and PZ performed the basic experiment. FL, and QZ wrote the manuscript, which was commented on by all authors.

Funding

This work was partially supported by the China Scholarship Council (CSC) (Grant number 201808080126), the Incubation Fund of Shandong Provincial Hospital (Grant number 2020FY019), the Young Scholars Program of Shandong Provincial Hospital, the Young Taishan Scholars Program of Shandong Province (Grant number tsqn201909183), the Academic Promotion Program of Shandong First Medical University (Grant number 2020RC008), and the Natural Science Foundation of Shandong Province (Grant numbers ZR2020QH072 and ZR2021QH307).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.782054/full#supplementary-material

References

Arthur, A., and Gronthos, S. (2020). Clinical Application of Bone Marrow Mesenchymal Stem/Stromal Cells to Repair Skeletal Tissue. Ijms 21 (24), 9759. doi:10.3390/ijms21249759

CrossRef Full Text | Google Scholar

Bendris, N., Loukil, A., Cheung, C., Arsic, N., Rebouissou, C., Hipskind, R., et al. (2012). Cyclin A2: A Genuine Cell Cycle Regulator. Biomol. Concepts 3 (6), 535–543. doi:10.1515/bmc-2012-0027

PubMed Abstract | CrossRef Full Text | Google Scholar

Boeuf, S., and Richter, W. (2010). Chondrogenesis of Mesenchymal Stem Cells: Role of Tissue Source and Inducing Factors. Stem Cel Res. Ther. 1 (4), 31. doi:10.1186/scrt31

CrossRef Full Text | Google Scholar

Bruedigam, C., Driel, M. v., Koedam, M., Peppel, J. v. d., van der Eerden, B. C. J., Eijken, M., et al. (2011). Basic Techniques in Human Mesenchymal Stem Cell Cultures: Differentiation into Osteogenic and Adipogenic Lineages, Genetic Perturbations, and Phenotypic Analyses. Curr. Protoc. Stem Cel Biol. 17, t1H–t3H. doi:10.1002/9780470151808.sc01h03s17

CrossRef Full Text | Google Scholar

Chen, Q., Shou, P., Zheng, C., Jiang, M., Cao, G., Yang, Q., et al. (2016). Fate Decision of Mesenchymal Stem Cells: Adipocytes or Osteoblasts. Cell Death Differ 23 (7), 1128–1139. doi:10.1038/cdd.2015.168

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, W., Chen, H., Su, H., Wu, X., Xie, Z., Wu, Y., et al. (2020). IL6 Receptor Facilitates Adipogenesis Differentiation of Human Mesenchymal Stem Cells through Activating P38 Pathway. Ijsc 13 (1), 142–150. doi:10.15283/ijsc19073

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreher, S. I., Fischer, J., Walker, T., Diederichs, S., and Richter, W. (2020). Significance of MEF2C and RUNX3 Regulation for Endochondral Differentiation of Human Mesenchymal Progenitor Cells. Front. Cel Dev. Biol. 8, 81. doi:10.3389/fcell.2020.00081

PubMed Abstract | CrossRef Full Text | Google Scholar

Galderisi, U., Peluso, G., and Di Bernardo, G. (2021). Clinical Trials Based on Mesenchymal Stromal Cells Are Exponentially Increasing: Where Are We in Recent Years. Stem Cel Rev. Rep. 16. doi:10.1007/s12015-021-10231-w

CrossRef Full Text | Google Scholar

Gari, M., Alsehli, H., Gari, A., Abbas, M., Alkaff, M., Abuzinadah, M., et al. (2016). Derivation and Differentiation of Bone Marrow Mesenchymal Stem Cells from Osteoarthritis Patients. Tissue Eng. Regen. Med. 13 (6), 732–739. doi:10.1007/s13770-016-0013-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hng, C. H., Camp, E., Anderson, P., Breen, J., Zannettino, A., and Gronthos, S. (2020). HOPX Regulates Bone Marrow-Derived Mesenchymal Stromal Cell Fate Determination via Suppression of Adipogenic Gene Pathways. Sci. Rep. 10 (1), 11345. doi:10.1038/s41598-020-68261-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, L., Yin, C., Zhao, F., Ali, A., Ma, J., and Qian, A. (2018). Mesenchymal Stem Cells: Cell Fate Decision to Osteoblast or Adipocyte and Application in Osteoporosis Treatment. Ijms 19 (2), 360. doi:10.3390/ijms19020360

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, C. A., and Jones, S. A. (2015). IL-6 as a keystone Cytokine in Health and Disease. Nat. Immunol. 16 (5), 448–457. doi:10.1038/ni.3153

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Zhang, P., Zhang, X., Lv, L., and Zhou, Y. (2021). Advances in Mesenchymal Stem Cell Transplantation for the Treatment of Osteoporosis. Cell Prolif 54 (1), e12956. doi:10.1111/cpr.12956

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, J., Shojaee, M., Mitchell Crow, J., and Khanabdali, R. (2021). From Mesenchymal Stromal Cells to Engineered Extracellular Vesicles: A New Therapeutic Paradigm. Front. Cel Dev. Biol. 9, 705676. doi:10.3389/fcell.2021.705676

CrossRef Full Text | Google Scholar

Kim, J., and Adachi, T. (2021). Cell-fate Decision of Mesenchymal Stem Cells toward Osteocyte Differentiation Is Committed by Spheroid Culture. Sci. Rep. 11 (1), 13204. doi:10.1038/s41598-021-92607-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Lefebvre, V., and Dvir-Ginzberg, M. (2017). SOX9 and the many Facets of its Regulation in the Chondrocyte Lineage. Connect. Tissue Res. 58 (1), 2–14. doi:10.1080/03008207.2016.1183667

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Ezzelarab, M. B., and Cooper, D. K. C. (2012). Do mesenchymal Stem Cells Function across Species Barriers? Relevance for Xenotransplantation. Xenotransplantation 19 (5), 273–285. doi:10.1111/xen.12000

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Dong, J., Zhou, D., and Zhang, Q. (2021). Identification of Key Candidate Genes Related to Inflammatory Osteolysis Associated with Vitamin E-Blended UHMWPE Debris of Orthopedic Implants by Integrated Bioinformatics Analysis and Experimental Confirmation. Jir Vol. 14, 3537–3554. doi:10.2147/JIR.S320839

CrossRef Full Text | Google Scholar

Loukil, A., Zonca, M., Rebouissou, C., Baldin, V., Coux, O., Biard-Piechaczyk, M., et al. (2014). High Resolution Live Cell Imaging Reveals Novel Cyclin A2 Degradation Foci Involving Autophagy. J. Cel Sci. 127 (Pt 10), 2145–2150. doi:10.1242/jcs.139188

CrossRef Full Text | Google Scholar

Mabuchi, Y., Okawara, C., Méndez-Ferrer, S., and Akazawa, C. (2021). Cellular Heterogeneity of Mesenchymal Stem/Stromal Cells in the Bone Marrow. Front. Cel Dev. Biol. 9, 689366. doi:10.3389/fcell.2021.689366

CrossRef Full Text | Google Scholar

Pierce, J. L., Begun, D. L., Westendorf, J. J., and McGee-Lawrence, M. E. (2019). Defining Osteoblast and Adipocyte Lineages in the Bone Marrow. Bone 118, 2–7. doi:10.1016/j.bone.2018.05.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Piotrowska, K., and Tarnowski, M. (2021). Bone Marrow Adipocytes-Role in Physiology and Various Nutritional Conditions in Human and Animal Models. Nutrients 13 (5), 1412. doi:10.3390/nu13051412

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, A. W., Marcon, B. H., Dallagiovanna, B., and Shigunov, P. (2020). Adipogenesis, Osteogenesis, and Chondrogenesis of Human Mesenchymal Stem/Stromal Cells: A Comparative Transcriptome Approach. Front. Cel Dev. Biol. 8, 561. doi:10.3389/fcell.2020.00561

PubMed Abstract | CrossRef Full Text | Google Scholar

Saeedi, M., Nezhad, M. S., Mehranfar, F., Golpour, M., Esakandari, M. A., Rashmeie, Z., et al. (2021). Biological Aspects and Clinical Applications of Mesenchymal Stem Cells: Key Features You Need to Be Aware of. Cpb 22 (2), 200–215. doi:10.2174/1389201021666200907121530

CrossRef Full Text | Google Scholar

Sulston, R. J., and Cawthorn, W. P. (2016). Bone Marrow Adipose Tissue as an Endocrine Organ: Close to the Bone. Horm. Mol. Biol. Clin. Investig. 28 (1), 21–38. doi:10.1515/hmbci-2016-0012

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamama, K., Fan, V. H., Griffith, L. G., Blair, H. C., and Wells, A. (2006). Epidermal Growth Factor as a Candidate for Ex Vivo Expansion of Bone Marrow-Derived Mesenchymal Stem Cells. Stem Cells 24 (3), 686–695. doi:10.1634/stemcells.2005-0176

PubMed Abstract | CrossRef Full Text | Google Scholar

Umehara, T., Kagawa, S., Tomida, A., Murase, T., Abe, Y., Shingu, K., et al. (2020). Body Temperature-dependent microRNA Expression Analysis in Rats: Rno-miR-374-5p Regulates Apoptosis in Skeletal Muscle Cells via Mex3B under Hypothermia. Sci. Rep. 10 (1), 15432. doi:10.1038/s41598-020-71931-w

PubMed Abstract | CrossRef Full Text | Google Scholar

van de Peppel, J., Strini, T., Tilburg, J., Westerhoff, H., van Wijnen, A. J., and van Leeuwen, J. P. (2017). Identification of Three Early Phases of Cell-Fate Determination during Osteogenic and Adipogenic Differentiation by Transcription Factor Dynamics. Stem Cel Rep. 8 (4), 947–960. doi:10.1016/j.stemcr.2017.02.018

CrossRef Full Text | Google Scholar

Wang, Z., Dai, Z., Pan, Y., Wu, S., Li, Z., and Zuo, C. (2017). E3 Ubiquitin Ligase DTX4 Is Required for Adipogenic Differentiation in 3T3-L1 Preadipocytes Cell Line. Biochem. Biophysical Res. Commun. 492 (3), 419–424. doi:10.1016/j.bbrc.2017.08.083

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Ma, C., Zhang, L., Zhang, L., Zhang, F., Ma, M., et al. (2019). MiR-449a-5p Mediates Mitochondrial Dysfunction and Phenotypic Transition by Targeting Myc in Pulmonary Arterial Smooth Muscle Cells. J. Mol. Med. 97 (3), 409–422. doi:10.1007/s00109-019-01751-7

CrossRef Full Text | Google Scholar

Zhang, Q., Dong, J., Zhang, P., Zhou, D., and Liu, F. (2021). Dynamics of Transcription Factors in Three Early Phases of Osteogenic, Adipogenic, and Chondrogenic Differentiation Determining the Fate of Bone Marrow Mesenchymal Stem Cells in Rats. Front. Cel Dev. Biol. 9, 768316. doi:10.3389/fcell.2021.768316

CrossRef Full Text | Google Scholar

Zhou, L., Qiu, M., Yang, L., Yang, L., Zhang, Y., Mu, S., et al. (2020). MicroRNA-1-3p Enhances Osteoblast Differentiation of MC3T3-E1 Cells by Interacting with Hypoxia-Inducible Factor 1 α Inhibitor (HIF1AN). Mech. Develop. 162, 103613. doi:10.1016/j.mod.2020.103613

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, B., Xue, F., Li, G., and Zhang, C. (2020). CRYAB Promotes Osteogenic Differentiation of Human Bone Marrow Stem Cells via Stabilizing β‐catenin and Promoting the Wnt Signalling. Cel Prolif 53 (1), e12709. doi:10.1111/cpr.12709

CrossRef Full Text | Google Scholar