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

The imbalance of osteogenic, adipogenic, and chondrogenic differentiation in bone marrow mesenchymal stem cells (BMSCs) occurred in multiple age-related degenerative diseases such as osteoporosis and osteoarthritis. In order to improve our understanding and control of multi-directional differentiation of BMSCs in rats, using high-throughput sequencing, we identified key gene regulatory events in the early stages of lineage commitment. Data analysis revealed two transcription factors (TFs, Tsc22d3, and Epas1) with elevated expression throughout the initiation of differentiation (3 h), lineage acquisition (12 h), and early lineage progression (72 h) of three-directional differentiation. For osteogenic differentiation, 792, 1,042, and 638 differentially expressed genes including 48, 59, and 34 TFs were identified at three time points, respectively. Moreover, the functional analysis demonstrated that 4, 12, and 5 TFs were only differentially expressed during osteogenic differentiation at 3, 12, and 72 h, respectively, and not during other two-directional differentiation. 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 Gbx2 expression occurred in chondrogenic differentiation. Thus, Hopx and other early responder TFs may control the osteogenic cell fate of BMSCs and participate in the development of osteoporosis. Gbx2 and other early responder TFs should be considered in mechanistic models that clarify cartilage-anabolic changes in the clinical progression of osteoarthritis.

Introduction

Bone mesenchymal stem/stromal cells (BMSCs) have the potentials of differentiating into several mesodermal derivatives, in particular osteogenic (osteogenic progenitors, osteocytes, osteoblasts, and bone lining cells), adipogenic and chondrogenic cells in the existence of physiological demands or pathological changes (Arthur and Gronthos, 2020; Jiang et al., 2021). Systematic cell-based therapy by administration of in vivo expanded BMSCs could serve as an excellent alternative for bone regeneration and tissue engineering (Arthur and Gronthos, 2020; Robert et al., 2020; Jiang et al., 2021), but it is not yet a standard clinical practice due to donor variation among patients and unpredictable direction for cell differentiation (Crapnell et al., 2013; O’Connor, 2019).

The adipogenesis and osteogenesis of BMSCs are competing and reciprocal, keeping homeostasis within the bone. The primary role of BMSCs-derived osteoblasts is to produce and secrete osteoid and mineralizing factors to remodel bone according to the metabolic and structural needs of the body (Yu and Wang, 2016). During aging and the development of bone remodeling disorders such as osteoporosis and osteonecrosis, the normal balance between osteogenic and adipogenic cell populations can shift toward the latter and result in increased formation of marrow fat (Zhang et al., 2018; Pittenger et al., 2019; Robert et al., 2020; Saeedi et al., 2021). Hyaline cartilage is responsible for bone formation in the embryo and can be found in the articular surface of the bones of adults (Carballo et al., 2017). Until now, the in vitro differentiation process is unable to produce cartilage similar to articular cartilage formed under normal physiological conditions (Huynh et al., 2019). Further understanding the mechanism for the three-directional differentiation of BMSCs is imperative.

The differentiation of BMSCs is a two-step process, lineage commitment (from MSCs to lineage-specific progenitors) and maturation (from progenitors to specific cell types) (Robert et al., 2020). The lineage commitment of MSCs to adipocytes, osteoblasts and chondrocytes share a common precursor (Chen et al., 2016; Yianni and Sharpe, 2020), and intensive studies in recent decades have demonstrated that a number of signaling pathways and transcription factors (TFs) are critical for this process (van de Peppel et al., 2017). RUNX2, SP7/Osterix, and SOX9 are among the leading TFs involved in the differentiation of MSCs into osteoblasts or chondrocytes, while PPARγ and cyclic AMP response element-binding protein (CREB) are required for adipogenic differentiation (Chen et al., 2016; Pittenger et al., 2019; Robert et al., 2020). Overexpression of RUNX2 in non-osteoblastic cells and adipose-tissue-derived MSCs elevates the level of osteoblastic markers, promotes mineralization and inhibits lineage commitment to adipocytes, indicating that a single TF is important for bone development and osteoblast differentiation (Robert et al., 2020). During adipogenesis, the production of cyclic AMP is augmented, resulting in phosphorylation of CREB and an increased expression of CEBPβ (Merrett et al., 2020). CEBPβ, together with the other two CEBP family members (CEBPα and CEBPδ), was reported to be implicated in adipogenic differentiation (Chen et al., 2016; Hu et al., 2018; Robert et al., 2020).

Although the enriched biological processes and TFs are mostly identical between these lineages within the first 4 days after induction (van de Peppel et al., 2017; Fang et al., 2020), the early stages of lineage commitment of human BMSCs could be further divided into three distinct phases with own mRNA dynamics: initiation of differentiation (0–3 h, phase I), lineage acquisition (6–24 h, phase II), and early lineage progression (48–96 h, phase III) (van de Peppel et al., 2017). Downstream analyses of four TFs (TP53, FOSl1, ESR1, and HOXA9) identified during the first phase of osteogenic differentiation are capable of regulating more than 50% of the differentially expressed genes (DEGs) in the second and third phases. However, the expression of TFs in chondrogenic differentiation during this stage was not systematically investigated. Meanwhile, this study used human BMSCs, while key TFs playing critical roles in lineage commitment of rat BMSCs were not considered.

BMSCs from 4 different species (human, pig, rat, guinea pig) fail to elicit a proliferative response from allogeneic lymphocytes and therefore can function across the species barrier (Li et al., 2012). Given the unique immunologic properties of rat BMSC, BMSC from rats is an attractive option for xenografts. In the current study, in vivo differentiation of rat BMSCs from a two-dimensional culture system was applied, and by using high-throughput transcriptome sequencing technique and multiple bioinformatics methods, we systematically investigated genes, especially TFs being regulated upon lineage commitment of rat BMSCs into adipocytes, chondrocytes, and osteocytes. One of our key findings was the characterization of genes and regulatory programs controlling the early stages of lineage commitment of rat BMSCs. Progress in this area will be of great value to understand age- and pathology-related fractures and metabolic syndromes, and to guide a better application of BMSCs in tissue engineering and regenerative medicine.

Materials and Methods

Cell Isolation and Induced Differentiation

This study was approved by the Animal Ethics Committee of Shandong Provincial Hospital affiliated to Shandong First Medical University. BMSCs were collected from femurs and tibia of 3 female Sprague Dawley (SD) rats (3 weeks old) by washing the marrow cavity using high glucose Dulbecco’s Modified Eagle Medium (DMEM) (Gibco, Rockville, MD, United States) until bones appear white. Obtained cell clumps were dispersed to cell suspension by using 21-gauge needles. The cells were cultured in DMEM supplemented with 10% fetal bovine serum (FBS) in an incubator of 5% CO₂ and 95% humidity at 37°C. BMSCs at passage 2 were seeded in 12-well cell culture plates and after 2 days, cells were induced to differentiate into osteoblasts using mesenchymal stem cell osteogenic differentiation medium (Promocell, Germany, C-28013), or to adipocytes using mesenchymal stem cell adipogenic differentiation medium (Promocell, Germany, C-28016), or to chondrocytes using mesenchymal stem cell chondrogenic differentiation medium (Promocell, C-28012). Differentiation media were replaced every 3–4 day.

Protein Concentration Assay and Bone Marrow Mesenchymal Stem Cells Staining

Protein concentration assay, Elisa assay, alizarin red staining, Oil Red O staining and Alcian blue staining were performed 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), resulting in a total of three experimental samples per time point (3 time points, namely, 3, 12, and 72 h) per lineage (osteogenic, adipogenic, or chondrogenic) to be used for gene expression profiling analyses. Plus, the control samples, there were a total of 30 samples. RNA was isolated as previously described (Bruedigam et al., 2011).

High-Throughput Sequencing

The quality of RNA was analyzed using NanoDrop^TM One (Thermo Fisher Scientific, United States). An mRNA library was constructed by NEBNext^® Ultra^TM II mRNA Library Prep Kit for Illumina^® (NEB, United States). Qubit^TM dsDNA HS Assay Kit (Invitrogen, United States), D1000 Screen Tape (Agilent, United States), and KAPA Library Quant kit (illumine, United States) universal qPCR Mix were used to assess the quality and yield of the constructed library. Finally, the mRNA was sequenced by NovaSeq 6000 (Illumina, United States).

Sequencing Result Analysis

After standardization and quality control of the sequencing data, the differential expression of mRNA was analyzed. | log2fold change (FC)| in mRNA expression ≥ 1 with an adjusted p-value < 0.05 was used as the threshold for determining gene upregulation or downregulation. Differentially expressed TFs were identified by using the AnimalTFDB 3.0 database.¹ Gene Ontology (GO) enrichment² and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment³ analysis was utilized to determine the biological function of the differentially expressed TFs and investigate its possible involvement in the lineage-commitment of rat BMSCs. The protein-protein interaction (PPI) network of identified genes was predicted using the multiple protein online tool in STRING database (Liu et al., 2021), a search tool for the Retrieval of Interacting Genes database (version 11.0b)⁴ and visualized in Cytoscape software (Version 3.6.2). Relationships between critical TFs and their differentially expressed target genes were predicted to construct the TF-target regulatory network based on the previous method (Zhang et al., 2017; Zhou et al., 2021).

Statistical Analysis

Statistical analyses were conducted by using SPSS Software, version 17.0 (IBM SPSS Statistics, United States). One-way analysis of variance (ANOVA) followed by LSD method was used to compare the protein concentration and ALP expression between different groups, which are expressed as the means ± standard deviation. Differences between means were considered statistically significant when p-values were < 0.05.

Results

Differentiation of Rat Bone Marrow Mesenchymal Stem Cells Into Osteoblasts, Adipocytes, and Chondrocytes

On day 21 after induced differentiation, histological staining for calcium in the osteogenic differentiation group showed that extracellular matrix (ECM) was mineralized (Figure 1A), Oil Red O staining for lipids in the adipogenic differentiation group showed that cells accumulated intracellular lipid vesicle (Figure 1B), and Alcian blue staining in the chondrocyte differentiation group showed cartilage formation (Figure 1C). Biochemical analyses of samples during differentiation revealed that total protein and alkaline phosphatase (ALP) concentration transiently increased during adipogenic, osteogenic, and chondrogenic differentiation (Figures 1D,E). Together, our observations established that rat BMSCs differentiated into three lineages, consistent with their expected multi-lineage potential.

FIGURE 1

Figure 1. Adipogenic, osteogenic and chondrogenic differentiation of rat BMSCs. Histochemical staining of calcium (A) with Alizarin red or adipocyte (B) with oil red O or cartilage (C) with Alcian blue staining after 21 days of osteogenic or adipogenic or chondrogenic differentiation. Scale bars, 100 and 50 μm; total protein (D) and ALP (E) concentration in osteogenesis, adipogenesis and chondrogenesis of rat BMSC. Error bars indicate SD. Mean and SD of three independent experiments with three technical replicates.

Identification of Differentially Expressed Transcription Factors in Three Phases

Data cross−comparability was evaluated via principal component analysis (PCA) for confirming biological variability between different samples. All samples were grouped separately (Figure 2A), indicating globally distinct expression profiles. Analysis of gene expression dynamics during rat BMSC differentiation revealed that transcript levels changed significantly 3 h upon induction of differentiation (Figure 2B). During induction of osteogenic differentiation, 792, 1,042, and 638 DEGs were identified at three time points, respectively. At least 48 TFs (21 upregulated and 27 downregulated ones) were significantly different at 3 h and this increased further to 59 (25 upregulated and 34 downregulated ones) at 12 h. The number of significantly modulated TFs decreased to 34 (22 upregulated and 12 downregulated ones) at 72 h. During adipogenic differentiation, the number of differentially expressed TFs was 58 at 3 h, 26 at 12 h, and 36 at 72 h. During chondrogenic differentiation, the number was 49 at 3 h, 38 at 12 h, and 39 at 72 h. Remarkably, at 3 and 12h after differentiation, the number of downregulated TFs was higher than the upregulated ones. Next, we calculated the number of differentially expressed TFs at each time point compared with the preceding time point and in three directions of differentiation, the number of TFs that changed per hour decreased as time extends (Figure 2C). Differentially expressed TFs in three time points during three-directional differentiation were summarized in Table 1.

FIGURE 2

Figure 2. Dynamic transcriptional changes upon osteogenic, adipogenic, and chondrogenic differentiation of rat BMSCs. (A) Principal component analyses of the gene expression profile in osteogenic, adipogenic, and chondrogenic groups at three time points; (B) number of significant differential TFs relative to the previous time point (based on three independent experiments); and (C) the number of differentially expressed TFs per hour compared to the previous time point (based on three independent experiments). AD, adipogenesis; OS, osteogenesis; CH, chondrogenesis.

TABLE 1

Table 1. The same differential expressed TFs at three time points after differentiation into osteoblasts, adipocytes, and chondroblasts.

The volcano plot and heatmap of differentially expressed TFs at 3, 12, and 72 h of osteogenesis were shown in Figure 3, while results of comparison between every two lineage-commitment at 3 h after inductive differentiation were presented in Figure 4. Supplementary Figures 1, 2 described the visualized results of differentially expressed TFs at three time points of adipogenic and chondrogenic differentiation.

FIGURE 3

Figure 3. Volcano plot and heatmap of differentially regulated TFs during osteogenesis at 3 h (A,B), 12 h (C,D), and 72 h (E,F) compared with undifferentiated cells (t = 0 h).

FIGURE 4

Figure 4. Volcano plot and heatmap of differentially regulated TFs between every two lineages at 3 h after inductive differentiation (t = 0 h). (A,B) comparison between osteogenesis and adipogenesis; (C,D) comparison between osteogenesis and chondrogenesis, and (E,F) comparison between chondrogenesis and adipogenesis. AD, adipogenesis; OS, osteogenesis; CH, chondrogenesis.

Identification of Specific Differentially Expressed Transcription Factors in Three Phases

In osteogenesis of BMSCs, there were 13, 20, and 14 TFs being identified to be differentially expressed only at 3, 12, and 72 h after inductive differentiation (Figure 5A). In adipogenesis, the numbers were 33, 3, and 14 (Figure 5B), and for chondrogenesis, the numbers were 20, 11, and 13 (Figure 5C). Notably, 11 (e.g., Hopx and Tsc22d3), 12 (e.g., Irf5 and Rab24), and 11 (e.g., Gbx2 and Epas1) TFs showed differentially expressed levels at all three time points in osteogenesis, adipogenesis, and chondrogenesis, respectively.

FIGURE 5

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

Comparison of Differentially Expressed Transcription Factors Between Different Lineages

At 3 h after induction of differentiation, it was noticed that there were 15, 24, and 17 TFs only differentially expressed in osteogenesis, adipogenesis, and chondrogenesis, respectively (Figure 5D). At 12 h, the numbers were 35, 6, and 18 (Figure 5E), and at 72 h, the numbers were 13, 12, and 19 (Figure 5F). Notably, 11 (e.g., Hmg20a and Hoxc6), 5 (e.g., Id2 and Id3), and 11(e.g., Tfcp2l1 and Scx) TFs showed differential expression at all three directions of differentiation at 3, 12, and 72 h after inductive differentiation, respectively.

Functional Enrichment Analysis of Differentially Expressed Transcription Factors

GO and KEGG pathway enrichment analysis of differentially expressed TFs at different time points of three directions of differentiation was performed to identify the most relevant biological processes (BPs), molecular functions (MFs), cellular components (CCs), and pathways. The top enriched terms in BP, CC, MF, and KEGG at 3 h after osteogenesis were presented in Figure 6A. For instance, the enriched BPs terms included “negative regulation of transcription from RNA polymerase II promoter,” “regulation of transcription, DNA-templated,” and “positive regulation of transcription from RNA polymerase II promoter.” The enriched terms at 12 and 72 h after osteogenesis were presented in Figures 6B,C, respectively, while results of the other two directions of differentiation were presented in Supplementary Figures 3, 4.

FIGURE 6

Figure 6. GO including biological process, cellular components and molecular function analysis, and KEGG enrichment analysis of differentially expressed TFs at 3 h (A), 12 h (B), and 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h).

The differentially expressed TFs between osteogenesis and adipogenesis at 3 h after inductive differentiation were also significantly enriched in BPs terms such as “negative regulation of transcription from RNA polymerase II promoter,” “regulation of transcription, DNA-templated.” Results of functional enrichment analysis of the differentially expressed TFs between every two lineage-commitment at 3 h of inductive differentiation were presented in Figure 7.

FIGURE 7

Figure 7. GO including biological process, cellular components and molecular function analysis, and KEGG enrichment analysis of differentially expressed TFs 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.

Protein Interaction Network Construction of Differentially Expressed Transcription Factors

The interactions between the proteins expressed from differentially expressed TFs at 3 h of osteogenesis consisted of 17 nodes and 13 edges, with Klf4 gene being identified to have a relatively high connectivity degree (Figure 8A). The PPI networks at 12 and 72 h after osteogenesis were presented in Figures 8B,C, respectively. Meanwhile, the PPI networks of the other two-direction differentiation at three time points were presented in Supplementary Figure 5.

FIGURE 8

Figure 8. Protein-protein interaction (PPI) network analysis of differentially expressed TFs. PPI network of differentially expressed TFs at 3 h (A), 12 h (B), and 72 h (C) of osteogenesis compared with undifferentiated cells (t = 0 h); the PPI network of differentially regulated TFs between osteogenesis and adipogenesis (D), between osteogenesis and chondrogenesis (E), and between chondrogenesis and adipogenesis (F) at 3 h after inductive differentiation (t = 0 h). Red indicates upregulated TFs, and green indicates the downregulated TFs.

The interactions between the proteins expressed from differentially expressed TFs between osteogenesis and adipogenesis at 3 h after inductive differentiation consisted of 6 nodes and 6 edges, with the Sox9 gene being identified to have a relatively high connectivity degree (Figure 8D). The interactions between the proteins expressed from differentially expressed TFs between osteogenesis and chondrogenesis, and between adipogenesis and chondrogenesis at 3 h after inductive differentiation were presented in Figures 8E,F.

Construction of Differentially Expressed Transcription Factors and Differentially Expressed Genes Regulatory Network

Based on the integrated analysis between TFs and DEGs, it could be noticed that during osteogenesis, the differentially expressed TFs at 3 h were able to regulate 111, 92, and 46 DEGs at 3, 12, and 72 h after differentiation (Figures 9A–C). Four differentially expressed TFs (Tead3, Cenpa, Epas1, and Klf4) identified at 3 h could function at all three time points. Meanwhile, during adipogenesis, the differentially expressed TFs at 3 h were able to regulate 113, 85, and 61 DEGs at 3, 12, and 72 h after differentiation (Figures 10A–C). Six differentially expressed TFs (Sox9, Snai2, Epas1, Id2, Cenpa, and Bhlhe40) occurred in three regulatory networks. For chondrogenesis, the differentially expressed TFs at 3 h were able to regulate 115, 33, and 60 DEGs at 3, 12, and 72 h after differentiation (Figures 11A–C). Six differentially expressed TFs (Runx3, Hic1, Scx, Sox9, Epas1, and Nfia) identified at 3 h could function at all three time points.

FIGURE 9

Figure 9. The regulatory network between differentially expressed TFs at 3 h and targeted differentially expressed genes at 3 h (A), 12 h (B), and 72 h (C) in osteogenesis compared with undifferentiated cells (t = 0 h). Blue indicates TFs, and green indicates targeted genes.

FIGURE 10

Figure 10. The regulatory network between differentially expressed TFs at 3 h and targeted differentially expressed genes at 3 h (A), 12 h (B), and 72 h (C) in adipogenesis compared with undifferentiated cells (t = 0 h). Blue indicates TFs, and green indicates targeted genes.

FIGURE 11

Figure 11. The regulatory network between differentially expressed TFs at 3 h and targeted differentially expressed genes at 3 h (A), 12 h (B), and 72 h (C) in chondrogenesis compared with undifferentiated cells (t = 0 h). Blue indicates TFs, and green indicates targeted genes.

The regulatory networks between TFs and DEGs at 12 and 72 h after differentiation were presented in Supplementary Figure 6. The regulatory networks of the differentially expressed TFs and DEGs between every two lineage-commitment at 3 h after inductive differentiation were presented in Supplementary Figure 7.

Discussion

Most previous studies investigating the differentiation process of BMSCs used primary cells isolated from humans. BMSCs isolated from rat functions across barriers and could be used as a source of regenerative therapy for bone metabolic diseases and bone defects. van de Peppel et al. (2017) discriminated the first 4 days of human BMSCs differentiation into three distinct phases and presented candidate genes for early regulation of osteogenic and adipogenic lineage-commitment. In the present study, we focused on the mechanism underlying the determination of osteogenesis, adipogenesis, and chondrogenesis of rat BMSCs. Primary rat BMSCs represent a heterogeneous cell population, but by using inductive differentiation mediums, the clock of the cells was synchronized to limit cell-cycle heterogeneity. Three time points (3, 12, and 72 h) in three phases of early lineage-commitment of rat BMSCs differentiation were chosen and the high-throughput sequencing data we generated discovered the key transcription factors guiding the differentiation direction. These processes occurred before time windows examined by most previous studies but were proved to play critical roles in the differentiation (Ng et al., 2008; Piek et al., 2010). Hopx and other early responder TFs may control the osteogenic cell fate of BMSCs and participate in the development of osteoporosis. Gbx2 and other early responder TFs should be considered in mechanistic models that clarify cartilage-anabolic changes in the clinical progression of osteoarthritis.

In the first phase of differentiation, gene expression analyses revealed a high number of differentially expressed TFs in all three directions of differentiation. In the regulatory model derived from our data, the first phase represents the initiation stage of the differentiation program and is characterized by expression changes of many transcription-related genes that regulate lineage commitment and set the stage for further differentiation toward a stable phenotype. Downstream analyses of the differentially expressed TFs identified at 3 h of osteogenic differentiation showed that they were capable of regulating 92 and 46 DEGs in the second and third phases, respectively. As for the second phase, differentiating lineages begin to deviate, as reflected by the fact that many transcriptional changes are direct targets of the differentially expressed TFs in the first phases. In the current study, TFs that changed per hour decreased as time extended, which suggested that the differentiating cells had reached a stable phenotype. This illustrates the transition from proliferation to differentiation and reflects the inverse correlation between these processes as has been described in various other differentiating cells.

Interestingly, we identified two TFs (Tsc22d3 and Epas1) that were upregulated in three lineages at all three time points, which may play an important role in the initiation of osteogenesis, adipogenesis, and chondrogenesis. Tsc22d3, also known as glucocorticoid−induced leucine zipper, is a glucocorticoid−responsive anti−inflammatory molecule that could regulate intracellular signaling pathways via hypoxia-induced factor-1α (HIF−1α) as well as AP−1 (Kanda et al., 2020). Tsc22d3 is also one of the most regulated genes in Cushing’s syndrome (CS) and has been demonstrated to regulate osteoblast and bone turnover. Epas1, also known as HIF-2, is an active nucleoprotein and transcription factor induced by hypoxia. It could regulate hundreds of genes including VEGF, Ang-2, and HO-1, acting to enhance the production of corresponding proteins, and launch the process of angiogenesis and osteogenesis (Zhou et al., 2020). Thus, transfection of Epas1 into BMSCs provides an alternative for treating bone defects in the clinic. Other noteworthy TF is Hopx, which is the smallest member of the homeodomain-containing protein family (Hng et al., 2020). Hopx manifested an up-regulated transcription level only under the condition of osteogenesis (3, 12, and 72 h). Unlike other typical homeobox proteins, Hopx binds to different protein partners and therefore recruits transcription factors to gene promoter regions (Trivedi et al., 2010; Hng et al., 2020). Hemming et al. (2016) have demonstrated a correlation between the Enhancer of zeste homolog 2 (EZH2) and Hopx, implicating a regulatory effect of Hopx during BMSC osteogenesis. Similarly, Gbx2 may participate in the decision of chondrogenesis of BMSCs. This TF is a LIF/STAT3 signaling downstream target that could allow long-term expansion of the undifferentiated embryonic stem cells (Wang et al., 2017). Sox9 was up-regulated in chondrogenesis but down-regulated in adipogenesis in tree time points, which was consistent with previous studies that this TF was essential for lineage differentiation of BMSCs (Lefebvre et al., 2019). Cooperation between GLI, JUN/FOSL2, and Sox9 is based on multiple binding sites near those of Sox9 and therefore transactivates many cartilage-specific genes (Liu and Lefebvre, 2015). Four TFs including Ybx3, Zbtb25, Hes6, and Zfp879 were inclusively upregulated at 3 h of osteogenesis. Further studies are needed to elucidate their role in the initiation of osteogenesis.

Then we used a variety of bioinformatics analysis methods to analyze these differentially expressed TFs to obtain a comprehensive understanding. As a result, GO terms of the differentially expressed TFs at 3 h of osteogenesis, adipogenesis, and chondrogenesis for biological process categories included “negative regulation of transcription from RNA polymerase II promoter,” “transcription, DNA-templated,” and “regulation of transcription, DNA-templated.” Furthermore, the KEGG pathway enrichment analysis revealed that the differentially expressed TFs at 3 h of osteogenesis and adipogenesis were predominantly associated with the “signaling pathways regulating pluripotency of stem cells” (Id2, Id3, and Klf4). Id2 and Id3 were downregulated in three lineages while Klf4 was upregulated in the process of adipogenesis and osteogenesis. This suggests that the initiation of BMSCs differentiation is similarly activated in both lineages, and changes of these signaling pathways and TFs are necessary to exit the immature multi-potent cell stage and acquire a specialized mesenchymal phenotype. While in the second phase, the enriched KEGG pathways begin to deviate. For osteogenesis, “signaling pathways regulating pluripotency of stem cells” was still significantly enriched, while for adipogenesis, “PI3K-Akt signaling pathway” and “MAPK signaling pathway” were significantly enriched. For chondrogenesis, “cAMP signaling pathway” was significantly enriched.

The regulatory network provided potential connections between identified DEGs and TFs. There were four, six, and six differentially expressed TFs identified at 3 h during osteogenesis, adipogenesis, and chondrogenesis possessing targets among DEGs identified at all three time points. It could be noticed that Cenpa and Epas1 overlap in osteogenesis and adipogenesis while Sox9 overlaps in adipogenesis and chondrogenesis. Like Epas1, Cenpa was also upregulated in both osteogenesis and adipogenesis at 3 h. Cenpa is a functionally conserved molecule about 17kDa that forms a centromere-specific nucleosome with H2A, H2B, and H4, whose overexpression is crucial for multiple kinds of cancers (Saha et al., 2020; Zhang et al., 2020). The role of Cenpa in lineage-commitment of BMSCs has not to be raised by former studies.

It could be noticed that there are some variations at the transcriptional level between the differentiation of human and rat BMSCs. At present, mechanism investigation regarding the lineage-commitment of rat BMSCs is still lacking. We performed high-throughput transcriptome sequencing and conducted a rigorous informatics analysis based on multiple databases, proposing multiple key transcription factors that may participate in the lineage-commitment of rat BMSCs. These findings provide research directions and potential molecular candidates for further studies investigating the underlying mechanisms of three-directional differentiation of rat BMSCs. However, this study has some shortcomings. First, this is an in vitro experiment based on isolated rat BMSCs. The differentiation processes may possess different characteristics in vivo. Findings obtained from in vitro experiments cannot simply be transposed to predict the reaction of an entire organism in vivo. Second, only three time points were selected in this study. These time points were separately locating in three phases and were used as delegates. Third, further in-depth studies are needed to establish whether MSCs from adipose tissue or other anatomical locations undergo a similar cascade of transcriptional events. Adipose-derived MSCs may respond very differently with respect to their clonogenic potential and doubling time in comparison with BMSCs (Russo et al., 2014). Fourth, the detailed functional mechanism of identified TFs was not clearly elucidated, and we could not attribute the lineage-commitment to a particular TF (or set of TFs). Last but not least, because there was a large and diverse amount of data, we could not validate these TFs by western blotting. Proteomics study may be needed in the future.

Conclusion

These transcription factors, functional terms and predicted targeted genes identified in the current study can help shed light on the molecular mechanism underlying the lineage-commitment of rat bone mesenchymal stem cells. Changes in TF activity during the early stage may play pivotal roles in the regulation of subsequent later phases of mesenchymal differentiation. Hopx and other early responder TFs may control the osteogenic cell fate of BMSCs, while Gbx2 and other early responder TFs should be considered in mechanistic models that clarify cartilage-anabolic changes. The latter finding exhibits possibilities for efficiently inducing osteoblast and chondrocytes differentiation, as part of a bone anabolic strategy for osteoporosis and osteoarthritis. In-depth experiments are required to further clarify the effect of these transcription factors.

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/query/acc.cgi?acc=GSE185140.

Ethics Statement

The animal study was reviewed and 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 No. 201808080126), the incubation fund of Shandong Provincial Hospital (Grant No. 2020FY019), the Young Scholars Program of Shandong Provincial Hospital, the Young Taishan Scholars Program of Shandong Province (Grant No. tsqn201909183), the Academic Promotion Program of Shandong First Medical University (Grant No. 2020RC008), and the Natural Science Foundation of Shandong Province (Grant Nos. ZR201911090016 and ZR202102180575).

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/fcell.2021.768316/full#supplementary-material

Supplementary Figure 1 | Volcano plot and heatmap of differentially regulated TFs during adipogenesis at 3 h (A,B), 12 h (C,D), and 72 h (E,F) compared with undifferentiated cells (t = 0 h).

Supplementary Figure 2 | Volcano plot and heatmap of differentially regulated TFs during chondrogenesis at 3 h (A,B), 12 h (C,D), and 72 h (E,F) compared with undifferentiated cells (t = 0 h).

Supplementary Figure 3 | GO including biological process, cellular components and molecular function analysis, and KEGG enrichment analysis of differentially expressed TFs at 3 h (A), 12 h (B), and 72 h (C) of adipogenesis compared with undifferentiated cells (t = 0 h).

Supplementary Figure 4 | GO including biological process, cellular components and molecular function analysis, and KEGG enrichment analysis of differentially expressed TFs at 3 h (A), 12 h (B), and 72 h (C) of chondrogenesis compared with undifferentiated cells (t = 0 h).

Supplementary Figure 5 | Protein-protein interaction (PPI) network analysis of differentially expressed TFs at 3 h (A), 12 h (B), and 72 h (C) of osteogenesis, and 3 h (D), 12 h (E), and 72 h (F) of chondrogenesis compared with undifferentiated cells (t = 0 h). Red indicates upregulated TFs, and green indicates the downregulated TFs.

Supplementary Figure 6 | The regulatory network between differentially expressed TFs and targeted differentially expressed genes at 12 and 72 h in osteogenesis (A,B), adipogenesis (C,D), and chondrogenesis (E,F) compared with undifferentiated cells (t = 0 h). Blue indicates TFs, and green indicates targeted genes.

Supplementary Figure 7 | The regulatory network between differentially expressed TFs and targeted 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. Blue indicates TFs, and green indicates targeted genes.

Footnotes

1. ^ http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/
2. ^ http://www.geneontology.org/
3. ^ https://david.ncifcrf.gov/summary.jsp
4. ^ http://string-db.org

References

Arthur, A., and Gronthos, S. (2020). Clinical application of bone marrow mesenchymal Stem/Stromal cells to repair skeletal tissue. Int. J. Mol. Sci. 21:21249759. doi: 10.3390/ijms21249759

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruedigam, C., Driel, M., Koedam, M., Peppel, J., van der Eerden, B. C., 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 Cell Biol. 1:17. doi: 10.1002/9780470151808.sc01h03s17

PubMed Abstract | CrossRef Full Text | Google Scholar

Carballo, C. B., Nakagawa, Y., Sekiya, I., and Rodeo, S. A. (2017). Basic science of articular cartilage. Clin. Sports Med. 36, 413–425. doi: 10.1016/j.csm.2017.02.001

PubMed Abstract | 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, 1128–1139. doi: 10.1038/cdd.2015.168

PubMed Abstract | CrossRef Full Text | Google Scholar

Crapnell, K., Blaesius, R., Hastings, A., Lennon, D. P., Caplan, A. I., and Bruder, S. P. (2013). Growth, differentiation capacity, and function of mesenchymal stem cells expanded in serum-free medium developed via combinatorial screening. Exp. Cell Res. 319, 1409–1418. doi: 10.1016/j.yexcr.2013.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, F., Li, Z., Zhao, Q., Xiong, C., and Ni, K. (2020). Analysis of multi-lineage gene expression dynamics during primordial germ cell induction from human induced pluripotent stem cells. Stem Cell Res. Ther. 11:100. doi: 10.1186/s13287-020-01620-y

PubMed Abstract | 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, 732–739. doi: 10.1007/s13770-016-0013-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hemming, S., Cakouros, D., Vandyke, K., Davis, M. J., Zannettino, A. C., and Gronthos, S. (2016). Identification of novel EZH2 targets regulating osteogenic differentiation in mesenchymal stem cells. Stem Cells Dev. 25, 909–921. doi: 10.1089/scd.2015.0384

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: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. Int. J. Mol. Sci. 19:360. doi: 10.3390/ijms19020360

PubMed Abstract | CrossRef Full Text | Google Scholar

Huynh, N., Zhang, B., and Guilak, F. (2019). High-depth transcriptomic profiling reveals the temporal gene signature of human mesenchymal stem cells during chondrogenesis. FASEB J. 33, 358–372. doi: 10.1096/fj.201800534R

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:e12956. doi: 10.1111/cpr.12956

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanda, A., Hirose, I., Noda, K., Murata, M., and Ishida, S. (2020). Glucocorticoid-transactivated TSC22D3 attenuates hypoxia- and diabetes-induced Muller glial galectin-1 expression via HIF-1alpha destabilization. J. Cell. Mol. Med. 24, 4589–4599. doi: 10.1111/jcmm.15116

PubMed Abstract | CrossRef Full Text | Google Scholar

Lefebvre, V., Angelozzi, M., and Haseeb, A. (2019). SOX9 in cartilage development and disease. Curr. Opin. Cell Biol. 61, 39–47. doi: 10.1016/j.ceb.2019.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Ezzelarab, M. B., and Cooper, D. K. (2012). Do mesenchymal stem cells function across species barriers? Relevance for xenotransplantation. Xenotransplantation. 19, 273–285. doi: 10.1111/xen.12000

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. F., and Lefebvre, V. (2015). The transcription factors SOX9 and SOX5/SOX6 cooperate genome-wide through super-enhancers to drive chondrogenesis. Nucleic Acids Res. 43, 8183–8203. doi: 10.1093/nar/gkv688

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. J. Inflamm. Res. 14, 3537–3554. doi: 10.2147/JIR.S320839

PubMed Abstract | CrossRef Full Text | Google Scholar

Merrett, J. E., Bo, T., Psaltis, P. J., and Proud, C. G. (2020). Identification of DNA response elements regulating expression of CCAAT/enhancer-binding protein (C/EBP) beta and delta and MAP kinase-interacting kinases during early adipogenesis. Adipocyte 9, 427–442. doi: 10.1080/21623945.2020.1796361

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, F., Boucher, S., Koh, S., Sastry, K. S., Chase, L., Lakshmipathy, U., et al. (2008). PDGF, TGF-beta, and FGF signaling is important for differentiation and growth of mesenchymal stem cells (MSCs): Transcriptional profiling can identify markers and signaling pathways important in differentiation of MSCs into adipogenic, chondrogenic, and osteogenic lineages. Blood 112, 295–307. doi: 10.1182/blood-2007-07-103697

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Connor, K. C. (2019). Molecular profiles of Cell-to-Cell variation in the regenerative potential of mesenchymal stromal cells. Stem. Cells Int. 2019:5924878. doi: 10.1155/2019/5924878

PubMed Abstract | CrossRef Full Text | Google Scholar

Piek, E., Sleumer, L. S., van Someren, E. P., Heuver, L., de Haan, J. R., de Grijs, I., et al. (2010). Osteo-transcriptomics of human mesenchymal stem cells: Accelerated gene expression and osteoblast differentiation induced by vitamin D reveals c-MYC as an enhancer of BMP2-induced osteogenesis. Bone 46, 613–627. doi: 10.1016/j.bone.2009.10.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Pittenger, M. F., Discher, D. E., Peault, B. M., Phinney, D. G., Hare, J. M., and Caplan, A. I. (2019). Mesenchymal stem cell perspective: Cell biology to clinical progress. NPJ Regen. Med. 4:22. doi: 10.1038/s41536-019-0083-6

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. Cell Dev. Biol. 8:561. doi: 10.3389/fcell.2020.00561

PubMed Abstract | CrossRef Full Text | Google Scholar

Russo, V., Yu, C., Belliveau, P., Hamilton, A., and Flynn, L. E. (2014). Comparison of human adipose-derived stem cells isolated from subcutaneous, omental, and intrathoracic adipose tissue depots for regenerative applications. Stem Cells Transl. Med. 3, 206–217. doi: 10.5966/sctm.2013-0125

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. Curr. Pharm. Biotechnol. 22, 200–215. doi: 10.2174/1389201021666200907121530

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, A. K., Contreras-Galindo, R., Niknafs, Y. S., Iyer, M., Qin, T., Padmanabhan, K., et al. (2020). The role of the histone H3 variant CENPA in prostate cancer. J. Biol. Chem. 295, 8537–8549. doi: 10.1074/jbc.RA119.010080

PubMed Abstract | CrossRef Full Text | Google Scholar

Trivedi, C. M., Zhu, W., Wang, Q., Jia, C., Kee, H. J., Li, L., et al. (2010). Hopx and Hdac2 interact to modulate Gata4 acetylation and embryonic cardiac myocyte proliferation. Dev. Cell 19, 450–459. doi: 10.1016/j.devcel.2010.08.012

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 Cell Rep. 8, 947–960. doi: 10.1016/j.stemcr.2017.02.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Tang, L., Liu, D., Ying, Q. L., and Ye, S. (2017). The transcription factor Gbx2 induces expression of Kruppel-like factor 4 to maintain and induce naive pluripotency of embryonic stem cells. J. Biol. Chem. 292, 17121–17128. doi: 10.1074/jbc.M117.803254

PubMed Abstract | CrossRef Full Text | Google Scholar

Yianni, V., and Sharpe, P. T. (2020). Epigenetic mechanisms driving lineage commitment in mesenchymal stem cells. Bone 134:115309. doi: 10.1016/j.bone.2020.115309

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, B., and Wang, C. Y. (2016). Osteoporosis: The result of an ‘Aged’ bone microenvironment. Trends Mol. Med. 22, 641–644. doi: 10.1016/j.molmed.2016.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Guo, F., Wei, J., Xian, M., Tang, S., Zhao, Y., et al. (2017). An integrated approach to identify critical transcription factors in the protection against hydrogen peroxide-induced oxidative stress by Danhong injection. Free Radic. Biol. Med. 112, 480–493. doi: 10.1016/j.freeradbiomed.2017.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q. Y., Li, Z. R., Gao, F. Q., and Sun, W. (2018). Pericollapse stage of osteonecrosis of the femoral head: A last chance for joint preservation. Chin. Med. J. 131, 2589–2598. doi: 10.4103/0366-6999.244111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Yang, L., Shi, J., Lu, Y., Chen, X., and Yang, Z. (2020). The oncogenic role of CENPA in hepatocellular carcinoma development: Evidence from bioinformatic analysis. Biomed. Res. Int. 2020:3040839. doi: 10.1155/2020/3040839

PubMed Abstract | 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 alpha inhibitor (HIF1AN). Mech. Dev. 162:103613. doi: 10.1016/j.mod.2020.103613

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R., Guo, F., Xiang, C., Zhang, Y., Yang, H., and Zhang, J. (2021). Systematic Study of Crucial Transcription Factors of Coptidis rhizoma Alkaloids against Cerebral Ischemia-Reperfusion Injury. ACS Chem. Neurosci. 12, 2308–2319. doi: 10.1021/acschemneuro.0c00730

PubMed Abstract | CrossRef Full Text | Google Scholar