Transcriptomic landscape regulated by the 14 types of bone morphogenetic proteins (BMPs) in lineage commitment and differentiation of mesenchymal stem cells (MSCs).

Mesenchymal stem cells (MSCs) are ubiquitously-existing multipotent progenitors that can self-renew and differentiate into multiple lineages including osteocytes, chondrocytes, adipocytes, tenocytes and myocytes. MSCs represent one of the most commonly-used adult progenitors and serve as excellent progenitor cell models for investigating lineage-specific differentiation regulated by various cellular signaling pathways, such as bone morphogenetic proteins (BMPs). As members of TGFβ superfamily, BMPs play diverse and important roles in development and adult tissues. At least 14 BMPs have been identified in mammals. Different BMPs exert distinct but overlapping biological functions. Through a comprehensive analysis of 14 BMPs in MSCs, we demonstrated that BMP9 is one of the most potent BMPs in inducing osteogenic differentiation of MSCs. Nonetheless, a global mechanistic view of BMP signaling in regulating the proliferation and differentiation of MSCs remains to be fully elucidated. Here, we conducted a comprehensive transcriptomic profiling in the MSCs stimulated by 14 types of BMPs. Hierarchical clustering analysis classifies 14 BMPs into three subclusters: an osteo/chondrogenic/adipogenic cluster, a tenogenic cluster, and BMP3 cluster. We also demonstrate that six BMPs (e.g., BMP2, BMP3, BMP4, BMP7, BMP8, and BMP9) can induce I-Smads effectively, while BMP2, BMP3, BMP4, BMP7, and BMP11 up-regulate Smad-independent MAP kinase pathway. Furthermore, we show that many BMPs can upregulate the expression of the signal mediators of Wnt, Notch and PI3K/AKT/mTOR pathways. While the reported transcriptomic changes need to be further validated, our expression profiling represents the first-of-its-kind to interrogate a comprehensive transcriptomic landscape regulated by the 14 types of BMPs in MSCs.

Introduction

Originally described by Friedentstein over 60 years ago as adherent cells with a fibroblast-like appearance, mesenchymal stem cells (MSCs) are a heterogeneous population of non-hematopoietic multipotent progenitors, which can undergo self-renewal and differentiate into multiple cell types, such as osteocytes, chondrocytes, adipocytes, tenocytes and myocytes upon appropriate stimulations.2, 3, 4, 5, 6, 7 Although originally isolated from bone marrow stromal cells, MSCs have been shown to exist in almost every type of tissue, including periosteum, brain, liver, bone marrow, adipose, skeletal muscle, amniotic fluid and hair follicles.4, 5, 8 In fact, MSCs represent one of the most commonly-used adult progenitor cells due to their wide availability, easiness of isolation and maintenance, and potential applications in regenerative medicine.5, 6 Furthermore, MSCs serve as excellent progenitor cell models for investigating lineage-specific differentiation induced by various cellular signaling pathways, such as bone morphogenetic proteins (BMPs).6, 9, 10, 11, 12, 13, 14

BMPs belong to the TGFβ superfamily, a group of homologous signaling proteins plays diverse and important roles in embryogenesis, organogenesis, cell proliferation and stem cell differentiation.15, 16, 17 Disruptions in BMP signaling can cause a variety of skeletal and extraskeletal anomalies.15, 17, 18 There are at least 14 types of BMPs in humans and rodents.6, 15, 19, 20 Different BMPs exert distinct but overlapping biological functions. For example, BMP7 is involved in proper kidney, eye, and limb development, while BMP4, BMP7 and BMP15 are important for proper reproductive tissue development.15, 17, 18 BMPs play important roles in regulating the proliferation and lineage-specific differentiation of mesenchymal stem cells.6, 17 BMP2 and BMP7 are shown to play a role in chondrogenic differentiation, whereas BMP12, BMP13, and BMP14 may be required for normal tendon healing.6, 17, 22 Furthermore, BMPs are considered as a group of the most potent osteoinductive factors.6, 17, 20 Through a comprehensive analysis of the 14 types of human BMPs in MSCs, we demonstrated that BMP9 (also known as growth differentiation factor 2, or GDF2) is one of the most potent BMPs among the 14 types of BMPs in inducing osteogenic differentiation of MSCs,14, 23, 24, 25 possibly by regulating several important downstream targets during BMP9-induced osteoblast differentiation.14, 26, 27, 28, 29, 30, 31, 32, 33, 34 Nonetheless, the global mechanistic role of BMP signaling in regulating the proliferation and differentiation of MSCs remains to be fully elucidated.

Here, we conducted a comprehensive transcriptomic profiling of the MSCs stimulated by the 14 types of BMPs at the early stage, which allows us to identify the BMP-initiated immediate early transcriptomic changes in MSCs. Hierarchical clustering analysis classifies 14 BMPs into three subclusters: an osteo/chondrogenic/adipogenic cluster, a tenogenic cluster, and BMP3 cluster. We also demonstrate that six BMPs (e.g., BMP2, BMP3, BMP4, BMP7, BMP8, and BMP9) can induce I-Smads effectively, while BMP2, BMP3, BMP4, BMP7, and BMP11 also up-regulate Smad-independent MAP kinase pathway. Furthermore, our transcriptomic analysis reveals that many BMPs can effectively regulate the expression levels of the signal mediators of Wnt signaling, Notch pathway, and PI3K/AKT/mTOR pathway. While further functional validations of the transcriptomic changes, our expression profiling analysis represents the first of its kind to interrogate a comprehensive early stage transcriptomic landscape regulated by the 14 types of BMPs in MSCs.

Materials and methods

Cell culture and chemicals

HEK-293 and C3H10T1/2 cell lines were obtained from the ATCC (Manassas, VA). HEK-293 derivatives 293pTP and RAPA cell lines overexpressing human Ad5 pTP and/or E1 genes were previously described.35, 36 All cell lines were maintained in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS, Hyclone, Logan, UT), 100 units of penicillin, and 100 μg of streptomycin at 37 °C in 5% CO2 as described.37, 38, 39, 40, 41 Unless indicated otherwise, all chemicals were purchased from Sigma–Aldrich (St. Louis, MO), or Thermo Fisher Scientific (Waltham, MA, USA).

Construction and amplification of recombinant adenoviral viruses expressing the 14 types of BMPs

Recombinant adenoviruses were constructed and generated using the AdEasy technology as described.42, 43, 44 Recombinant adenoviruses expressing the 14 types of BMPs were previously characterized and described.23, 31, 32, 45, 46, 47, 48 Briefly, the coding regions of the 14 types of human BMPs were PCR amplified and subcloned into the adenoviral shuttle vector pAdTrack-CMV, and used to generate recombinant adenoviral plasmids in bacterial BJ5183 cells with the Ad5 backbone vector, resulting in pAd-BMPs, which were subsequently used to generate recombinant adenoviruses in HEK-293. The resultant Ad-BMPs adenoviruses were further amplified in 293pTP or RAPA cells.32, 33, 49 These Ad-BMPs also co-express enhanced green fluorescent protein (eGFP) as a tracking marker of viral infection efficiency. The analogous Ad-GFP was used as a mock virus control.50, 51, 52 For all adenoviral infections, polybrene (4–8 μg/ml) was added to enhance infection efficiency as previously described. Detailed information about vector constructions is available upon request.

Total RNA isolation and microarray hybridizations

Exponentially growing, subconfluent C3H10T1/2 cells seeded in 75 cm2 cell culture flasks were infected with a pre-determined optimal titer of Ad-BMPs or Ad-GFP, and maintained in complete medium supplemented with 0.5% FBS. At 30 h after infection, total RNA was isolated using the RNAgent Total RNA Isolation kit (Promega, Madison, WI) by following the manufacturer's instructions.

Microarray analysis was carried out as previously described.26, 27, 28, 29, 31, 54, 55 Briefly, the purity and integrity of total RNA was verified using an Agilent 2100 Bioanalyzer and a GeneSpec III. Fully characterized RNA samples were then used for target preparation by following the Affymetrix GeneChip Expression Analysis Manual (Santa Clara, CA), and subjected to hybridizations to Affymetrix Mouse Genome 430A 2.0 Array GeneChips, which were carried out at the Functional Genomics Facility of The University of Chicago.

Normalization and filtering of the microarray dataset

The acquisition and initial quantification of array images were first examined with the Affymetrix Microarray Suite Version 5.0 (MAS 5.0) with the default parameters (Alpha 1: 0.04; Alpha 2: 0.06; Tau: 0.015; Global scaling target signal: 500). The quality of hybridizations was assessed by examining the MAS 5.0 Report file for housekeeping gene hybridization, Spike control hybridization, percentage of genes called present, 5′ to 3′ ratio, signal to background ratio and scale factor ratio, and then by using DNA-Chip Analyzer (dCHIP) analysis to eliminate regional image contamination and/or sample contamination.

The acquired microarray raw data were then filtered and normalized to remove noise through the following two-step strategy. The first step was to remove the genes in all samples with ≤100 intensity units. The rationale for choosing 100 as the first-step cutoff intensity was based on our observation that the final concentration of spike control BioB in the hybridization mix was 1.5 pM, which was equivalent to 1–3 RNA molecules per cell, but the signal intensity of BioB was normally above 100 when the global scaling target signal was set as 500 (Affymetrix Microarray Suit default setting). Over 99% of the genes with signal intensity ≤100 were called absent by MAS 5.0. The second step was to remove the genes that receive an “absent” call for all chips. Thus, such data filtration ensured that only genes that were considered as significantly present at least in one of the samples were used for further analysis. All data were scaled to a target signal of 500 and therefore they were comparable among samples. In our studies, each sample was normalized to medium signal intensity prior to further analysis. The full list of all transcripts analyzed in this microarray study is shown in Supplemental Table 1.

Microarray data analysis and hierarchical clustering analysis

After normalized and filtering, the microarray data were subjected to SAM analysis. The hierarchical clustering analysis was carried out by using the dCHIP 1.3 software. Specifically, the DNA-Chip Analyzer (dCHIP) was used to analyze the *.CEL files obtained from GCOS (GeneChip Operating Software, Affymetrix). A PM/MM difference model was used to estimate gene expression level and the invariant set approach to normalize data. The thresholds for selecting significant genes were set at a relative difference > 2-fold and absolute difference > 100 signal intensity units at p < 0.05. About 519 genes met the three criteria and considered as significant changes. The dCHIP Clustering Analysis was further carried out to develop the gene clustering tree to visually display high similarity in their standardized expression across all tested samples. Such model-based approach allowed probe-level analysis and facilitated automatic probe selection in the analysis stage to reduce errors caused by outliers, cross-hybridizing probes, and image contamination.

For higher level analysis, we first identified a list of variable genes across pre-defined 3 groups (BMP-2, 6 and 9 as osteogenic group, GFP as control group, and the rest BMPs as unknown group) by ANOVA filtration at p < 0.05. We used this gene list for clustering analysis and Principal Component Analysis (PCA). The dCHIP default setting was used for both analyses. A full list of the genes selected for clustering analysis is shown in Supplemental Table 2.

RNA isolation and quantitative real-time PCR (qRCR) analysis

Exponentially growing C3H10T1/2 cells were infected with a pre-determined optimal titer of Ad-BMPs or Ad-GFP, and maintained in complete medium supplemented with 0.5% FBS. At 30 h after infection, total RNA was isolated using the TRIZOL Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions and subjected to reverse transcription reactions using hexamer and M-MuLV Reverse Transcriptase (New England Biolabs, Ipswich, MA, USA) as previously described.58, 59, 60, 61, 62 The cDNA products were used as PCR templates. The qPCR primers for mouse Smad6 and Smad7 were designed by using the Primer3 Plus program: Smad6, 5′-ATC ACC TCC TGC CCC TGT-3′ and 5′-CTG GGG TGG TGT CTC TGG-3′; and Smad7, 5′-AAG ATC GGC TGT GGC ATC-3′ and 5′-CCA ACA GCG TCC TGG AGT-3′. The qPCR was carried out using our previously optimized TqPCR protocol. Briefly, the 2× SYBR Green qPCR reactions (Bimake, Houston, TX) were set up according to manufacturer's instructions. The cycling program was modified by incorporating 4 cycles of touchdown steps prior to the regular cycling program. Gapdh was used as a reference gene. Each reaction condition was done in triplicate. All sample values were normalized to Gapdh expression by using the 2−ΔΔCt method. The specificity of each qPCR reaction was verified by melting curve analysis and further confirmed by resolving the PCR products on 1.5% agarose gels.

MEV heatmap analysis

To effectively visualize the effects on the 14 BMPs on the expression of downstream target genes and/or crosstalk with other signaling pathways, we employed the cloud-based application WebMev platform for microarray data visualization. Specially, the relative fold changes for a given gene were calculated by dividing the microarray signal units from the BMP groups with that from the GFP control's, and subsequently used to run the MeV analysis for individual genes. The relative expression levels of the genes selected for heatmap analysis are shown in Supplemental Figure 1.

Results

The 14 BMPs exhibit significantly different transcription regulatory effects on I-Smads in MSCs

To explore whether the primary structures are related to their abilities to regulate gene expression, we performed a phylogenetic analysis of the 14 types of human BMPs based on their amino acid sequences (Fig. 1A), and found that several BMPs shared close phylogenetic similarities and formed subclusters, such as BMP2 and BMP4, BMP5 and BMP7, BMP10 and BMP9, BMP13 and BMP14, and BMP3 and BMP15. Nonetheless, such phylogenetic relationships may not necessarily correlate with their biological function relationships.

Figure 1

Phylogenetic analysis of the BMP family and the impact of different BMPs on I-Smad expression in MSCs. (A) Phylogenetic analysis of the 14 types of human BMPs based on their primary amino acid sequences. (B) Effect of individual BMPs on Smad6 and Smad7 expression in MSCs. Subconfluent C3H10T1/2 cells were infected with adenoviruses expressing the 14 types of BMPs or GFP control. At 36 h post infection, total RNA was isolated and subjected to qPCR analysis of Smad6 and Smad7 expression. All qPCR reactions were done in triplicate. Gapdh was used as the reference gene. The % change was calculated by comparing the relative expression in BMP-stimulated samples with that of GFP control group.

Since it is well established that inhibitory Smads (or I-Smads) Smad6 and Smad7 are important feedback inhibitors of the BMP signaling pathways, it's conceivable that expression levels of these I-Smads may reflect individual BMPs' ability to activate the BMP signaling cascade. When the MSCs were stimulated with various BMPs, we found that BMP7 was shown to induce the highest expression levels of Smad6 and Smad7, while BMP8, BMP9, BMP2, BMP3, BMP4, BMP6, BMP10, and BMP12 induced modest to high expression levels of Smad6 and/or Smad7 (Fig. 1B). On the contrary, BMP13 was shown to suppress the expression of Smad6 and Smad7, while BMP14, BMP5, BMP11, and BMP15 were shown to weakly induce the expression of Smad6 and Smad7 in MSCs (Fig. 1B). These results suggest that the 14 BMPs may exhibit drastically distinct abilities to regulate downstream target genes through either Smad-dependent and/or Smad-independent pathways in MSCs.

The 14 BMPs display distinct transcriptomic landscapes, which can be classified into three subclusters in MSCs

We sought to comprehensively profile the gene expression patterns in MSCs upon the stimulation of the 14 types of BMPs. In order to obtain the expression profiles of BMP-regulated immediate early target genes, the exponentially growing MSCs were infected with the Ad-BMPs and control Ad-GFP for 30 h at a low FBS growth condition and harvested total RNAs for microarray analysis. After the acquired dataset was normalized, filtered and processed with dCHIP analysis under high stringencies, 519 significant genes were identified and subjected to hierarchical clustering analysis. The clustering results indicate that, based on the transcriptomic patterns in MSCs, the 14 BMPs could be roughly divided into three subclusters: Cluster 1 that contains BMP9, BMP2, BMP6, BMP7 and BMP4 (Fig. 2A); Cluster 2 that contains BMP11, BMP14, BMP13, BMP15, BMP12, and BMP5 (Fig. 2B); and Cluster 3 that contains BMP10, BMP8 and BMP3 (Fig. 2C).

Figure 2

Hierarchical clustering analysis of the transcripts regulated by 14 BMPs in MSCs. The microarray raw data were first normalized and filtered with MAS 5.0. A list of 519 genes was identified using dCHIP data filtration default settings and was used for hierarchical clustering analysis by dCHIP. Three subclusters are identified: an osteo/chondrogenic/adipogenic cluster (A), a tenogenic cluster (B), and BMP3 cluster (C). The expression level matrix is shown in a log ratio representing normalized values from −3 (blue, below the mean) to +3 (red, above the mean). The mean (0 value) is represented by the white color.

Such clustering results are interesting as there have been reports about the functional similarities of the BMPs within the clusters. For example, all five BMPs in Cluster 1 have been shown to induce osteogenic, chondrogenic and adipogenic differentiation of MSCs although BMP9 is the most potent osteogenic factor, followed by BMP2 and BMP6.14, 17, 20, 23, 24, 25 Thus, Cluster 1 should be considered as an osteo/chondrogenic/adipogenic cluster. For the Cluster 2, BMP13, BMP14 and/or BMP15 were shown to promote tendon/ligament injury repair and significantly improve biomechanical characteristics of injured tendon,22, 66 so this cluster could be tenogenic although the biological functions for BMP11, BMP12 and BMP5 in MSCs or tendon/ligament tissues are not well understood. In the Cluster 3, BMP3 has been well-established as a negative regulator of bone formation,67, 68 although the biological functions of BMP10 and BMP8 in MSCs are currently unclear.

The Principal Component Analysis (PCA) further illustrates the essential transcriptomic or functional relationships among the 14 BMPs (Fig. 3). The three potent osteogenic BMPs (BMP9, BMP2 and BMP6) are more closely related to each other, while the tenogenic BMPs (BMP13, BMP14, and to lesser extent, BMP15) are also shown to closely related. On the contrary, as the only known negative regulator of osteogenic differentiation, BMP3 is distantly related to other BMPs. Thus, to the best of our knowledge, our results are the first to depict a comprehensive transcriptomic landscape of the 14 BMPs in MSCs.

Figure 3

Principal component analysis (PCA) of the transcriptomic correlations regulated by the 14 BMPs in MSCs. The normalized microarray dataset was further filtered by ANOVA at p < 0.005, and the resulted gene list was used for PCA. The color codes indicate closely correlated transcriptomic profiles.

The 14 BMPs exhibit different transcription regulatory effects on the known important mediators of BMP signaling in MSCs

Similar to the qPCR results shown in Fig. 1B,C, the 14 BMPs exhibited distinctly different abilities to induce Smad6 and Smad7 expression in MSCs. BMP2, BMP4, BMP6, BMP7 and BMP9 were shown to up-regulate Smad6 and Smad7 most effectively, while BMP11, BMP12, BMP13, BMP14 and BMP15 displayed weak abilities to induce Smad6 and Smad7 expression (Fig. 4A). Furthermore, BMP2 and BMP9, to lesser extents, BMP3, BMP5, BMP8 and/or BMP10, were shown to modestly induce naturally occurring BMP antagonists follistatin (Fst) and noggin, and regulatory protein latent TGFβ-binding protein 2 (Ltbp2) (Fig. 4A).

Figure 4

Transcriptomic landscape of known BMP downstream targets regulated by the 14 BMPs in MSCs. (A) Microarray-based expression of the feedback inhibitors of the TGFβ/BMP signaling pathway. (B) Microarray-based expression of the known targets of the TGFβ/BMP signaling pathway. The heatmaps were generated by using the web-based MeV (Multiple Experiment Viewer) data visualization program. Fst, follistatin; Ltbp2, latent TGFβ-binding protein 2; Gh, growth hormone, Ctgf, connective tissue growth factor.

We also analyzed the effects of the 14 BMPs on the expression of the currently identified downstream targets of BMP signaling, such as Id2, Id3, Hey1, Hif1α, growth hormone (Gh), and connective tissue growth factor (Ctgf).26, 27, 28, 29, 31, 34 We found that BMP2, BMP4, BMP9, and to lesser extent, BMP6, BMP7, effectively induced the expression of these target genes, although BMP3, BMP5, BMP8, BMP12 and BMP13 were shown to induce some of these target genes (Fig. 4B). Collectively, the above results indicate that the osteogenic BMPs are in general capable of inducing Smad6/Smad7 and/or the known BMP downstream target genes.

Based on the fold of changes and detection p-values, we further identified the top 50 transcripts that were up-regulated by osteogenic BMPs BMP2, BMP6 and BMP9 (Table 1), and the top 50 down-regulated transcripts by BMP2, BMP6 and BMP9 (Table 2). As BMP4 and BMP7 also have certain osteogenic activity, we included their expression levels in these tables. AS indicated in Table 1, the previously identified BMP9 target gene Hey1 is one of the top regulated genes by all osteogenic BMPs, although many target genes on the list have not been characterized. On the other hand, most of the top 50 down-regulated transcripts have yet to be fully characterized, including several members of the heat shock protein family (Table 2). Thus, further studies on the osteogenic BMP-downregulated genes are highly warranted.

Table 1

List of the top 50 transcripts up-regulated by BMP2, BMP6 and BMP9 in MSCs.

Gene	Accession no	BMP2	BMP6	BMP9	BMP4	BMP7	Ave fold change
Hey1	NM_010423.1	17.52	19.83	16.79	20.02	14.15	18.05
Rbm5	NM_148930.3	19.55	12.84	20.37	21.45	18.85	17.59
Mob1a	NM_145571.2	16.59	3.32	17.63	15.87	31.30	12.51
Epb4.1l2	NM_001199265.1	8.09	14.82	7.99	9.02	16.15	10.30
Ahcyl1	NM_145542.3	12.49	2.62	11.95	9.17	15.86	9.02
Bzw1	NM_025824.3	9.66	5.39	8.78	9.19	17.97	7.94
EST	BI696323	7.88	4.15	11.69	13.03	13.42	7.91
Efnb2	NM_010111.5	12.41	5.18	5.77	14.06	3.16	7.79
Il6st	NM_010560.3	8.39	2.68	9.69	6.51	4.91	6.92
Msn	NM_010833.1	6.89	7.06	6.15	7.69	9.36	6.70
Eif2s3x	NM_012010.1	8.34	3.45	7.22	7.78	11.42	6.33
Osmr	NM_011019.3	8.14	2.06	8.34	8.40	2.98	6.18
Gdi2	NM_008112.4	6.50	5.21	6.80	6.37	6.82	6.17
Sptbn1	NM_175836.2	6.35	7.63	4.32	4.92	11.57	6.10
Ghr	NM_010284.2	8.12	2.37	7.72	8.14	4.15	6.07
Angptl2	NM_011923.1	8.45	4.34	5.36	6.56	5.21	6.05
Vcl	NM_009502.1	7.27	0.95	8.87	7.91	10.64	5.69
Golph3	NM_025673.1	5.17	5.81	5.64	6.16	9.36	5.54
Fzd1	NM_021457.1	6.27	3.11	7.16	8.35	3.54	5.51
EST	BM239368	3.95	5.51	6.86	8.32	8.43	5.44
Chp1	NM_019769.1	6.20	4.40	5.70	7.72	10.30	5.43
Ptgs2	NM_011198.3	4.98	5.86	5.36	5.47	20.03	5.40
Dynlt3	NM_025975.2	5.70	4.38	6.01	5.47	6.09	5.36
Smad6	NM_008542.3	4.88	4.03	7.12	4.69	5.74	5.35
Adam17	NM_009615.5	5.66	3.76	6.13	8.39	7.70	5.18
Trim35	NM_029979.3	6.28	2.44	6.73	7.16	11.44	5.15
Sdc1	NM_011519.2	6.89	2.47	6.06	5.36	3.13	5.14
Rad21	NM_009009.4	4.73	4.29	6.38	4.80	5.25	5.13
Myo1c	NM_008659.1	4.94	6.09	4.07	5.08	9.37	5.03
Tpp2	NM_009418.2	5.10	5.14	4.84	4.13	8.00	5.03
Slc12a2	NM_009194.3	6.11	2.00	6.75	3.66	2.54	4.95
Serinc1	NM_019760.3	4.97	3.71	5.97	5.78	6.08	4.88
Cd44	NM_009851.2	5.32	4.37	4.93	5.62	6.23	4.87
Nr3c1	NM_008173.1	5.30	1.24	7.97	7.28	7.32	4.84
Rnf4	NM_011278.4	6.19	1.63	6.67	6.48	9.52	4.83
Ncoa4	NM_019744.1	4.53	4.92	4.99	5.03	7.06	4.81
Serp1	NM_030685.3	5.34	4.34	4.74	2.32	7.42	4.80
G6pd2	NM_019468.1	4.15	6.46	3.75	5.18	6.05	4.79
Thbs1	NM_011580.3	6.16	3.83	4.37	7.08	6.34	4.79
Rab18	NM_011225.1	4.83	3.91	5.62	5.13	5.47	4.78
Tmod3	NM_016963.2	4.81	4.52	4.88	5.16	6.54	4.74
Ahcyl1	NM_145542.3	4.23	5.08	4.77	3.43	5.98	4.70
Caprin1	NM_016739.3	5.85	1.47	6.64	8.88	9.28	4.65
Tmed2	NM_019770.1	4.70	3.97	5.09	5.06	7.66	4.59
Slc1a3	NM_148938.3	5.29	1.53	6.81	5.82	1.94	4.54
Sec61a1	AF145253.1	5.54	3.20	4.86	3.42	7.33	4.53
Adam10	NM_007399.3	5.21	3.82	4.28	5.38	4.66	4.44
EST	BB648848	5.13	1.47	6.68	9.15	8.05	4.43
Timp2	NM_011594.1	5.92	1.60	5.68	5.19	4.46	4.40
Lims1l	NM_026148.1	4.74	3.12	5.20	5.86	7.95	4.35

Table 2

List of the top 50 transcripts down-regulated by BMP2, BMP6 and BMP9 in MSCs.

Gene	Accession no.	BMP2	BMP6	BMP9	BMP4	BMP7	Ave fold change
Hspb1	NM_013560.2	38.87	246.46	15.08	7.99	57.41	45.25
Hsp25	NM_013560.1	4.06	6.66	5.28	2.86	4.59	15.84
Tna	NM_011606.1	4.26	12.75	3.43	2.40	10.66	10.30
Otor	NM_020595.2	5.55	7.63	3.07	2.50	6.67	9.21
EST	AI391218	5.26	5.48	3.42	2.61	6.35	10.26
Apod	NM_007470.1	4.91	8.68	2.79	2.52	7.85	8.37
Hspa1a	NM_010479.2	3.22	4.13	5.43	2.57	0.21	16.30
Tm4sf2	NM_019634.1	3.67	18.84	1.98	1.91	12.38	5.94
Angptl4	NM_020581.1	3.36	27.02	1.99	2.98	39.80	5.97
Arhn	NM_009708.1	2.91	4.29	3.68	4.34	1.74	11.03
Cd248	NM_054042.1	4.50	2.12	6.12	3.11	3.21	18.35
Ssbp2	NM_024186.5	3.33	5.36	2.56	2.95	4.66	7.67
Psmb10	NM_013640.1	3.16	7.03	2.26	3.14	6.14	6.78
Sugp2	NM_001168290.1	3.13	3.22	3.41	4.13	3.40	10.24
Palm	NM_023128.4	3.00	2.94	3.88	2.45	3.73	11.64
Ugcg	NM_011673.3	3.59	4.71	2.18	2.78	4.17	6.54
Vat1	NM_012037.2	3.70	2.26	4.20	6.18	3.64	12.59
Mt1	NM_013602.1	2.56	13.25	1.99	3.24	17.17	5.97
Vamp2	NM_009497.3	3.19	2.53	3.74	2.45	2.85	11.22
Synpo	NM_177340.2	2.36	8.05	2.14	3.02	11.72	6.42
Ppnr	NM_012022.1	4.24	2.70	2.41	3.34	4.25	7.22
Igf1	NM_010512.4	2.59	4.69	2.35	1.32	5.16	7.04
Tle2	NM_019725.2	2.95	5.52	1.91	1.65	4.72	5.74
Safb2	NM_001029979.2	3.85	1.90	3.84	3.65	4.22	11.51
Abi3bp	NM_178790.4	3.14	15.04	1.50	1.45	7.74	4.50
Actl6a	NM_019673.2	3.05	2.40	3.23	2.03	1.75	9.70
Ifitm1	NM_026820.3	2.44	12.37	1.75	2.85	5.86	5.26
Sepw1	NM_009156.2	2.64	2.93	2.88	2.92	4.30	8.64
Psmc4	NM_011874.1	3.26	3.26	2.19	3.23	1.55	6.57
Fxyd1	NM_019503.1	2.44	5.98	1.98	1.43	4.73	5.94
Nol12	NM_133800.3	2.88	2.31	3.26	1.98	2.14	9.78
Pafah1b1	NM_013625.4	3.11	7.15	1.60	2.30	2.41	4.80
Clu	NM_013492.2	2.71	9.68	1.61	2.08	9.25	4.83
Sfxn3	NM_053197.1	2.14	3.52	2.87	2.45	1.72	8.61
Cnp	NM_001146318.1	2.73	2.58	2.89	1.55	2.09	8.67
Egr1	NM_007913.1	4.73	2.76	1.90	2.66	2.55	5.70
Ets2	NM_011809.3	2.50	3.21	2.50	2.62	2.74	7.51
H2-Ke6	NM_013543.2	2.54	2.73	2.72	2.50	2.73	8.16
Fkbp8	NM_010223.1	2.25	3.63	2.44	2.91	4.25	7.33
Il11ra1	NM_010549.3	2.52	3.89	2.08	2.64	2.53	6.24
Ephb3	NM_010143.1	2.31	3.35	2.46	2.03	5.69	7.37
Sepp1	NM_009155.3	2.66	11.30	1.48	1.60	11.60	4.44
Hagh	NM_001159626.1	2.23	3.00	2.77	3.49	1.76	8.32
Cdc42ep5	NM_021454.1	1.82	6.32	2.14	2.03	5.78	6.43
Tap2	NM_011530.3	2.19	5.85	1.82	2.54	4.25	5.45
Ssbp4	NM_133772.1	2.26	2.64	2.76	1.85	2.14	8.27
C3	NM_009778.2	2.40	38.87	1.35	1.50	62.50	4.04
EST	AV110855	2.35	1.78	4.87	3.56	7.76	14.62
Thap7	NM_026909.2	2.51	2.10	3.14	2.75	1.75	9.42
Fntb	NM_145927.2	2.23	3.09	2.34	3.01	2.17	7.02

We further analyzed the differentially regulated transcripts by tenogenic BMPs, e.g., BMP13, BMP14 and BMP15. As shown in Table 3, most of the top 50 up-regulated genes and EST transcripts are yet to be fully studied. Interestingly, one of the top 50 down-regulated genes/transcripts by the three tenogenic BMPs also belongs to the heat shock protein family (Table 4) although the biological functions for the vast majority of the top down-regulated genes/transcripts need to be investigated, especially in the context of tenocytes and/or ligament tissue.

Table 3

List of the top 50 transcripts up-regulated by BMP13, BMP14 and BMP15 in MSCs.

Gene	Accession no.	Fold change by BMP13	Fold change by BMP14	Fold change by BMP15	Ave detection p-value	Ave fold increase
Rtp4	NM_023386.1	2.16	11.34	2.18	0.0022	5.23
Epb4.1l2	NM_001199265.1	1.89	4.78	6.76	0.0029	4.48
EST	BI696323	1.52	3.63	6.90	0.0007	4.01
Vcl	NM_009502.1	0.68	3.17	6.98	0.0217	3.61
Il6st	NM_010560.3	1.62	3.81	5.08	0.0012	3.50
Osmr	NM_011019.3	1.77	3.50	5.00	0.0161	3.42
Runx1	NM_009821.1	2.37	3.01	4.60	0.0024	3.33
Eif2s3x	NM_012010.1	1.74	2.59	5.65	0.0024	3.32
Fzd1	NM_021457.1	1.19	3.30	5.32	0.0017	3.27
Sgpp1	NM_030750.1	1.82	2.84	4.81	0.0063	3.15
Vps35	NM_022997.4	1.23	2.78	5.23	0.0269	3.08
Gdi2	NM_008112.4	1.24	3.05	4.92	0.0039	3.07
Rnf4	NM_011278.4	1.71	2.62	4.79	0.0024	3.04
Nedd4a	NM_010890.1	0.96	2.93	4.87	0.0017	2.92
Hnrnpr	NM_028871.1	1.24	2.48	4.92	0.0029	2.88
Rad21	NM_009009.4	1.06	3.03	4.52	0.0300	2.87
Kpnb1	NM_008379.3	1.27	1.97	5.27	0.0017	2.84
Ifit1	NM_008331.1	2.04	5.44	0.98	0.0007	2.82
Atf2	NM_001025093.1	1.92	3.15	3.37	0.0007	2.81
Nr3c1	NM_008173.1	0.98	2.55	4.85	0.0034	2.80
Ccl9	NM_011338.2	2.57	2.97	2.63	0.0044	2.72
Mgea5	NM_023799.1	1.41	2.90	3.78	0.0078	2.70
Sec61a1	AF145253.1	1.06	2.63	4.39	0.0100	2.70
P4ha1	NM_011030.2	1.40	2.64	3.99	0.0051	2.68
EST	BM239368	1.30	2.75	3.96	0.0012	2.67
Ahcyl1	NM_145542.3	1.46	2.97	3.52	0.0007	2.65
Trim35	NM_029979.3	1.17	2.45	4.29	0.0330	2.64
Zfp313	NM_030743.1	1.37	2.84	3.48	0.0095	2.56
Adam10	NM_007399.3	1.01	2.95	3.66	0.0056	2.54
Dynlt3	NM_025975.2	0.89	2.70	4.00	0.0249	2.53
Sh3bgrl	NM_019989.1	1.14	2.77	3.67	0.0007	2.52
Srfbp1	NM_026040.1	2.52	2.19	2.80	0.0425	2.50
Man1a	NM_008548.4	2.41	3.18	1.90	0.0051	2.49
Ccng1	NM_009831.2	1.96	2.45	3.06	0.0007	2.49
Golph3	NM_025673.1	0.85	2.18	4.32	0.0007	2.45
EST	AK005023.1	1.02	1.99	4.24	0.0024	2.42
Saa3	NM_011315.1	3.07	2.01	2.16	0.0022	2.41
EST	BG066678	2.00	3.13	2.10	0.0129	2.41
Timp2	NM_011594.1	0.83	2.48	3.87	0.0007	2.39
Myo1c	NM_008659.1	0.96	2.22	3.94	0.0222	2.38
Scya5	NM_013653.1	1.91	3.32	1.74	0.0034	2.32
Aspn	NM_025711.3	1.75	2.40	2.82	0.0007	2.32
Psmg2	NM_134138.1	1.57	2.97	2.42	0.0012	2.32
Serp1	NM_030685.3	1.43	2.04	3.47	0.0007	2.32
Tnfrsf6	NM_007987.1	2.42	2.17	2.32	0.0198	2.30
Angptl2	NM_011923.1	1.26	2.40	3.24	0.0007	2.30
Caprin1	NM_016739.3	0.96	2.01	3.92	0.0063	2.30
Slc1a3	NM_148938.3	1.65	2.26	2.93	0.0083	2.28
Nbr1	NM_008676.3	1.09	2.28	3.47	0.0012	2.28
Msn	NM_010833.1	1.01	1.91	3.81	0.0034	2.24

Table 4

List of the top 50 transcripts down-regulated by BMP13, BMP14 and BMP15 in MSCs.

Gene	Accession no.	Fold change by BMP13	Fold change by BMP14	Fold change by BMP15	Ave detection p-value	Ave fold decrease
Hspb1	NM_013560.2	9.22	35.21	2.82	0.7710	6.10
Hspa1a	NM_010479.2	5.24	3.95	2.42	0.1130	3.50
EST	BC004738.1	3.07	2.82	2.58	0.4104	2.81
Hmga2	NM_010441.1	2.18	2.52	2.63	0.0012	2.43
Hsp25	NM_013560.1	3.24	2.78	1.71	0.0012	2.39
Fam101b	NM_029658.1	1.98	2.05	3.56	0.0649	2.36
Rpl37a	NM_009084.1	5.51	3.12	1.14	0.0022	2.18
Tmx2	NM_025868.3	1.88	2.22	2.53	0.0801	2.18
Efnb2	NM_010111.5	2.09	4.63	1.45	0.0227	2.17
Tnrc6a	NM_144925.3	2.21	2.45	1.75	0.0029	2.09
Hmgcs1	NM_145942.4	1.89	2.52	1.82	0.0007	2.03
Sugp2	NM_001168290.1	1.58	2.99	1.92	0.0063	2.01
Ehd1	NM_010119.1	1.90	1.99	2.11	0.0356	2.00
Nol12	NM_133800.3	1.63	2.47	2.02	0.0674	1.98
Fntb	NM_145927.2	1.81	2.08	2.02	0.0034	1.96
Slc25a37	NM_026331.3	1.40	2.30	2.49	0.0247	1.93
Zfp777	NM_001081382.1	1.61	1.51	3.65	0.1377	1.93
Fdft1	NM_010191.2	1.93	2.17	1.70	0.0007	1.91
Rbx1	NM_019712.3	2.24	1.62	1.97	0.0017	1.91
Klf13	NM_021366.3	1.56	1.99	2.19	0.0022	1.87
Ebf1	NM_007897.1	2.33	1.77	1.63	0.0007	1.87
Nop56	NM_024193.2	1.58	2.57	1.66	0.0022	1.84
Vat1	NM_012037.2	1.51	1.91	2.22	0.0061	1.83
Actl6a	NM_019673.2	1.55	1.47	2.96	0.0039	1.80
L1cam	NM_008478.3	2.29	2.13	1.30	0.0007	1.79
Ccm2	NM_146014.3	1.55	2.33	1.64	0.1814	1.78
Dag1	NM_010017.3	2.33	1.75	1.46	0.0007	1.78
Lonp1	NM_028782.2	1.62	1.85	1.90	0.0022	1.78
Cnn2	NM_007725.2	1.80	1.97	1.59	0.0007	1.77
Atp13a2	NM_029097.2	2.25	1.39	1.85	0.0083	1.76
Sfxn3	NM_053197.1	1.37	2.01	2.12	0.0007	1.76
Qars	NM_133794.2	1.51	1.80	1.94	0.0007	1.73
Ctdspl	NM_133710.1	1.67	1.45	2.16	0.0012	1.71
Ppt2	NM_019441.4	1.36	2.12	1.80	0.0107	1.70
Plec	NM_011117.2	1.53	2.00	1.62	0.0007	1.69
Tpi1	NM_009415.2	1.58	2.07	1.51	0.0007	1.69
Tpm2	NM_009416.2	1.52	1.59	2.02	0.0007	1.68
Hspa1b	NM_010478.2	2.19	1.88	1.24	0.0007	1.67
Hn1	NM_008258.1	1.62	1.59	1.79	0.0034	1.66
Fam111a	NM_026640.1	1.79	1.63	1.57	0.0007	1.66
Ptms	NM_026988.2	1.67	1.80	1.53	0.0007	1.66
Palm	NM_023128.4	1.65	1.33	2.21	0.0081	1.66
Bptf	NM_176850.2	1.75	1.76	1.48	0.0159	1.65
Col6a2	NM_146007.2	1.39	1.77	1.87	0.0007	1.65
Cox7a2l	NM_009187.1	1.57	1.43	2.05	0.0012	1.65
Cnp	NM_001146318.1	1.58	2.15	1.36	0.0034	1.64
Cdk2ap1	NM_013812.2	1.46	1.65	1.83	0.0042	1.64
Crlz1	NM_023054.1	1.93	1.08	2.53	0.0034	1.63
Smpdl3b	NM_133888.1	1.44	1.64	1.83	0.0286	1.62
Cd151	NM_009842.3	1.35	1.95	1.65	0.0007	1.61

The 14 BMPs differentially regulate lineage-specific regulators and markers in MSCs

We next analyzed how the 14 BMPs regulated lineage-specific differentiation of MSCs. For osteogenic lineage, the microarray analysis indicated that BMP2, BMP4, BMP9, BMP11 and BMP15 induced Runx2 expression, while BMP2, BMP3, BMP4, BMP7 and BMP11 up-regulated Sost expression (Fig. 5). Interestingly, BMP6 and BMP7 did not significantly induce Runx2 expression. Chondrogenic markers Sox9, Sox5 and Sox6 were upregulated by BMP4, and/or BMP11, BMP6, BMP7, BMP9, BMP13 or BMP14 (Fig. 5). BMP8, BMP9 and BMP12 were shown to induce Pparγ expression while BMP5, BMP8, BMP9, BMP10, BMP13, and BMP14 upregulated Lpl expression (Fig. 5). For myogenic lineage, BMP2, BMP4, BMP6, BMP7, and BMP8 were shown to induce MyoD expression (Fig. 5). Lastly, the tenogenic marker Pax3 was induced by BMP5, BMP10, BMP12, and BMP15, while Scx expression was repressed by BMP2 through BMP10, and BMP14 (Fig. 5). Interestingly, some of the lineage markers, known to be regulated by certain BMPs, were not significantly impacted in this microarray analysis. One possible explanation is that the microarray analysis was conducted in the MSCs stimulated by BMPs at the early time point (i.e., 30 h post Ad-BMP infection). Based our prior experience, many of the lineage-specific regulators and markers are usually robustly induced at 48–72 h after Ad-BMP infection.30, 45, 69, 70 Nonetheless, the microarray analysis offered a comprehensive view about the effects of BMPs on lineage-specific differentiation of MSCs.

Figure 5

Transcriptomic effect of the 14 BMPs on lineage-specific regulators and markers in MSCs. Microarray-based expression analysis of the osteogenic (Runx2 and Sost), chondrogenic (Sox9, Sox5 and Sox6), adipogenic (Pparγ and Lpl), myogenic (MyoD1), and tenogenic (Scx and Pax3) regulators and/or markers affected by the 14 BMPs in MSCs. The heatmaps were generated by using the web-based MeV program. The genes with an increased overall expression by some BMPs are labeled in red, and the ones with decreased expression labeled in blue, while those without significant changes are labeled in black, when compared with that of the GFP control's. Sost, sclerostin; Pparg, peroxisome proliferator-activated receptor gamma; Lpl, lipoprotein lipase; Scx, scleraxis BHLH transcription factor.

The 14 BMPs differentially affect MAPK signaling as a potential Smad-independent (non-canonical) signaling pathway in MSCs

It has been well documented that BMPs may transduce signaling through Smad-independent, non-canonical signaling pathways, mostly through the MAP kinase pathway.71, 72 We found that at least 9 of 17 members of the detected MAPK family were significantly up-regulated by several BMPs, such as BMP2, BMP3, BMP4, BMP7, BMP11, and to much lesser extent, BMP5, BMP8, BMP9, BMP10, BMP12, BMP13, BMP14 and BMP15 (Fig. 6). Mapk6, Mapk10, and Mapk14 were shown most significant up-regulation, while BMP3, BMP4, BMP7, and BMP11 were shown to upregulate most members of MAPK family (Fig. 6). Conversely, Mapk4k4 expression was repressed by most BMPs except BMP13 and BMP14, while Mapk3k11 expression was repressed by 10 of the 14 types of BMPs (Fig. 6). Furthermore, BMP13, BMP14 and BMP15 were shown to exert the least effect on MAPK expression (Fig. 6). Taken together, these results indicate that most BMPs can regulate MAPK expression and thus transduce signaling through the MAPK signaling pathway.

Figure 6

Transcriptomic effect of the 14 BMPs on the MAP kinase pathway in MSCs. Microarray-based expression analysis of the members of the MAP kinase pathway regulated by the 14 BMPs in MSCs. The heatmaps were generated by using the web-based MeV program. The genes with an increased overall expression by some BMPs are labeled in red, and the ones with decreased expression labeled in blue, while those without significant changes are labeled in black, when compared with that of the GFP control's. Significantly up-regulated transcripts are labeled in red.

The 14 BMPs exhibit different crosstalk capabilities with the Wnt signaling pathway in MSCs

The BMP signaling pathway can crosstalk with several major signaling pathways.6, 17, 73, 74 We first examined whether BMPs regulate the expression of Wnts and/or the important components of the Wnt signaling pathway. We found that 10 Wnts (e.g., Wnt1, Wnt3a, Wnt5b, Wnt7a, Wnt7b, Wnt8b, Wnt11, Wnt14, Wnt14b, and Wnt16) were up-regulated by various BMPs, while Wnt8a expression was repressed by most BMPs except BMP3, BMP7, BMP11, and BMP12 (Fig. 7A). Among the 14 BMPs, BMP4 and BMP7, to lesser extent, BMP2, BMP3, BMP5 and BMP6, were shown to up-regulate the expression of the most numbers of Wnts, while BMP9 and BMP12 only up-regulated one Wnt expression (Fig. 7A).

Figure 7

Transcriptomic effect of the 14 BMPs on the Wnt signaling pathway in MSCs. Microarray-based expression analysis of the members of the Wnt signaling pathway regulated by the 14 BMPs in MSCs. The heatmaps were generated by using the web-based MeV program. The genes with an increased overall expression by some BMPs are labeled in red, and the ones with decreased expression labeled in blue, while those without significant changes are labeled in black, when compared with that of the GFP control's. Significantly up-regulated transcripts are labeled in red, while down-regulated genes are labeled in blue.

We further analyzed the effects of BMP signaling on the expression of Wnt receptors, and other critical components of Wnt signaling pathway, and found that most of the detected Fzd receptors (e.g., Fzd1, Fzd3, Fzd6, and Fzd7) were up-regulated by various BMPs, especially by BMP2, BMP4, BMP8, BMP9, and BMP11(Fig. 7B). Furthermore, we found that β-catenin expression was up-regulated by BMP2, BMP8, BMP9, and BMP15, while BMP2, BMP8, BMP9 and BMP14 were shown to up-regulate Lef1 expression (Fig. 7B). Surprisingly, the expression of Tcf3 and Tcf4 was repressed by most BMPs, except BMP12, BMP13, BMP14 and BMP15 (Fig. 7B). Taken together, these results indicate that at least one of the molecular mechanisms underlying BMP-Wnt crosstalk may occur at the expression levels of Wnt ligands, Fzd receptors, β-catenin and Lef1 in MSCs.

The 14 BMPs differentially regulate the components of Notch signaling in MSCs

We also analyzed the effect of BMPs on the expression of Notch receptors and ligands. As shown in Fig. 8, Notch2, Notch4, Dll1, Dll3, Dll4 and Jag1 were up-regulated by various BMPs, especially BMP2, BMP3, BMP4, BMP5, BMP7, and BMP10. Interestingly, Notch1 expression was repressed by most BMPs except BMP5, BMP6, BMP7 and BMP14 (Fig. 8). Furthermore, BMP12 and BMP15, to lesser extend BMP6 and BMP13, did not significantly affect the expression of Notch receptors and ligands (Fig. 8). Interestingly, BMP9 was shown to upregulate the expression of Notch2 and Dll1 (Fig. 8). Nonetheless, these results suggest that the expression of Notch receptors and/or ligands may be effectively regulated in MSCs by some BMPs, such as BMP4, BMP2 and BMP7.

Figure 8

Transcriptomic effect of the 14 BMPs on the Notch signaling pathway in MSCs. Microarray-based expression analysis of the members of the Notch signaling pathway regulated by the 14 BMPs in MSCs. The heatmaps were generated by using the web-based MeV program. The genes with an increased overall expression by some BMPs are labeled in red, and the ones with decreased expression labeled in blue, while those without significant changes are labeled in black, when compared with that of the GFP control's. Significantly up-regulated transcripts are labeled in red, while down-regulated genes are labeled in blue.

The 14 BMPs differentially regulate PI3K/AKT/mTOR signaling axis in MSCs

BMP signaling has been shown to regulate PI3K/AKT/mTOR signaling axis.75, 76, 77, 78, 79 Here, we analyzed the effects of the 14 BMPs on the expression of the major components of PI3K/AKT/mTOR axis. We found that Pik3c2g, Pik3cd, Pik3r1, Akt2 and Akt3 were up-regulated by several BMPs, notably by BMP2, BMP3, BMP4, BMP7, BMP8, BMP9, and BMP11, while BMP13 and BMP14 did not significantly affect the expression level of PI3K/AKT/mTOR (Fig. 9). Interestingly, the expression of Pik3r2 and Akt1 was repressed by most BMPs. Furthermore, mTor expression was only marginally affected by BMP3, BMP5, and BMP15 (Fig. 9). Taken together, these results suggest that most of the 14 types of BMPs may up-regulate Pik3c2g, Pik3cd, Pik3r1, Akt2 and Akt3 but have little effect on Akt1 and mTor expression in MSCs.

Figure 9

Transcriptomic effect of the 14 BMPs on the PI3K/AKT/mTOR pathway in MSCs. Microarray-based expression analysis of the members of the PI3K/AKT/mTOR pathway regulated by the 14 BMPs in MSCs. The heatmaps were generated by using the web-based MeV program. The genes with an increased overall expression by some BMPs are labeled in red, and the ones with decreased expression labeled in blue, while those without significant changes are labeled in black, when compared with that of the GFP control's. Significantly up-regulated transcripts are labeled in red, while down-regulated genes are labeled in blue.

Discussion

BMPs play important roles in regulating the proliferation and lineage-specific differentiation of mesenchymal stem cells. In this study, we conducted a comprehensive transcriptomic profiling of the MSCs stimulated by the 14 types of BMPs at the early stage, which allows us to identify the immediate early transcriptomic changes in MSCs. The feedback inhibition mechanism-mediated through Smad6/7 indicates that several BMPs (e.g., BMP2, BMP3, BMP4, BMP7, BMP8, and BMP9) induce the Smad-mediated signaling pathway more effectively, while BMP5, BMP11, BMP13, BMP14, and to a lesser extent BMP12 and BMP15 do not significantly impact Smad signaling in MSCs. Nonetheless, since BMPs can induce signaling activities through the Smad-independent, non-canonical signaling pathways, mostly through the MAP kinase pathway,71, 72 our results reveal that at least 9 of 17 members of the detected MAPK family were significantly up-regulated by most BMPs, specially BMP2, BMP3, BMP4, BMP7, and BMP11, although BMP13, BMP14 and BMP15 exerted the least effect on MAPK expression.

Transcriptomic clustering analysis indicates that the biological functions of different BMPs can be classified into three large groups and are not always related to their phylogenetic relationships in MSCs. For example, although BMP2 and BMP4, BMP6 and BMP7 are phylogenetically distant from BMP9, their expression profiles classified them as a large subcluster of osteogenic BMPs, which is consistent with the results we previously reported previously in MSCs.17, 20, 23, 24, 25, 45, 73 Conversely, the phylogenetically related BMP12, BMP13, BMP14, to lesser extent BMP15, form a subcluster of tenogenic factors in MSCs.22, 66 As expected, BMP3 exhibits a transcriptomic profile almost opposite to that of BMP9, which is consistent with the bone forming inhibitory role of BMP3 in MSCs. Interestingly, it has been reported that BMP9 and BMP10 may be functionally redundant in vascular endothelial cells.80, 81 Yet, our transcriptomic analysis indicates that BMP9 is distantly related to BMP10 in MSCs. Thus, it is conceivable that BMPs regulate transcriptomic landscapes may be cell-type and/or context-dependent. Nonetheless, our gene profiling analysis identified the lists of the most up-regulated and down-regulated genes by the osteogenic BMPs and tenogenic BMPs, respectively.

A comprehensive transcriptomic analysis also allows us to examine the scopes and extents of the crosstalk between BMPs and other major signaling pathways.6, 17, 73, 74 Our results revealed that at least 10 Wnts (e.g., Wnt1, Wnt3a, Wnt5b, Wnt7a, Wnt7b, Wnt8b, Wnt11, Wnt14, Wnt14b, and Wnt16) were up-regulated by various BMPs (especially BMP2, BMP3, BMP4, BMP5 and BMP6, and BMP7), while Wnt8a expression was repressed by most BMPs. Similarly, we found most of the detected Fzd receptors (e.g., Fzd1, Fzd3, Fzd6, and Fzd7) were up-regulated by various BMPs, especially by BMP2, BMP4, BMP8, BMP9 and BMP11, while β-catenin expression was up-regulated by BMP2, BMP8, BMP9 and BMP15. Furthermore, the Lef1 expression was up-regulated by BMP2, BMP8, BMP9 and BMP14 although the expression of Tcf3 and Tcf4 was repressed by most BMPs except BMP12, BMP13, BMP14 and BMP15.

Accordingly, our results further demonstrated many members of the Notch signaling pathway (such as Notch2, Notch4, Dll1, Dll3, Dll4 and Jag1) were up-regulated by several BMPs, especially BMP2, BMP3, BMP4, BMP5, BMP7 and BMP10, while Notch1 expression was repressed by most BMPs. Nonetheless, several BMPs, such as BMP6, BMP12, BMP13 and BMP15, exerted little effect on the expression of Notch receptors and ligands. Lastly, we analyzed the effects of the 14 BMPs on the expression of the major components of PI3K/AKT/mTOR axis, and found that Pik3c2g, Pik3cd, Pik3r1, Akt2 and Akt3 were up-regulated by several BMPs, most notably by BMP2, BMP3, BMP4, BMP7, BMP8, BMP9, and BMP11. However, our results revealed that most BMPs exerted little effect on Akt1 and mTor expression in MSCs.

It is noteworthy that our transcriptomic profile only reflects the early stage of BMP-regulated gene expression in MSCs as a thorough analysis should include multiple time points post BMP stimulations. This may explain why some of the lineage-specific regulators and/or markers (e.g., Sox9, MyoD, and Scx) were not significantly induced in our microarray analysis. We previously reported that these lineage-regulators were induced in the intermediate stage (usually 2–4 days after BMP stimulation) in MSCs.22, 24, 45, 66, 82 Another limitation of the current study is that our microarray analysis was only carried out a single MSC line. It is conceivable that BMP-regulated transcriptomic landscapes would be more comprehensive if multiple MSC lines were examined. The highly differentially expressed genes should be further validated and functional characterized. With the rapid development of the next-generation sequencing technologies, large scales of genomewide transcriptomic analyses are made possible. Thus, the employment of new genome-based technologies should further expand our understanding of BMP signaling mechanisms in MSCs.

Conflict of interest

The authors declare no conflict of interests.