Epigenetic Plasticity Drives Adipogenic and Osteogenic Differentiation of Marrow-derived Mesenchymal Stem Cells*

Terminal differentiation of multipotent stem cells is achieved through a coordinated cascade of activated transcription factors and epigenetic modifications that drive gene transcription responsible for unique cell fate. Within the mesenchymal lineage, factors such as RUNX2 and PPARγ are indispensable for osteogenesis and adipogenesis, respectively. We therefore investigated genomic binding of transcription factors and accompanying epigenetic modifications that occur during osteogenic and adipogenic differentiation of mouse bone marrow-derived mesenchymal stem cells (MSCs). As assessed by ChIP-sequencing and RNA-sequencing analyses, we found that genes vital for osteogenic identity were linked to RUNX2, C/EBPβ, retinoid X receptor, and vitamin D receptor binding sites, whereas adipocyte differentiation favored PPARγ, retinoid X receptor, C/EBPα, and C/EBPβ binding sites. Epigenetic marks were clear predictors of active differentiation loci as well as enhancer activities and selective gene expression. These marrow-derived MSCs displayed an epigenetic pattern that suggested a default preference for the osteogenic pathway; however, these patterns were rapidly altered near the Adipoq, Cidec, Fabp4, Lipe, Plin1, Pparg, and Cebpa genes during adipogenic differentiation. Surprisingly, we found that these cells also exhibited an epigenetic plasticity that enabled them to trans-differentiate from adipocytes to osteoblasts (and vice versa) after commitment, as assessed by staining, gene expression, and ChIP-quantitative PCR analysis. The osteogenic default pathway may be subverted during pathological conditions, leading to skeletal fragility and increased marrow adiposity during aging, estrogen deficiency, and skeletal unloading. Taken together, our data provide an increased mechanistic understanding of the epigenetic programs necessary for multipotent differentiation of MSCs that may prove beneficial in the development of therapeutic strategies.

The mesenchymal stem cell (MSC) 2 is a multipotential unit with the ability to differentiate along osteogenic, adipogenic, chondrogenic, myogenic, and neurogenic pathways (1). A common source of MSCs is derived from bone marrow. In addition to bone marrow, however, MSCs can also be isolated from adipogenic tissue sources (aMSC or ASC), umbilical cord tissue, pancreas, and most other vascularized tissues (1). Regardless of the source, these cells retain their ability to differentiate into a number of common lineages when appropriately stimulated, displaying cell markers such as CD29, CD44, CD70, CD73, CD90, CD105, and Stro-1, while lacking markers for early hematopoietic lineages (CD34 and others) (2,3). Although the markers for identification of a "true" MSC remain somewhat controversial, multipotential differentiation morphology remains the hallmark of the MSC. These phenotypes are driven by distinct and essential sets of cell signaling pathways and their respective target transcription factors. Adipogenic differentiation relies on the expression of peroxisome proliferator-activated receptor ␥ (PPAR␥), and agonists of the PPAR␥ pathway promote the formation of adipocytes at the expense of osteoblasts (4,5). Osteogenic differentiation, on the other hand, cannot proceed without expression of Runt-related transcription factor 2 (RUNX2) (6,7) and its downstream effects on osterix (OSX, Sp7) expression. Wnt and Hedgehog signaling pathways favor the osteogenic versus the adipogenic lineage (8,9). Likewise, bone morphogenic proteins and TGF␤ direct MSCs toward both osteogenic and chondrogenic phenotypes (10). Interleukins and STAT signaling uniquely promote specific chondrogenic differentiation (11). In addition to chemical signals, MSCs are able to respond to cues from physical interactions with intranuclear actin to up-regulate RUNX2 as well as to repress adipogenesis through down-regulation of C/EBP␤ and PPAR␥ in response to mechanical loading that favors bone (12)(13)(14). In all, the MSC is a versatile cell capable of responding to various signals for varied terminal differentiation programs.
One of the major transcription factors active in preosteogenesis is the vitamin D receptor (VDR), a ligand-activated second-ary activator of MSCs. VDR is able to regulate the expression of genes that are vital to osteoblastic identity in early osteoblastic cells as well as in MSCs through coordination with the actions of RUNX2 and C/EBP␤ (15,16). The VDR has a pervasive genomic footprint in preosteoblastic cells that is reduced as these cells undergo further osteoblastic differentiation, but VDR is still an active regulator of transcription in both differentiation states (15). Likewise, recent efforts have been made to establish a role for the VDR in the adipocyte, where, as in the osteoblast, VDR has a somewhat conflicted role throughout differentiation. Global VDR knock-out in mice results in adipocyte atrophy, and mice with the VDR ablated specifically in adipose tissue show increased mammary epithelial density (17,18). The active metabolite of vitamin D, 1,25(OH) 2 D 3 , has been shown to prevent adipocyte differentiation in vitro and simultaneously up-regulate cytokine and leptin activity in fat tissue (19,20). Despite the specifics of these activities, VDR is present in cells of all mesenchymally derived lineages and remains an important transcription factor in vivo for both osteoblasts and adipocytes.
The regulatory activities of these transcription factors are accompanied (or preceded) by epigenetic modifications to either DNA or to histone protein tails to allow additional transcription factor access to genes vital for lineage determination. The epigenetic marks in osteoblastic cells (15,16,21) and in adipocyte cell culture models (22) have revealed distinct patterns for gene expression in terminally differentiated cell types. Previous research has determined that the histone-modifying enzymes responsible for many of these patterns, including KDM5B, EZH2, and others, play a clear role in skeletal development and differentiation of the MSC through modulation of RUNX2 (23,24). A recent study found distinct open regions of DNase I-hypersensitive DNA in a human temperature-dependent model of osteoblast differentiation that favored genes known for the osteoblastic phenotype (25). Chromatin marks are dynamic in nature and can shift location based upon the stimulus received (26). In fact, the plasticity of histone and chromatin modifications is believed to allow for memory context in the brain (27), circadian rhythm (28,29), and even the potential for trans-differentiation from one cell type to another, as exemplified by smooth muscle cell transitions into macrophages in atherosclerosis (30). Trans-differentiation of mesenchymal cells can be used to create neural stem cells (31), insulin-producing cells (32), auditory cells (33), or other therapeutically relevant cell lines. Their accessibility and availability have made MSCs (stromal or adipose-derived) an attractive cell type for studies in personalized regenerative medicine.
Given the diverse differentiation potential of marrow-derived MSCs as well as their dynamic nature, we sought to investigate the transcription factor cascades as well as histone modification patterns that characterize these differentiating cells. We differentiated MSCs in culture for 7 or 15 days using osteogenic medium and 7 days using adipogenic medium. At these time points, we performed RNA-seq analysis to establish changes in the transcriptome and ChIP-seq analysis to investigate the cistromes for a number of transcription factors (RUNX2, C/EBP␤, C/EBP␣, PPAR␥, VDR, and RXR) as well as several epigenetic marks (H3K9ac, H4K5ac, H3K4me1, H3K4me3, and H3K36me3) that are associated with these transitions. We found that chromatin state profiles for the undifferentiated MSC closely resembled those of osteogenic cells differentiated for 7 or 15 days, indicating a potentially predetermined epigenetic pathway for these cells. These chromatin state profiles were dramatically altered, however, following adipogenic differentiation for 7 days, especially near genes vital for adipogenic function. These profiles accurately predicted the transition and gene output activity at the locus in which they resided. Finally, we found that despite numerous genomic, chromatin, and gene expression changes, MSCs were able to trans-differentiate from an adipogenic to an osteogenic lineage and vice versa. These two cell types displayed the appropriate cell morphology, transcriptomic, and genomic changes of the emerging trans-differentiated phenotype and demonstrated a plasticity in the histone modifications that allowed these transitions to occur in vitro. These observations provide a potential mechanism during pathological conditions that lead to skeletal fragility and increased marrow adiposity as in aging, estrogen deficiency, and in skeletal unloading. These expression patterns and profiles will be beneficial in understanding the mechanisms that underlie the cellular and genomic transitions that occur during differentiation and may lead to insights that will be required for the successful development of therapies or drug discovery models in the future.

Marrow-derived MSCs Differentiate into Mature Adipocytes and Osteoblasts-Marrow-derived
MSCs present an interesting model of differentiation because they are capable of transitioning into a wide variety of cell types, most notably osteoblasts, adipocytes, and chondrocytes. In the marrow space, stromal cells differentiate into these cell types to maintain bone integrity and health. We chose to investigate the genetic (transcription factors) and epigenetic (histone modifications) pathways involved in the differentiation of the MSCs to osteoblasts and adipocytes. As outlined in Fig. 1 and detailed under "Experimental Procedures," marrow-derived cells were isolated from C57bl/6 mice, and selected MSCs were cultured and validated as described previously (34,35). MSCs were then expanded for differentiation, staining, ChIP-seq, and RNA-seq analyses. They were then induced for 7 days in adipogenic medium (MSC-F7) or for 7 (MSC-B7) and 15 (MSC-B15) days with osteogenic medium. At these time points, the cells were collected for staining, ChIP-seq, and RNA-seq analysis. Oil Red O (ORO) staining was performed to examine lipid accumulation in adipocytes, and alizarin red (AR) staining was used to detect calcium deposition in mineralizing osteoblasts. The differentiation states displayed the appropriate cell morphology and staining as in Fig. 1.
Having established the anticipated morphology, we then confirmed appropriate transcriptomic changes as well. As can be seen in Fig. 2, differentiated MSCs were easily distinguished as either osteoblasts (MSC-B7) or adipocytes (MSC-F7) from their gene expression profiles. These patterns identified in MSC, MSC-B7, MSC-B15, and MSC-F7 cells are displayed as heat maps from a hierarchical clustering analysis (using Euclidean distance) of the top 200 up-regulated (top) and down-reg-ulated (bottom) genes contrasting MSC-B7 versus MSC ( Fig.  2A) and MSC-F7 versus MSC (Fig. 2B). A large contingent of the genes up-regulated in MSC-B7 cells was also up-regulated in MSC-B15 cells, as can be observed by the gene clusters. This relationship is better visualized through the Venn diagrams documented in Fig. 2C. Here, genes that were up-regulated Ͼ2-fold within a 95% confidence interval (moderated t test) versus MSCs are shown on the left. A total of 373 genes were up-regulated in MSC-B7 cells (versus MSCs), 658 genes were up-regulated in MSC-B15 (versus MSC), and 1,312 genes were up-regulated in MSC-F7 cells (versus MSCs). Similarly, an analysis of the down-regulated genes revealed that there were many more genes suppressed in each condition: MSC-B7 with 9,152, MSC-B15 with 8,005, and MSC-F7 with 12,711 genes. As expected, both the up-and down-regulated genes of MSC-B7 and MSC-B15 cells displayed more closely related gene expression patterns than those with MSC-F7 cells; however, many genes were unique to both the MSC-B7 and MSC-B15 differentiation states.
We then examined these gene networks through the Gene-MANIA (Gene Multiple Association Network Integration Algorithm) pathway and literature analysis tool (36). As shown in Fig. 2D, the top osteogenic genes shared in common with both the MSC-B7 and MSC-B15 cells are displayed. The genes within the black circles are those explored in GeneMANIA. These were computed to exhibit connections between each other as well as with genes predicted to be within a gene regulatory network (gray circles). We found that these highly regu-lated genes were tightly correlated with not only the algorithm prediction (orange lines) but also genes known to be co-expressed (purple) and/or co-localized (blue). Furthermore, this gene network can be interrogated by gene ontology (GO) analysis through the DAVID (Database for Annotation, Visualization, and Integrated Discovery) functional analysis tool (37,38). The gene network associated with osteogenic differentiation was highly enriched for GO terms related to skeletal formation, ossification, and extracellular structure organization (Fig. 2E). In fact, osteocytic genes shared with MSC-B7 and MSC-B15 cells were further increased in the MSC-B15 cells, suggesting that by day 15, the MSCs were either progressing toward or had become osteocytes, as recently suggested (39). Interestingly, gremlin-1 (Grem1) expression was identified following osteoblast differentiation; Grem1-positive MSCs are hypothesized to be unable to differentiate into adipocytes (40). Similarly, genes up-regulated in MSC-F7 cells versus MSCs were also interrogated using the GeneMANIA pathway analysis (Fig. 2F) as well as DAVID functional analysis tools (Fig. 2G). The interactions within the adipogenic network of genes were tighter than those within MSC-B7 and MSC-B15 cells. These genes were highly enriched in GO categories involved in adipogenesis. Thus, the MSC-F7 transition was the more extreme differentiation process, although we suspect that associated transcriptomic changes in MSCs of fat origin could well be different. Taken together, the transcriptome analyses indicate that MSCs express the necessary gene profiles for osteoblasts, osteocytes, and adipocytes and that larger transcriptomic changes are observed during marrow-derived MSC differentiation into adipocytes than into bone cells.
MSC Differentiation Is Driven by Specific and Unique Transcription Factor Complexes-Transcriptomic changes (gene expression) are a direct consequence of genomic regulation. Transcription factors that respond to cellular stimuli are able to bind to the genome and activate the transcription of lineage-or stimuli-dependent genes. Here, we sought to investigate lineage-dependent transcription factors and to define a role for a secondary effector in hormone-dependent VDR activities. Therefore, MSCs were differentiated into MSC-B7, MSC-B15, or MSC-F7 cells and then treated with either vehicle (ethanol) or 100 nM 1,25(OH) 2 D 3 for 3 h, as indicated. ChIP-seq analysis was performed at each of these differentiation time points for VDR, RXR, RUNX2, C/EBP␤, C/EBP␣, and PPAR␥ in MSC, MSC-B7, and MSC-F7 cells. Table 1 provides direct genomic binding quantitation for these factors as overlapping or unique ChIP-seq transcription factor DNA binding peaks after ChIPseq interrogation (bioinformatically determined genomic locations where protein-DNA binding was enriched ϳ200 -300 bp in length) between each differentiation state. Tabular Venn diagrams depict the total number of peaks found in each data set for each antibody. For example, there were 1,752 total peaks for VDR in the unstimulated (vehicle) condition in MSCs. This number was dramatically reduced to only 245 peaks in MSC-B7 cells. 241 sites were found to overlap in each condition, leaving 1,511 peaks exclusively in MSCs and only four peaks unique to MSC-B7 cells. This pattern was similar for the remaining transcription factors. C/EBP␣ and PPAR␥ were not investigated in MSC-B15 cells.  Recently, we found that VDR/RXR, RUNX2, and C/EBP␤ were highly correlated to sites near genes that were important for the differentiation of preosteoblastic MC3T3-E1 cells into more mature mineralizing osteoblasts and described this as an osteoblast enhancer complex (15,16). In MSCs, we found a similar relationship. As can be seen in Fig. 3A, peak locations for VDR/RXR (9,599 1,25(OH) 2 D 3 -stimulated peaks) overlapped those with RUNX2 (33,266 peaks) and C/EBP␤ (46,232 peaks After differentiation, whereas the cistromes for each of these factors were quantitatively lowered, especially in the case of the VDR, there were 177 genomic locations in which VDR/RXR, RUNX2, and C/EBP␤ were found to overlap. In fact, the 15,000 RUNX2 and C/EBP␤ binding sites that overlapped in MSCs were reduced to ϳ6,000 upon differentiation (MSC-B7).
Similarly, there was a high degree of overlap in the factors responsible for adipogenic differentiation. In Fig. 3B, whereas C/EBP␣, C/EBP␤, and PPAR␥ were found to overlap at 11,542 locations in the undifferentiated MSC, this number was significantly reduced to 6,494 regions in MSC-F7 cells, led by strong decreases in C/EBP␤ binding activity post-differentiation (from 46,232 to 22,986). As Fig. 3B and Table 1 illustrate, although the total number of C/EBP␣ peaks remained constant near 40,000 in both the MSCs and MSC-F7 cells, the locations of these sites changed dramatically, resulting in only 22,181 overlapping peaks with nearly 20,000 peaks unique to both MSCs and MSC-F7 cells. The same observation was made for PPAR␥. However, in both cases, the primary drivers of adipogenesis, C/EBP␣, and PPAR␥ remain highly correlated (41). Finally, because genomic programs remain present in MSCs to drive subsequent multilineage differentiation, the overlap between VDR/RXR, RUNX2, C/EBP␤, and PPAR␥ was examined. These factors assembled collectively at 4,173 sites across . The resulting RNA was reverse-transcribed and analyzed by sequencing. The top 200 genes that were up-and down-regulated were clustered using hierarchical clustering followed by Euclidean distance metric and displayed in a heat map for those osteogenic genes (MSC-B7) (A) and adipogenic genes (MSC-F7) (B) on a log 2 scale as indicated. C, gene expression was compared for each differentiation state versus MSC (day 0) and displayed as -fold change that was greater than both 2-fold and 95% moderated t test confidence interval (CI). Overlapping genes are represented in the Venn diagram (up-regulated genes, left, red; down-regulated genes, right, blue) with total regulated genes listed under the differentiation state. D, gene association was investigated through GeneMANIA interactions, where input genes (black circles) were predicted to interact with each other as well as proposed networked genes (gray circles) through algorithm prediction (yellow lines), co-expression (pink lines), co-localization (blue lines), or shared protein domains (green) from the literature. Higher correlation results in a denser network with the gene connection lines being shorter. E, gene association was also interrogated for osteogenic genes through DAVID functional analysis tools, which led to the enriched GO categories as listed with correlated p values (right). The top GO terms are listed for the up-regulated (left) and down-regulated genes (right). F, genes up-regulated by adipogenic differentiation were interrogated by Gene-MANIA as in D and G and DAVID GO term correlation as in E.

TABLE 1 Venn diagram data with differentiation state (left) with total peaks indicated in parentheses versus the indicated differentiation state (right) for each specific antibody or treatment (VDR and RXR have vehicle and 1,25(OH) 2 D 3 conditions; the remainder are vehicle-treated cells)
Those peaks found only in the left condition are shown under "Only" on the left; those that are found in both the left and right are shown under "Overlapping"; and those found only in the right condition shown under "Only" on the right. the genome, although many additional unique and subcombinatorial sites also exist. The overlapping peak regions for these and perhaps other factors suggest that these sites represent transcription factor "hot spots" for the regulation of this mesenchymal-derived lineage as recently proposed (41,42). A prime example of overlapping ChIP-seq peaks is evident at the Mmp13 and Adipoq genomic regions. In Fig. 3D, ChIP-seq data are presented where individual data tracks are overlaid (track hubs) for MSCs (pink), MSC-B7 (blue), MSC-B15 (yellow), and MSC-F7 (green) cells normalized by sequence tag number and to input sample. ChIP-seq data tracks are presented for VDR (treated for 3 h with 100 nM 1,25(OH) 2 D 3 ), RXR, C/EBP␤, RUNX2, C/EBP␣, and PPAR␥ antibodies. Our recently published work on the Mmp13 gene locus defined a distinct set of enhancers located upstream of the transcriptional start site (TSS) at Ϫ10, Ϫ20, and Ϫ30 kb for VDR, RUNX2, and C/EBP␤ (43). These sites and factor occupancy are identical in the MSCs, although several of these enhancers also contain additional C/EBP␣ and PPAR␥ binding activities as well. For the Mmp13 gene, however, significant increases in  . E, ChIP-seq peaks are correlated to those genes that were up-or down-regulated. ChIP-seq peak locations for VDR/RXR, RUNX2, and C/EBP␤ (top, 4,567 regions) or C/EBP␣, C/EBP␤, and PPAR␥ (bottom, 11,542 regions) were converted into nearby genes using the GREAT genomic association rules (4,452 genes and 7,869 genes, respectively). These genes were then compared with the genes that were found to be up-regulated (red) or down-regulated (blue) by RNA-seq. The resulting gene lists were then compared using DAVID functional analysis tools, which led to the enriched GO categories as listed with correlated p values (right).
factor binding were not apparent for any particular cell differentiation state because the Mmp13 gene is expressed and upregulated in each of these cell types (supplemental Table 1). This was in stark contrast to the regulation of the Adipoq gene and transcription factor binding activity within its locus. The Adipoq gene is highly up-regulated in MSC-F7 cells (115-fold) and controlled in part through the binding of C/EBP␣, PPAR␥, and its promiscuous RXR heterodimer binding partner. These findings are apparent in the ChIP-seq data set documented in Fig. 3D (right). There were, however, increases in PPAR␥ and RXR binding near and immediately upstream of the Adipoq TSS specifically in MSC-F7 cells, as indicated by increases in peak amplitude (enrichment of green color). Despite these increases in MSC-F7 cells versus MSCs (115-fold), there appear to be low expression levels of Adipoq in MSCs and MSC-B7 cells. These residually bound transcription factors may facilitate basal expression. Finally, when ChIP-seq analyses and RNA-seq gene expression data were combined, genes were enriched for osteoblastor adipocyte-specific pathways. In Fig. 3E, ChIP-seq peaks are linked to nearby genes through the Genomic Regions Enrichment of Annotations Tool (GREAT) (44). We used the GREAT algorithm and the "basal plus extension" standard rule wherein each gene is assigned to a basal regulatory domain of 5 kb upstream and 1 kb downstream of the TSS and the gene regulatory domain is then extended in both directions to the nearest gene's basal domain by no more than 1 Mb in a single direction. 4,567 genomic locations contained overlapping peaks of VDR/ RXR, RUNX2, and C/EBP␤ and were localized to 4,452 genes. Thus, many genes contain multiple binding sites within their genomic loci as determined by GREAT. The intersection of these peak-associated genes with those that were actually regulated yields only 172 up-regulated (or 46% of the total) and 1,018 (or 11%) down-regulated genes through differentiation (independent of any vitamin D activity). These gene groups were then evaluated by the DAVID functional analysis tool and displayed as enriched GO group terms (38,45). Genes up-regulated in MSC-B7 cells were highly correlated to skeletal formation, with groups such as skeletal system development, ossification, and cartilage development being enriched and containing specific p value scores Ͻ Ͻ0.05. Conversely, peaks associated with the down-regulated genes were involved in broad generic categories with lower enrichment scores. 11,542 peaks located at the intersection of C/EBP␣, C/EBP␤, and PPAR␥ were found to be associated with 7,869 genes. These peak-associated genes were coordinated with 46% of the upregulated genes and only 25% of the down-regulated genes. Upregulated genes were associated with adipogenic processes, such as energy creation, and down-regulated genes were found to be involved in cell adhesion and migration. Up-regulated genes were better correlated with the gene regions that contained peaks of transcription factor binding. Although this does not define a role for transcription factor peaks near genes that are not regulated, many of these genes are still expressed within the different lineages (linear RNA-seq RPKM value Ͼ1). Therefore, basal (not statistically inducible) expression may still require transcription factor proximity.
Contribution to Gene Regulation by the Vitamin D Hormone-VDR is widely expressed in many tissues and can be up-regulated at both the transcriptional and post-translational levels upon ligand binding, the latter probably through a stabilization mechanism imposed by 1,25(OH) 2 D 3 (46). The VDR is a ligand-inducible steroid hormone receptor and may have a mechanism of DNA binding and gene activation different from that of factors such as RUNX2 or C/EBP␤ (15,16), which are generally prebound to DNA and activated through post-translational modification. VDR appears to be restrictive for differentiation (15) yet exhibits pervasive binding activity throughout the genome. Recently, we found that the levels of VDR decreased in differentiated osteoblasts, although lowered concentrations of VDR were still sufficient to activate gene transcription (15). We therefore assessed these behaviors of the VDR during the course of MSC differentiation, as can be seen in Fig. 4. A similar number of genes were statistically up-regulated in MSCs (173 genes) and in MSC-B7 (311 genes), MSC-B15 (218 genes), and MSC-F7 (293 genes) cells after a 24-h treatment with 100 nM 1,25(OH) 2 D 3 . However, the number of genes down-regulated by 1,25(OH) 2 D 3 after differentiation decreased dramatically, from 337 to 95, 56, and 99, respectively. The overlaps of these regulated genes appear in Venn diagram form in Fig. 4A. Surprisingly, although only 11 genes were up-regulated (Ͼ2-fold, Ͼ95% confidence interval) in all three differentiation states (Cyp24a1, Prrg4, Sema3e, Fgl2, Inmt, Gpdp5, Neu3, Enpp3, Mdfi, Klhl33, and Grk5), only a single gene was statistically down-regulated (Ibsp). The complete gene lists can be found in supplemental Table 1.
Similar to the observations made in MC3T3-E1 cells (15), the genomic binding activity of the VDR is decreased dramatically after differentiation. Protein levels of VDR, however, were suppressed to a lesser degree in MSC-B7 and MSC-B15 cells (data not shown) than in MC3T3-E1 cells. In Fig. 4B, the ChIP-seq peak data for VDR binding in both vehicle (Veh)-and 1,25(OH) 2 D 3 (1,25D 3 )-treated cells is summarized. VDR binding increased in MSCs from 1,752 peaks to 10,557 peaks, which is consistent with the ligand-activated nature of VDR DNA recruitment. Although the numbers of VDR binding sites were not identical, much of the MSC data overlapped that seen in MC3T3-E1 cells (vehicle, 947; 1,25(OH) 2 D 3 , 7,007) at Ͼ90% of all genomic locations (data not shown). Despite the decrease in genomic binding of VDR in each of the differentiation states, however, this binding remained generally ligand-dependent, although a modest but significant VDR binding activity at selected sites on the genome in the absence of 1,25(OH) 2 D 3 is present in multiple cell types and was evident in MSCs as well (15,47,48). 1,408 of the 1,752 VDR peaks seen in the MSCs in the absence of 1,25(OH) 2 D 3 were shared with the 10,557 peaks detected after 1,25(OH) 2 D 3 stimulation, leaving only 344 peak regions occupied by VDR across the genome in the absence of ligand. Finally, the VDR DNA binding partner RXR was found at greater than 90% of the VDR peak locations, and the remaining 10% exhibited a level of RXR enrichment that was evident but did not reach statistical significance (data not shown, available in raw ChIP-seq reads). Importantly, the most common de novo motif found by HOMER sequence overrepresentation in each differentiation state belonged to a canonical DR3 type   VDRE (Fig. 4D). However, a further position weight matrixbased search of these ChIP-seq peaks suggested that 92% of the VDR peaks contained a DR3 type VDRE, as found by Genomatix (VDR_RXR.01 and VDR_RXR.04, core similarity 0.75, matrix similarity optimum Ϫ0.10) (49). Therefore, there is limited evidence for an abundance of alternative DNA binding modes (homodimer, tethering, non-canonical VDRE) for the VDR aside from the canonical VDR/RXR heterodimer. Although the VDR is present and ligand-activated in all cell types, the genomic impact of the VDR appeared to be greatest in undifferentiated MSCs. GREAT association of VDR genomic ChIP-seq peaks with neighboring genes resulted in 9,599 peaks being associated with 7,167 genes. This highlights the likelihood that many genes have multiple peaks of VDR/RXR binding, as was demonstrated for PPAR␥, RUNX2, C/EBP␤, and C/EBP␣ binding in Fig. 3E. There is a higher correlation with ChIP-seq peaks (76%) and genes that are up-regulated by 1,25(OH) 2 D 3 . The most enriched GO term categories for the up-and down-regulated genes are shown in the right panel of Fig. 4E. According to these enriched categories, 1,25(OH) 2 D 3 -mediated gene expression does not favor osteoblast-specific categories. Finally, ChIP-seq data tracks are presented for three genes that were up-regulated by gene expression in each of the differentiation states: MyoD family inhibitor (Mdfi) (Fig. 4F), proline-rich Gla protein 4 (Prrg4) (Fig. 4G), and glycerophosphodiester phosphodiesterase domain-containing 5 (Gdpd5) (Fig. 4H). The overlaid ChIPseq data tracks for VDR in the vehicle (Veh, yellow) and 1,25(OH) 2 D 3 (1,25D 3 , blue) conditions (where overlapping tracks appear green) are displayed for each of the differentiation states of the MSCs as well as for RXR in MSCs only for comparison. There was a strong enrichment of VDR following treatment with 1,25(OH) 2 D 3 , as indicated in the ChIP-seq data hub tracks. Whereas Mdfi and Prrg4 loci display single peaks of VDR/RXR binding activity at Ϫ8 and Ϫ15 kb, respectively, it appears as though the Gdpd5 locus contains several peaks of VDR/RXR binding activity that may comprise enhancers for the control of expression of this gene. This latter model was found to be important in the regulation of many vitamin D-responsive genes, such as Mmp13, Vdr, c-Fos, Tnfsf11, Trpv6, and others (15, 43, 47, 50 -53).
Epigenetic Transitions Predict Crucial Regulatory Loci during Differentiation-Transcription factors bind to the genome in a punctate pattern and are vital mediators of external cellular stimuli. These events can be preceded by, or are an effect of, the chromatin landscape in which they reside; therefore, histone modifications can be accurate predictors of the genomic transitions (26). We investigated five different chromatin modification marks that specify important genomic locations, such as transcriptional start sites (H3K4me3) or enhancer signatures (H3K4me1), or genomic events such as transcriptional elongation (H3K36me3) or gene activation (H4K5ac and H3K9ac) (26). In Fig. 5A, ChIP-seq peaks are tabulated for these chromatin modification marks in the absence or presence of 1,25(OH) 2 D 3 for each differentiation state. From a cis-regulatory element annotation system (CEAS) analysis in Fig. 5B, we found that these histone marks followed the expected genomic distribution profiles when all gene bodies throughout the genome are averaged together to display the ChIP-seq tag densities as reported previously (54). During differentiation, peaks shifted in amplitude (greater or lesser peak reads) but not necessarily in location, similar to recent observations (55). This was demonstrated in Fig. 5C, where the peak overlap of the H3K9ac mark in MSCs versus MSC-B7 cells, as one example, represents 99 and 92% of all peak locations, respectively. This was in contrast to the H3K9ac mark, where 76 and 87% of peaks were found to overlap following differentiation of MSCs into MSC-F7 cells, leaving 7,543 sites unique to MSCs and 2,760 sites unique to MSC-F7 cells. These relationships were observed in the remaining histone modification marks as well (ChIP-seq raw reads). Therefore, MSC-B7 cells exhibited an epigenetic profile that is similar to that of undifferentiated MSCs, whereas the profile for MSC-F7 cells is dramatically different.
Accompanying changes in the transcriptome, there were also differences in the amplitude of the epigenetic marks when comparing MSCs with either MSC-F7 cells or MSC-B7 cells, as seen in Fig. 5D. Here, TSSs were averaged across the genome (from ϩ3 to Ϫ3 kb) using the CEAS analysis tool and compared using ChIP-seq tag densities, shown as "average profile." Averaged H3K9ac data that span the TSSs of all genes appear as a black line and display a typical H3K9ac profile that is enriched before and immediately after the TSS (valley of activity directly over the TSS). However, an examination of the TSS region from genes that are regulated during differentiation into MSC-B7 cells (up-regulated, blue; down-regulated, yellow) and MSC-F7 cells (up-regulated, green; down-regulated, dark green) revealed specific enrichments of the H3K9ac mark at genes regulated in MSC-B7 compared with MSC-F7 cells (left). For genes regulated in MSC-F7 cells (right), H3K9ac enrichment was much greater than that seen in MSC-B7 cells. These genome-wide observations are best illustrated at the gene level at two loci that surround the Adipoq (adiponectin) (Fig. 5E) and Lipe (hormone-sensitive lipase) genes (Fig. 5F), which are highly and specifically up-regulated in MSC-F7 cells. ChIP-seq data overlaid track hubs for each differentiation state (MSC (pink), MSC-B7 (blue), MSC-B15 (yellow), and MSC-F7 (green)) for H3K9ac, H4K5ac, H3K4me1, H3K4me3, and H3K36me3 revealed an increase in each histone mark specifically in MSC-F7 cells (enrichment of green). Although enrichment at the Adipoq locus was highly punctate for H4K5ac and H3K4me1 directly upstream of the TSS, overlapping PPAR␥ and C/EBP␣ binding sites, these histone marks were more broadly distributed across the Lipe gene locus.
Histone modifications are also predictive of activated loci after treatment with 1,25(OH) 2 D 3 . In Fig. 5G, CEAS-averaged TSS analysis was applied to genes that were up-or down-regulated by 1,25(OH) 2 D 3 treatment. Genes that were activated or up-regulated (red) by 1,25(OH) 2 D 3 exhibited enrichments for both H3K9ac (left) and H3K4me3 (right) that were specific for 1,25(OH) 2 D 3 (solid red line). Genes that were down-regulated (gold lines; 1,25(OH) 2 D 3 (solid) and vehicle (dashed)) showed no enrichment for these histone modifications, which are typically associated with active chromatin. These data are highlighted through the ChIP-seq overlaid data tracks visualized in Fig. 5H

Gene expr -1,25(OH) 2 D 3 vs Veh Gene expr -1,25(OH) 2 D 3 vs Veh H3K9Ac
H3K4me3 played near the Mdfi gene, which was recently profiled for its transcription factor binding (15). Mdfi is shown to contain potential enhancers near the TSS and upstream at Ϫ8 kb and is up-regulated by 1,25(OH) 2 D 3 in virtually all of the differentiation states examined (MSC, 11-fold; MSC-B7, 13-fold; MSB-B15, 4.5-fold; MSC-F7, 12-fold; supplemental Table 1). The ChIP-seq data hub tracks, which contrast cells treated with either vehicle or 1,25(OH) 2 D 3 , reveal 1,25(OH) 2 D 3 -induced enrichment of H3K9ac (blue) at the Ϫ8 kb region in all differentiation states. H3K4me3 was selectively enriched only immediately downstream of the Mdfi TSS. Therefore, undifferentiated MSCs exhibit an epigenetic signature that resembles that seen in MSC-B7 cells, and histone modifications appear to be predictive of activated gene loci following both differentiation and hormone-induced treatments.
Epigenetic Patterns in MSCs Predict Lineage Preference for Osteoblastic Differentiation-There is growing evidence for an epigenetic preference that favors the differentiation of marrowderived MSCs toward the osteoblast lineage in genes that are up-regulated by osteogenic differentiation. This stands in stark contrast to genes in MSCs that are selectively up-regulated following adipogenic differentiation. An excellent example of this contrast can be seen at the gene locus that encodes osteopontin (OPN, Spp1) in Fig. 6A as compared with the gene locus for Cebpa in Fig. 6B. ChIP-seq data hub tracks are presented for RXR, C/EBP␤, RUNX2, C/EBP␣, PPAR␥, H3K9ac, H4K5ac, H3K4me1, H3K4me3, and H3K36me3 (top) in each differentiation state (MSC (pink), MSC-B7 (blue), MSC-B15 (yellow), and MSC-F7 (green)) as well as VDR and RXR (bottom) in the absence or presence of 1,25(OH) 2 D 3 . Spp1 is up-regulated following osteogenic differentiation (2-4-fold) and by 1,25 (OH) 2 D 3 (4-fold) yet down-regulated during adipogenic differ-entiation (2-fold). The Spp1 genomic locus retains several putative enhancers located upstream of the TSS that are occupied by VDR/RXR, RUNX2, and C/EBP␤, as documented previously (15). The Spp1 gene is expressed in a wide variety of cell types, however, and may also manifest components responsible for cell type-selective expression. Thus, the more proximal enhancers contained ligand-activated VDR/RXR recruitment and histone modifications that remain unchanged relative to undifferentiated MSCs. However, enhancers at the Ϫ27 and Ϫ43 kb regions are also seen to bind additional transcription factors that include C/EBP␤ and RUNX2 as well as PPAR␥ and C/EBP␣. At these sites, the histone profiles for H3K9ac, H4K5ac, and H3K4me1 are unexpectedly more robust in MSC-F7 cells as compared with MSCs (green increased). The separation of these epigenetic marks to more proximal or more distal regions may provide distinct mechanisms of upregulation (proximal) and down-regulation (distal) at the Spp1 gene in different cell types. The locus near the Cebpa gene, on the other hand, which is up-regulated in MSC-F7 cells ϳ3-fold compared with its down-regulation in MSC-B7 cells (4-fold down-regulation) and MSC-B15 cells (11-fold down-regulation), is dominated by an increase in histone marks in MSC-F7 cells (green) and an increase in PPAR␥ at sites downstream of Cebpa. There is also a ligand-dependent binding of VDR/RXR downstream of Cebpa, despite the fact that Cebpa does not appear to be statistically up-or downregulated by 1,25(OH) 2 D 3 in any of these differentiated cell types (supplemental Table 1). We are therefore able to use chromatin signature marks similar to those found near Spp1 and Cebpa to accurately predict important genomic loci that might be essential for differentiation. Transcriptomic and Epigenetic Plasticity in Trans-differentiated Cells-Following the genome-wide observation that MSC and MSC-B7 epigenetic profiles are similar, whereas MSC to MSC-F7 profiles exhibit more dramatic differences, we hypothesized that adipogenic differentiation might represent a more committed terminal state of differentiation than that seen in the osteoblast, a subset of which are known to differentiate further into bone lining cells or osteocytes. To investigate this hypothesis, we conducted a series of trans-differentiation experiments wherein cells were grown in osteogenic (OM), adipogenic (AM), or standard growth medium lacking any added compounds (Basal) (Fig. 7) for 3 days and then cultured for an additional 7 days in medium that induced the alternative state of differentiation. We then stained the cells with either AR or ORO, isolated RNA for qPCR gene expression analysis, and subjected the cells to ChIP-qPCR analysis at days 0, 3, and 10, as summarized in Fig. 7A. As documented in Fig. 7B, we found that after 3 days in osteogenic medium, AR staining revealed, as expected, an extensive calcium deposition that resulted in dense bone nodules after 10 days in culture. When transferred to adipogenic medium, however, the cells began to form lipid droplets, as assessed by ORO staining as seen in Fig. 7B (right), indicating a change in morphology toward adipogenic cells. ORO staining was positive only in cells that had been exposed to the adipogenic medium and was highly efficient. In contrast, cells exposed to adipogenic medium (Fig. 7B) displayed lipid droplet formation and positive staining for ORO both at 3 days and at 10 days. However, adipogenic cells that were transferred back to basal (control) medium displayed a reduction in lipid droplet formation, whereas this stain was absent in the adipogenic cells that were transferred into osteogenic medium for 7 days. Indeed, the latter cells began to form dense bone mineral nodules and stained positive for AR.
To confirm potential transcriptomic signatures, RT-qPCR was performed on samples acquired at 0, 3, and 10 days for genes relevant to osteogenic biology (Spp1, Mmp13, and Runx2) as well as those important for adipogenic identity (Cebpb, Pparg, Adipoq, Cebpa, and Lipe), displayed in Fig. 7C. Spp1 expression was increased 3-fold in control medium (tan bounding box) after 3 days (first tan bar), unchanged after an

Spp1 (OPN)
+7D +7D . Cellular trans-differentiation was accompanied by morphological and transcriptomic changes. A, experimental design and overview. Marrowderived cells were isolated, and MSCs were preserved through culturing. MSCs were then differentiated in osteogenic (OM) or adipogenic medium (AM) for 3 days. At this time, the cells were then switched to the opposing differentiation mixture for an additional 7 days (10 days total). Gene expression, ChIP assay, and staining were performed at days 0, 3, and 10. B, differentiated cells were stained with AR (left) and with ORO (right). All images were at 4ϫ magnification. C, RT-qPCR was performed on cells that were differentiated throughout the time course focused on Spp1, Mmp13, Runx2, Cebpb, Pparg, Adipoq, Cebpa, and Lipe genes. MSCs without differentiation (D0, blue) are compared with cells cultured for 3 days in basal medium (Basal 3D, tan bounding box and first tan bar), osteogenic medium (OM 3D, gray bounding box and first gray bar), or adipogenic medium (AM 3D, green bounding box and first green bar). Cells were then differentiated for an additional 7 days as indicated with basal medium (B, tan), osteogenic medium (O, gray), or adipogenic medium (A, green). Data were completed in triplicate analysis Ϯ S.E. Data are displayed as relative quantitation (RQ) normalized to ␤-actin (⌬⌬C T ). a, p Ͻ 0.05, 3 days versus 0 days by paired t test. *, p Ͻ 0.05, 7-day treatment versus 3 days within medium group by paired t test. additional 7 days (second tan bar) in control medium, and further increased in osteogenic medium (black bar). The transition into adipogenic medium decreased Spp1 expression 4-fold (first tan bar versus green bar). In cells initially exposed to osteogenic medium (gray bounding box), Spp1 was elevated in response to a change to basal medium and elevated even further when subjected to osteogenic medium for a further 10 days. As with basal medium, the further exposure of osteogenic cells to adipogenic medium resulted in a strong repression of Spp1. Finally, Spp1 expression in adipogenic cells was only slightly elevated when the cells were further exposed to either basal or osteogenic medium for additional periods of time. Continual adipogenic medium exposure resulted in sustained suppression of Spp1 at 3 and 10 days. These relationships were also observed for the Mmp13 and Runx2 genes. Expression profiles for Cebpb, Pparg, Adipoq, Cebpa, and Lipe were dependent upon the presence of adipogenic medium and signaling and were increased regardless of when the adipogenic medium was applied.
In a final study, we utilized ChIP-qPCR analysis to evaluate the effects of trans-differentiation on key histone modifications. The experimental design and display were identical to those seen in Fig. 7C. Cells were subjected to a ChIP assay at days 0, 3, and 10 and immunoprecipitated with antibodies to H3K9ac, H4K5ac, H3K4me1, H3K4me3, and H3K36me3 as well as with a nonspecific IgG control. As seen in Fig. 8A, the ChIP-qPCR primers used were located near the Spp1 TSS (Ϫ200/ϩ100 bp), and the data displayed represent ChIP analysis for H3K9ac (top) and IgG control (bottom) at this site. The ChIP-qPCR results seen in Fig. 8A mirror the results of Spp1 gene expression seen in Fig. 7C. H3K9ac levels remained unchanged in cells that were treated with either basal (tan bounding box) or osteogenic (gray bounding box) medium for 3-10 days. Cells cultured in adipogenic medium (green bounding box), however, displayed a 4 -5-fold decrease in histone acetylation, indicating a loss of transcriptional activity. The nonspecific IgG immunoprecipitation data demonstrate that these activities were specific because they remained unchanged during treatment conditions and were significantly lower than all immunoprecipitation samples. The IgG for the remaining regions was as low or lower than the data displayed for Spp1; therefore, this was omitted from the figure for clarity. In Fig. 8B, ChIP-qPCR analysis of marks for both H3K9ac and H3K4me3

Spp1 TSS -H3K9ac
Basal 3D  was also performed using primers located near a region ϩ500 bp downstream of the Lipe gene TSS. As seen in Fig. 8B, the data derived from this analysis also mirrored that seen for the expression of the Lipe gene in Fig. 7C. Histone marks were increased only in cells cultured in adipogenic medium, and these marks were increased even after mineral matrix began to form (OM 3D). Histone marks at the loci for the Cebpa (ϩ2 kb) and Adipoq (TSS ϩ0.5 kb and ϩ11 kb) genes were also increased when adipogenic medium was applied regardless of the differentiation state. Taken together, the flexibility or plasticity of the epigenetic profiles enables MSCs to respond to alternative stimuli rather than to commit to a terminal differentiation process. Thus, MSCs are able to display cellular characteristics appropriate to mature mineralizing osteoblasts even after displaying adipogenic features that include the formation of ORO-positive lipid droplets.

Discussion
Recent investigations into the differentiation of the osteoblast have revealed the presence of an intricate relationship between the transcription factors RUNX2, C/EBP␤, and VDR/ RXR that is formed around genes important for osteoblast function and differentiation (15,16). These studies were initially performed in the MC3T3-E1 cell model system in which the transition of preosteoblasts into mature mineralizing osteoblast cells was investigated. Similarly, there have been several studies utilizing the preadipocyte cell line 3T3-L1 to explore adipocyte differentiation with PPAR␥, C/EBP␣, and numerous additional transcription factors, including the VDR, as targets (22,56,57). The current model system of marrow-derived MSCs presented herein displays distinct advantages over both the MC3T3-E1 and the 3T3-L1 cell models. First, the multipotential nature of the MSCs allows an investigation of cellular transition into several different cell types using the same starting material. Second, marrow-derived MSCs can be isolated from transgenic or knock-out animal models, permitting investigations of the contributions of the genetic manipulations to the growth and differentiation of the cells (34). Finally, due to their accessibility from patients, the differentiation and reprogramming of MSCs ex vivo have become a focus in clinical therapeutics and drug discovery. For these reasons, we sought a thorough investigation of the mechanisms of gene activation, expression, genome binding, and epigenetic transitions in the differentiating MSCs and to provide a dense archive of ChIPseq and RNA-seq that can be further deciphered by the research community because data on bone tissues are absent or limited from ENCODE data sets.
As we observed previously in the MC3T3-E1 cells (16) and as evidenced from the 3T3-L1 literature (22,56,57), transcription factors essential for osteogenic (RUNX2 and C/EBP␤) and/or adipogenic (PPAR␥ and C/EBP␣) differentiation appear poised or are already present at many locations through the genome in undifferentiated MSCs regardless of their origin. During the process of differentiation, these factors either increase/decrease in abundance at specific enhancers or, with less frequency, bind at de novo sites. The peak overlap table (Table 1) demonstrates that although there are unique peaks of activity for each differentiation state, the majority of transcription fac-tors are bound to the same enhancer locations, and if changes do occur, they are largely of amplitude. Examples of this behavior were apparent at the Adipoq gene locus seen in Fig. 3D. At this gene, however, PPAR␥ and RXR binding both increased from an existing residual level at specific sites (adjacent to the TSS and Ϫ3 kb) and occupied additional de novo sites (Ϫ1 kb) following MSC-F7 differentiation (green). We also found that many additional factors bound to the Ϫ3 kb enhancer upstream of Adipoq. Binding sites for factors such as VDR, RXR, C/EBP␤, and C/EBP␣ did not change in response to MSC differentiation, and all cell differentiation states were equally overlapping. Similar observations were made at the Mmp13 (Fig. 3D), Spp1 (Fig.  6A), and Cebpa (Fig. 6B) loci. Further studies using ChIP-re-ChIP, where ChIP is performed for one factor, followed by a second immunoprecipitation for a co-bound factor, would help reveal which of these factors are present at the same time and physically interacting on the same piece of DNA to further our understanding of the mechanism (58). However, it is important to remember that these factors may not be directly related to activity of genes in the same genomic vicinity. Interactions of chromatin modifiers (readers and writers for the histone marks) and other coactivators also operate during the activation. Perhaps other transcription factors and/or co-repressor presence are necessary for repression. Nonetheless, the presence and location of these essential factors provide the first level of mechanistic understanding, and their activities highlight the complexity of regulation and the potential coordination of several transcription factors that are frequently observed at thousands of locations throughout the genome.
Ligand activation of the VDR, a secondary activator and signal-dependent factor, and genomic binding with its RXR partner create a different gene activation dynamic in differentiated cells compared with parental MSCs. Although the genomic footprint of ligand-activated VDR is dramatically reduced (from 10,557 to 823 (MSC-B7), 1,365 (MSC-B15), and 616 (MSC-F7)), 1,25(OH) 2 D 3 -activated VDR is able to regulate similar numbers of genes in each cell type. These findings are consistent with data generated in the MC3T3-E1 cell line (15,16) when these preosteoblast precursor cells were differentiated into mature osteoblastic mineralizing cells in osteogenic medium. The gene networks activated by 1,25(OH) 2 D 3 in each differentiated cell type derived from the MSCs were distinct, with only 11 genes being up-regulated and one gene (Ibsp) down-regulated in common in each differentiation state. Finally, examination of the most common or overrepresented sequences in these peaks revealed that a typical VDRE was found to be the most abundant sequence located in 20 -35% of the total peaks, which lends credence to the fidelity of this data set. These analyses, however, should not be misinterpreted to suggest that a VDRE was not found in the remaining 80 -65% of the ChIP-seq peaks. A position weight matrix-based search found that Ͼ90% of all ChIP-seq peak sequences contained a DR3 type VDRE near the peak maxima or pairs of VDREs flanking the apparent peak center. These data are consistent with our previous conclusions regarding VDR cistromes in human LS180 colonic cells, mouse IDG-SW3 osteocytic cells, and mouse MC3T3-E1 osteoblastic cells (15,47,59). Interestingly, these regions are also enriched in sequences for AP1/Jun, RUNX2, and C/EBP␤ motifs in the osteogenic differentiated cells (15). However, in the MSC-F7, sequence motifs for the transcription factors NKX2-5 (NK2 homeobox 5) and IRF5/6 (interferon regulatory factor 5 and 6) motifs were present. The NKX2-5 protein is involved in heart development (60). IRF5/6 transcription factors are involved in inflammatory pathway signal transduction (61). Despite the overrepresentation of these motifs in our data sets, each pathway will require extensive investigation to determine their contribution to VDR hormone signaling in adipocytes.
From a broad level of expression, epigenetic changes were a predictor of the important gene loci for differentiation. These histone marks result from the recruitment and activity of histone-modifying enzymes, such as histone acetyl-and methyltransferases, that are the readers and writers of the genome (62,63). Detailed CEAS analysis of ChIP-seq enrichment near gene TSSs correlated well with genes that were regulated. The epigenetic profiles for MSC-B7 (7 days in osteogenic medium) were nearly indistinguishable from those seen in undifferentiated MSCs at genes that were up-regulated following osteoblast differentiation, including Mmp13 and Spp1, and near the Dmp1 (115-fold increase in MSC-B7 versus MSC) locus (ChIP-seq raw data). This was in contrast to observations made following adipogenic differentiation, where histone marks were specifically enriched around adipogenic genes, such as Adipoq, Lipe, Cidec, Cebpa, Fabp4, Klf15, Plin1, Pparg, Setd8, and many others (supplemental Table 1 and ChIP-seq raw data). In these studies, we focused on the process of activation in response to differentiation. Each histone mark provided insight into different aspects of chromatin regulation, such that ChIP-seq analysis of candidate transcription factors may not be necessary initially to understand the importance of the gene locus for differentiation. The H3K9ac and H4K5ac marks associated with active enhancers were duplicative in nature (activation); however, each mark revealed unique profiles. At the Adipoq gene locus, the H3K9ac mark exhibited its highest enrichment near the TSS coinciding with H3K4me3, whereas H4K5ac, while enriched at the TSS, was also increased at the upstream enhancers that also coincided with the H3K4me1 mark. Coupled with our factor occupancy data, it is clear that these regions are important for PPAR␥ and C/EBP␣ binding. Nevertheless, the importance of this region as an enhancer is implied even in the absence of these types of data. The dynamics of acetylation profiles (e.g. the punctate patterns at the Adipoq and Cebpa genes versus the broad histone profiles at the Lipe gene) facilitate the categorization of gene activities across the genome and create a hierarchy of genes involved in differentiation. It is possible that a diverse set of multiple enhancers and broad active chromatin may indicate a wider portfolio of responses to cellular stimuli and differentiation. The reverse of these activation marks are those that are involved in repression (H3K27me3) and silencing (H3K9me3). Although we did not examine repressive marks in these studies, their involvement would further illuminate the patterns necessary for the multipotential nature of the MSCs.
We believe that the activation dynamics identified in this study enhance an already attractive cell model system in which to study the epigenetic plasticity of MSC regulation. Whereas there is strong evidence for trans-differentiation of cells in vivo (30 -33), we also sought to define this process in response to the adipogenic or osteoblastic signaling growth mixtures at a mechanistic level. The environment, in our case the differentiation mixture, determined the phenotypic outcome, epigenetic restructuring, and therefore function analogous to macrophages that are deposited in different tissues and appear to acquire unique genetic and epigenetic features from their environment that determine their specific functions (64). This concept could be applicable to MSCs relative to their tissue origins. The epigenetically preprogrammed pathway for marrow-derived MSCs appears to be osteogenic; however, the potential for adipocyte differentiation remains intact, despite expression of the reported adipocyte defying Grem1 in the MSCs and MSC-B7 cells (40). This differentiation potential is highlighted in the pathology associated with skeletal fragility, including increased marrow adiposity due to less physical exertion, aging, and estrogen deficiency (65)(66)(67). The marrow microenvironment and resulting MSCs may react and respond to the increased adiposity, causing an increase in the propensity for MSCs to become adipose or even for existing osteogenic cells to be trans-differentiated into fat. This in turn results in a cascade of signaling pathways being activated and changes in the expression and elaboration of Wnts/␤-catenin, bone morphogenic proteins, and other growth factors that are capable of driving cells toward a clinically undesired outcome (68 -70).
Our original hypothesis that the adipogenic cells might be terminally and irreversibly committed appeared to be inaccurate because we were successful in transitioning adipocytes back into mineral matrix-depositing osteoblast-like cells. The cell morphology and staining patterns in these trans-differentiated cells were consistent with individual osteoblast versus adipocyte phenotypes. These transitions did not occur randomly, however, because differentiation was uniform across the plates, and the gene expression profiles were unequivocal for the appropriate cell type. For example, there were no lingering remnants of osteoblast expression in adipocyte trans-differentiated cells. The most interesting results found during the trans-differentiation process were in relationship to the profiles of histone marks. Despite examining these transitions with ChIP-qPCR rather than ChIP-seq analysis, we still found that epigenetic profiles mirrored the gene expression results at many locations (guided by the ChIP-seq for primer design). Because there was high expression of Adipoq in cells with an adipogenic medium mixture, the histone modifications decreased after the adipogenic signals were removed. In fact, this was true for all genes up-regulated by adipogenesis (Pparg, Cebpa, Lipe, and others (ChIP-seq raw data)), thereby demonstrating that epigenetic plasticity is required for these transitions. Importantly, these patterns of the transitioning MSC may contain unique epigenetic signatures that are dependent upon cellular origin (adipose, osteogenic, or chondrogenic), disease status (healthy or tumorigenic), or even their final site of localization (bone matrix, fat pad, or brain).
The identification of histone marks and accompanying signals, whether in vitro or in vivo, is vital to understanding gene regulation in a broader context. Genome-wide association studies undertake the difficult task of identifying regions of the genome, and more precisely SNPs, that have a higher correla-tion with certain disease risk groups. A majority of the SNPs that have been identified occur in the vastness of the non-coding portions of the genome (71). These SNPs can act as genomic beacons to highlight crucial enhancers, regulators, or structure for gene regulation and disease progression. Many enhancers may be connected through intergenic (or interchromosomal) loops, whereby SNPs may have far reaching, difficult-to-predict effects on gene expression. Our recent data utilizing CRISPR (clustered regularly interspaced short palindromic repeats) in the Mmp13 gene locus showed that disruption of the Ϫ30 kb RUNX2-bound enhancer was responsible for the altered binding activity of several additional transcription factors at the remaining distal enhancers and ultimately inactivation of transcription of the gene itself (43,72). The presence of SNPs in regions that contain histone marks may be important in determining their specific functions and may help to determine mechanisms and perhaps guide therapeutic intervention.
In conclusion, the multipotential nature of MSCs creates a desirable model in which to study the formation of many different tissue types. We found that transcription factors required for each differentiation state were present in the undifferentiated MSCs and were poised for regulation. These genetic responses were the cause of, or reaction to, the epigenetic changes that drive the cellular transitions. The epigenetic profiles of the marrow-derived MSC closely resemble the osteogenic differentiated cells and indicate a default bone fate. The plasticity of these epigenetic marks allows the MSC to maintain its multipotential role and provide precursors for the various differentiated cell types as needed.