This is the first high throughput microarray study on BMP6 induced transcriptional program in human MSC. It covers the whole process from early to late stage osteoblast differentiation and mineralization. We conducted a comprehensive gene set analysis to identify relevant regulatory mechanisms and functional groups. We inferred a series of significant KEGG pathways, GO terms and experimental sets at different stages of BMP6 induction process. We not only showed which pathways or gene sets are significant, but also when and how they are involved in the osteoblast differentiation and mineralization. Different from common pathway analyses [13, 14, 16], our work further captures the interconnections among individual pathways or functional groups and integrate them into a whole system. Taken together, this work provides clearer mechanistic picture of osteoblast differentiation and function.
We inferred novel and coherent sets of regulatory mechanisms downstream of BMP6 signaling during osteoblast differentiation and mineralization. First, the same set of KEGG pathways are constantly involved in BMP6 induction. Their roles in osteogenic induction are clarified based on their perturbation patterns and connected to relevant discoveries in literature. These significant KEGG pathways are not separated but rather they work as a unified super regulatory system, and the pathway perturbation patterns we derived reflect a dynamic transmission process of the regulatory signal at transcriptional level along the super system. Second, a varying set of GO processes and functional groups are involved at different stage of BMP6 induced osteoblast differentiation and mineralization. These suggest novel yet plausible regulatory mechanisms, which are connected to but have not directly and explicitly introduced in literature works. Third, the most significant experimental sets suggest novel transcriptional regulators including MYB and BAF57, and regulatory pathways consistent with predictions based on KEGG and GO gene sets above.
Connections between KEGG pathways are evident as shown in the super regulatory network of pathways (Figure 4 and Table 2). Perturbations propagate along the super network at two levels: at protein level, the phosphorylation, binding, activation/inhibition events relay along pathways and transmit into interconnected pathways, as stated by KEGG graphs (Figure 3); at transcriptional level, gene expression perturbation propagates through auto-regulatory (feedback and feed forward) loops within pathways, and bridges into its neighbor pathways through the multiple component genes in common, as suggested by our pathway analysis results (Figure 4 and Table 2). Protein level transmission is faster, but transcriptional level transmission lasts longer hence ensure the long term biological effects. These longer lasting transcriptional effects are clearly evident in our data as the gene expression levels were perturbed long after the withdrawal of BMP6 treatment. The hard wired KEGG pathways and interconnections between them define how BMP6 signal triggers downstream programs (Figure 4).
There are also connections between BMP6 signal and the processes/groups represented by GO terms and experimental sets. For example, Notch signal and IGF signal are involved in the whole induction process (Figure 5a) like all significant KEGG pathways (Figure 2). It follows that these two signals should be also part of the super regulatory system and interconnected with multiple significant KEGG pathways of the network (Figure 4). This hypothesis is well supported by our data and literature: (1) Notch signal directly interacts with BMP signal. SMAD1 and NIC synergize to induce expression of HEY1 and other Notch targets [44–46]. Indeed up-regulation of HEY1 (and HEY2) requires continued BMP6 treatment by 96 hours (Additional file 1: Supplementary Figure 4c). Besides this binding-synergy at protein level, Notch ligand JAG1 expression is also up-regulated directly by BMP6 (Additional file 1: Supplementary Figure 4c). Notch also interacts with Wnt signal . (2) Growth hormone (GH) signal  as part of cytokine-cytokine receptor interaction (KEGG pathway 04060) and Jak-STAT pathways (KEGG pathway 04630) are activated and by BMP signal (Figure 4) with GHR up-regulated (Additional file 1: Supplementary Figure 4a-b), GH signal up-regulates (IGF1 and IRS1, Additional file 1: Supplementary Figure 4d) and activates IGF signal in turn . Connections between BMP signal and the two predicted transcriptional regulators, MYB and BAF57, are described above in the Results.
When examined together, we find a consistent picture emerging from the lists of significant KEGG pathways, GO terms, and experimental sets. For example, KEGG focal adhesion and ECM receptor interaction pathways (Table 1), GO homophilic cell adhesion (Table 3) and extracellular matrix structural constituent (not shown) groups consistently show the relevance of cell adhesion and extracellular matrix in osteoblast differentiation and mineralization. GO immune response groups (Table 3) echoes KEGG cytokine-cytokine receptor interaction pathway (Table 1). Similarly, significant experimental target genes sets (Table 4) closely reflected changes in the regulatory KEGG pathways (Table 1) or GO processes/groups (Table 3).
Interestingly, there are discrepancies among the clustering and the significant KEGG pathways, GO terms, and experimental sets. For example, Notch signaling is defined as both a KEGG pathway and a GO process. This KEGG pathway is not significant (not shown) but this GO process is (Figure 5). This discrepancy arises from two sources: (1) different definitions, i.e. KEGG pathways contain partially different set of genes from corresponding GO processes. While KEGG pathways tend to cover the whole homeostatic signal transmission systems even across multiple transcriptional cycles, GO usually covers one or multiple discrete steps or functional groups for a process. KEGG and GO definitions can be considered complementary and both provide valuable gene sets for our analysis. (2) GAGE  treats KEGG pathways and GO term gene sets differently: genes under a GO term are taken as a group coregulated towards a single direction, either all up or all down regulated, whereas genes in a KEGG pathway are frequently not coregulated and expression changes in both directions are counted. Timing discrepancies exist between experimental sets and corresponding KEGG pathways. For example, IFN positive target sets (only IFN beta shown) are not significant at 24-96 hours and MYB target set not at 96 hours (Figure 6) while Jak-STAT pathway and Wnt signaling pathway are significantly perturbed all the time (Figure 2). This can be explained by the fact that two-directional perturbation treatment for KEGG pathway does not account for direction or net effect of the perturbation, whether inhibited, activated or no overall effect. In the other hand, GO term analysis has no such issue, and IRS negative set and corresponding GO IGF receptor binding group are both significant all the time.
In this work, we took a systems approach in studying MSC differentiation. We combined experimental and computational work to reconstruct a unified picture of BMP6 signaling. The same set of experimental design and computational approach could be used to study other physiological and pathological processes, such as the differentiation of other cell types or tumorogenesis.