DEGs and KEGG pathway of OS induction
To systematically explore gene expression changes during the early stage of hMSCs differentiated into OS, we identified the DEGs by comparing the OS_2d, OS_3d, OS_4d groups with the group OS_0d. There were 452 DEGs in the OS_2d group, including 262 up-regulations and 190 down-regulations while 656 DEGs in the OS_3d group, of which 350 were up-regulated and 306 were down-regulated (Figure 1A and B). In addition, a total of 687 genes with 378 up-regulations and 309 down-regulations were detected to be differentially expressed in the OS_4d group (Figure 1C). Moreover, the common upregulations were nearly twice compared to downregulations in these three time points, which might mean that the OS was the result of the upregulation of pro-differentiation genes and activation of related pathways (Figure S1A and B).To explore the potential mechanisms in the OS initiative differentiation, KEGG pathway enrichment analysis were performed with DAVID database. The top 10 pathways sorted by the number of DEGs enriched in every term were displayed via scatter plot (Figure 1D-F). Obviously, PI3K-Akt signaling pathway was the most significant pathway in each time point during the OS differentiation.
DEGs and KEGG pathway of AD induction
There were a total of 729 DEGs in AD_2d group, 870 DEGs in AD_3d group and 905 DEGs in AD_4d group, including 364, 441 and 461 up-regulations while 365, 429 and 444 down-regulations in the AD_2d, AD_3d and AD_4d group, respectively (Figure 2A-C). The count of DEGs in AD induction showed little change in each time point as well, but increased significantly compared to DEGs OS, which might be due to hMSCs favor differentiate into AD rather than OS [35, 36]. As shown in the Venn diagrams, the number of the down-regulations was similar to the up-regulations in the AD, but was more than twice that in the OS, meaning the down-regulation play more significance in the AD. (Figure S1A and B). As the same as OS induction, the top 10 KEGG pathways enrichment at each time point were visualized in the scatter plot (Figure 2D-F). According to numbers of enriched DEGs, the PI3K-Akt signaling pathway was also ranked at first, ignoring the pathways in cancer usually enriched most genes in the DAVID. Therefore, we surmised that the PI3K-Akt signaling pathway also exerted a crucial effect on AD from BMSCs.
Key KEGG pathway enrichment
We next investigated the key signaling pathway in these two lineages of differentiation. The PI3K-Akt signaling pathway, FoxO signaling pathway, Focal adhesion and Cell cycle were the common signaling pathways in OS differentiation (Figure 3A). In AD induction, PI3K-Akt signaling pathway, Rap1 signaling pathway, Focal adhesion, Regulation of actin cytoskeleton, Cell cycle and HTLV-I infection were all enriched in the three groups (Figure 3C). Three common pathways including PI3K-Akt signaling pathway, Focal adhesion and Cell cycle were all existed in the OS and AD induction. Therefore, we hypothesized that the PI3K-AKT signaling pathway, as well as Focal adhesion, and Cell cycle were the crucial mechanism for initiating the differentiation of hMSCs. For the differences between the two lineages differentiation, the FoxO signaling pathway was specific for OS induction and the Rap1 signaling pathway belonged only to AD induction, indicating that FoxO signaling pathway and Rap1 signaling pathway were the specific signaling pathway for the OS and AD induction, respectively. And the expression level of DEGs enriched in the FoxO signaling pathway and Rap1 signaling pathway were showed at heatmap (Figure 3B and D).
Key genes in the OS induction
Interactions between multiple genes could be well understood by PPI analysis. The DEGs OS were used to perform the PPI analysis with Between and Stress algorithm of cytoHubba, from which we choose the top 20 hub genes to construct the sub-network (Figure 4A-C and Figure S2). As shown in these six PPI networks, the common genes with high degree and clustering coefficients were FoxO3, IL-6, JUN and CAT, which might play important roles in the OS induction. At the same time, candidate genes enriched in the FoxO signaling pathway in each group were also used to perform the PPI network (Figure 4D-F). Similarly, FoxO1, IL-6 and CAT were localized at the center of network. These results indicated that FoxO3, IL-6, JUN and CAT maybe the key genes for early OS differentiation. Since FoxO1, IL-6 and CAT were all enriched in the FoxO signaling pathway, further indicating that FoxO signaling pathway may be the specific signaling pathway for early OS differentiation.
Key genes in the AD induction
DEGs AD were used to perform the PPI networks, and sub-networks were formed by the top 20 hub genes. VEGFA, FGF2, MYC and PTEN had the highest degree and clustering coefficients in each group computed by Between algorithm (Figure 5A-C), while VEGFA, FGF2, MYC, CCND1 and PTEN were the core computed by Stress algorithm (Figure S3), which meant that these genes played core roles in the AD induction. Similarly, the DEGs related to Rap1 signaling pathway were selected to construct PPI network. VEGFA, FGF2 and PIK3R1 were identified as the core genes in the network (Figure 5D-F). These interactions relied strongly on VEGFA and FGF2, which may determine the tendency of BMSCs differentiate into AD. Therefore, hub genes VEGFA, FGF2 and Rap1 signaling pathway might be the crucial mechanisms for early AD differentiation.
GO analysis of OS induction
Gene ontology enrichment analyses were performed by the DEGs from each group, and the significant top 15 of the GO terms including BP, CC and MF were displayed as the bar diagrams (Figure 6A-C). In the MF category, receptor binding, protein binding, growth factor activity were the common GO terms in the groups of OS_2d, OS_3d and OS_4d, while cytoplasm, cytosol, extracellular exosome were the common GO terms in the CC category. In the BP category, signal transduction, positive regulation of transcription and positive regulation of cell proliferation were all enriched in the three time points of OS induction. To further investigate the relationship of the key genes and GO functional categories, the key genes from the FoxO signaling pathway and PI3K-Akt signaling pathway were put to perform the chord plots (Figure 6D-F). Key genes including FoxO3, IL-6, CAT and PIK3R1 mainly clustered into protein binding, membrane, cytosol, nucleus, extracellular space, plasma membrane, cytoplasm, negative and positive regulation of apoptotic process. Suggesting that GO terms in OS induction were a series of biological responses that initiated by ligand-receptor binding and transcriptional information into the nucleus.
GO analysis of AD induction
The top 15 of the GO terms performed by the DEGs from each group of AD were shown in the bar diagrams (Figure 7A-C). The results of the AD differentiation were similar to that in the OS, implied that there were similar biological effects on the two lineages differentiations from hMSCs. Chord plots of AD induction (Figure 7D-F) performed by the DEGs of Rap1 signaling pathway and top 15 of GO terms in each group, showed that FGF2 was mainly enriched in the protein binding, growth factor activity, cytoplasm, nucleus, extracellular space, extracellular region, signal transduction, positive regulation of cell proliferation and positive regulation of cell proliferation. Moreover, another target gene VEGFA was mainly related to protein binding, growth factor activity, membrane, extracellular space, extracellular region, cytoplasm and positive regulation of cell proliferation. This was consistent with previously reported that VEGFA was an extracellular signal molecule, which regulated Rap1 signaling pathway by combining some intracellular signal factors [37, 38].
Osteo-adipogenesis relative gene expression
Expression level of OS (ALP) and AD (LPL) specific genes at day 2-4 was increased comparing to the day 0, suggesting that hBMSCs were successfully induced into OS and AD (Figure 8A and C). The expression levels of the key genes for OS (IL6 and CAT) or AD (PIK3R1, FGF2 and VEGFA) induction were changed dramatically at the first 2 days, but no obvious change was founded after that, indicating that the first 2 days maybe the critical period for the lineage commitment determinant. However, FoxO3 at day 2-4 was significantly increased in the OS induction, revealing that it may be key genes for osteogenic differentiation at a steady period of lineage-continuation (Figure 8A). In addition, the expression of PIK3R1 was higher in both OS and AD, suggesting PIK-AKT signaling pathway is a promoter in the differentiation of BMSCs.
Signaling pathway of OS and AD induction
The expression of the key proteins in the signaling pathway relative to the OS and AD induction were detected using inflorescence staining (Figure 8 B and D). In consistent with the results of PCR, the expression of FoxO3 and CAT in OS induction at day 2-4 was higher compared to the day 0, while IL6 was conversely lower, suggesting that FoxO3 signaling pathway was involved in the osteogenic differentiation (Figure 8B). Additionally, the downregulation of VEGFA and FGF2 in AD induction indicated that Rap1 signaling pathway maybe crucial for AD differentiation. Meanwhile, the expression of PIK3R1 protein was also increased both in OS and AD.