Genomic occupancy of Runx2 with global expression profiling identifies a novel dimension to control of osteoblastogenesis

1The development of bone tissue, new bone formation in the adult, mineral homeostasis and maintenance of bone mass are mediated by cells of the osteoblast lineage. These bone forming cells progress through a highly regulated differentiation program, with each subpopulation of cells acquiring stage-specific phenotypes that are characterized by distinct profiles of expressed genes. Osteoblast commitment and differentiation are dependent on the appropriate expression of Runx2 (Runt-related transcription factor 2), the master regulator of bone formation. Ablation of Runx2 in mice results in the absence of a mineralized skeleton, and disruption of Runx2 function causes bone defects in the human disorder cleidocranial dysplasia. Runx2 controls a complex gene-regulatory network during osteoblastogenesis. It upregulates a variety of osteoblast lineage-specific genes, including Osx (osterix), Ocn (osteocalcin), and Bsp (bone sialoprotein), and represses the expression of non-osteoblast genes such as PPAR-γ (peroxisome proliferator-activated receptor gamma) and MyoD (myogenic differentiation), which are required for adipogenic and myogenic commitment, respectively. Studies have demonstrated that Runx2 controls gene expression by interacting with multiple classes of co-regulatory factors. In addition to the traditional transcriptional mechanisms, several epigenetic mechanisms have been identified for Runx2-mediated gene expression. For example, Runx2 is retained on mitotic chromosomes as an epigenetic bookmarking factor to maintain cellular identity after cell division. A significant body of evidence substantiates the contribution of Runx2 towards the epigenetic regulation of gene expression during osteoblast differentiation, through interactions with histone deacetylases, histone acetyltransferases, and SWI/SNF complex components.

2Runx2 has been shown to activate or repress gene expression through binding to cis-regulatory DNA elements, the Runx motif (TGTGGT) located in or near gene promoter regions. Genome-wide binding profiles of the transcriptional factor CTCF and lineage-specific transcription factors, such as PPAR-γ, MyoD, and GATA3 (GATA binding protein 3), via ChIP-Seq (chromatin immunoprecipitation followed by high-throughput sequencing) have helped to delineate cis-regulatory networks that are critical for cell lineage control in adipocytes, myocytes and T cells, respectively. These studies have highlighted the regulatory importance of long-range interactions and binding of phenotypic transcription factors to non-promoter genomic elements. These findings have greatly expanded existing paradigms of transcriptional regulation. In contrast, the entire scope of Runx2 binding elements remains largely unknown, hindering a comprehensive understanding of the cis-regulatory network through which Runx2 regulates the transcription program for bone formation.

3Hence, we have characterized the genome-wide occupancy of Runx2 by ChIP-Seq in MC3T3-E1 preosteoblasts, a well-studied in vitro model for osteoblastogenesis. Runx2 occupancy was determined at three hallmark stages of osteoblast differentiation: proliferation, matrix deposition, and mineralization. Analyses of these data demonstrated that Runx2 occupancy on genes is differentiation stage-dependent in both binding intensities and binding regions, indicating a shift in the regulatory mechanisms required for the entire program of osteoblastogenesis. By coupling genome-wide Runx2 binding with gene expression profiling, we have identified new Runx2 targets that were validated for functional activities of Runx2 binding in both promoter and non-promoter regions. Our study of Runx2 genome-wide occupancy establishes a foundation for future investigation of the Runx2-controlled regulatory network during bone formation and homeostasis.

4Runx2 is a known master activator of bone formation, but thus far only a small number of osteoblast-specific target genes have been characterized. To identify genome-wide occupancy of Runx2 in osteoblast lineage cells, we performed ChIP-Seq using a Runx2-specific antibody at stages during the in vitro differentiation of MC3T3-E1 preosteoblasts (Figure 1). This in vitro cell model recapitulates in vivo osteoblast differentiation and therefore was used for ChIP-Seq studies. Alkaline phosphatase activity, an early osteoblast marker, increased as cells proceeded from proliferation (day 0) to matrix deposition (days 9 and 21) and decreased upon mineralization (day 28), visualized by Von Kossa staining (Figure 1A). Runx2 mRNA and protein levels significantly increased during the initial stage of differentiation. While the mRNA levels reached steady state, Runx2 protein levels declined during late mineralization (Figure 1B). Osteoblast phenotypic markers, including the transcription factor Osx/Sp7 (osterix/Sp7 transcription factor), a marker of committed osteoprogenitors, the extracellular matrix protein Col1a1 (collagen type I alpha 1), and the specialized mineral binding proteins Bsp/Ibsp (bone sialoprotein/integrin-binding sialoprotein) and Ocn/Bglap2 (osteocalcin/bone gamma-carboxyglutamate (gla) protein 2), displayed expression patterns consistent with the progression of osteoblastogenesis (Figure 1C).

5We next isolated and sequenced DNA from chromatin bound by Runx2. Sequence reads from Runx2 ChIPs and input controls were mapped to the mouse genome. Statistically significant enrichments of Runx2 were identified by MACS (Model-Based Analysis of ChIP-Seq) (Additional file 1). RefSeq annotations were used to assign Runx2 enrichments to categories of genomic locations (with non-overlapping definitions for transcriptional start site (TSS), promoter, exonic, intronic, transcription termination site (TTS), upstream and intergenic regions (see Figure 2 for details; Additional file 2). As Runx2 protein levels changed, the total number of Runx2 peaks changed accordingly (Figure 2A, left panel). The overall distribution of Runx2 binding among the categories of genomic locations was relatively unchanged during differentiation (Figure 2A, right panel). Among genomic locations, Runx2 occupancy at promoters showed the greatest variation from 17.6% at day 0 to 8.8% at day 9 (Figure 2A, right panel). The majority of Runx2 binding occurred at intergenic and intronic regions (Figure 2A,B). However, when compared with a randomly sampled background distribution of genomic intervals (Additional file 3), Runx2 binding displayed preferential enrichment in a genic context, particularly at promoters and exons. One exception was under-represented Runx2 binding at intergenic regions when compared to nonspecific or random binding (Figure 2B). This relative enrichment of Runx2 occupancy in genic contexts was also observed in comparison with Runx2 motifs (Figure 2B). A de novo Runx2 motif (Figure 2C, top) was discovered using MEME on the highest-ranked (by P-value) 500 ChIP-Seq peaks and was then scanned using FIMO (version 4.7.0) over the mouse genome. Notably, analysis by MEME did not discover a significantly enriched secondary motif associated with Runx2 peaks. Despite the nearly random distribution of Runx2 motifs throughout the genome, Runx2 occupancy in differentiating osteoblasts is characterized by associations with promoters, exons, introns and other genic elements. These associations are perhaps due to epigenetic factors, including chromatin conformation and accessibility, along with co-factors, and suggest that the presence of a Runx motif does not necessarily indicate the physical association of Runx2. Interestingly, the distribution of Runx2 occupancy among classes of genomic elements is similar to that of CTCF (Figure 2B), a ubiquitous transcription factor that exhibits a broad spectrum of DNA binding in many cell lines.

6The de novo Runx2 motif (Figure 2C, top) was scanned over all ChIP-Seq peaks called by MACS from ChIP-Seq reads collected from day 9 MC3T3 cells, and a histogram of the distance between the peak summit and the highest-scoring motif instance was collected. The distribution in Figure 2D is characterized by a sharp peak for Runx2 motifs at the summits of ChIP-Seq peaks. Using TOMTOM it was determined that the de novo Runx2 motif was similar (E-value = 2.6 × 10-4, q-value = 2.6 × 10-4) to the JASPAR MA0002.2 RUNX motif (Figure 2C, bottom) and both motifs share the core TGTGGT sequence with the known Runx2 binding consensus.

7We grouped Runx2 peaks based on the presence or absence of Runx2 binding at specific time points during differentiation (Figure 3A, off/on). The resulting seven distinct clusters reflected the dynamics of Runx2 binding in relation to the progression of osteoblastogenesis (Figure 3A). In Figure 3B, the mean ChIP-Seq read densities (that is, peak intensities) of the clustered Runx2 peaks were plotted to compare their relative enrichments.

8The seven clusters differed in the numbers and intensities of Runx2 peaks (Figure 3A). The largest group was cluster 6, which exhibited the presence of peaks primarily at day 9. This finding is consistent with day 9 representing committed osteoblasts with the highest amount of cellular Runx2 protein (Figure 1B) and therefore the greatest number of Runx2 peaks (Figure 3B). The strongest Runx2 peak intensities were found in cluster 1 reflecting regions that were bound by Runx2 constitutively throughout osteoblast differentiation (Figure 3B). The weakest Runx2 peak intensities among the three time points occurred at the peaks in cluster 7, with marginal enrichment of Runx2 at day 0. It is noteworthy that the peaks in cluster 4 have the second highest peak intensities at both days 9 and 28 (matrix deposition and mineralization stages; Figure 3B). Cluster 5, like clusters 1 and 4, exhibited the highest peak intensities on day 28, indicating their importance in maintaining osteoblast phenotype.

9Runx2 binding in each cluster was further examined for distribution preferences of peaks in different genomic regions, in contrast to the genome-wide distribution of random 100 bp DNA fragments (detailed in Additional file 3; Figure 3C). The random DNA fragments (grey bars in Figure 3C) are distributed mainly in intergenic, intronic, and upstream regions. When compared to the random control, we observed that the distribution of Runx2-enriched peaks was biased towards gene regions (exons, introns, promoters, TTS regions, upstream regions). In contrast, Runx2 binding in the intergenic regions is lower than the random control (Figure 3C). In promoters and exons, all clusters showed the highest enrichment over random binding, suggesting strong regulation by Runx2 at these genomic regions.

10To further explore the relationship between Runx2 peaks and peak-associated genes in osteogenic differentiation, we performed functional annotations for the peaks in the seven clusters using GREAT (Genomic Regions Enrichment of Annotations Tool; Figure 3D). Gene Ontology (GO) terms associated with the largest clusters are shown in Figure 3D. Runx2 binding in cluster 1 yielded GO terms of general biological processes such as protein folding and RNA metabolism (detailed in Additional file 4). Cluster 4 peaks at days 9 and 28 were frequently related to the GO terms of osteoblast differentiation, bone developmental processes, and osteogenic signaling pathways. These terms often included differentiation-related and well-known Runx2 target genes, such as Runx2, Bsp/Ibsp, and Osx/Sp7 (Additional file 4). Similarly, cluster 6 (day 9) peaks often associated with bone formation and extracellular matrix organization. The annotation of other smaller clusters is shown in Additional file 5. For examples, cluster 3 peaks were associated with apoptosis, programmed cell death, and DNA damage; cluster 7 containing peaks found only in day 0 was linked with negative regulation of cell cycle control and the phenotypes of non-osseous mesenchyme-derived cells (Additional files 4 and 5). These functional annotations associated with Runx2 peaks are generally consistent with the progression of MC3T3-E1 differentiation, supporting a temporal transcription network programmed by Runx2.

11To discover previously unknown Runx2 target genes, we first determined the Runx2 binding patterns of well-known Runx2 target genes found in differentiation cluster 4; for example, Runx2, Osx/Sp7, and Ocn/Bglap2 (Figure 4A; Additional file 6), and Bsp/Ibsp (Additional file 7). For these genes Runx2 binding was distributed in promoters, the gene body (introns/exons) and sites distal from the gene body. Furthermore, Runx2 enrichment increased in the loci of these genes during osteoblast differentiation at matrix and mineralization stages. We then examined the genes associated with cluster 1 peaks (Additional file 8), and identified genes, including Ezh2 (enchancer of zeste homolog 2) (Figure 4B), with Runx2 binding profiles that displayed a ubiquitous but increasing level of Runx2 occupancy. Ezh2 is a component of PRC2 that epigenetically regulates gene expression by methylating histone H3 lysine 27 and was recently found to be involved in commitment of mesenchymal stem cells towards the osteoblast lineage. The 5′ proximal promoter region of Ezh2 is bound by Runx2 (Figure 4B; Additional file 6) and shows an increase in reads (occupancy) during differentiation. When Ezh2 mRNA levels were measured during osteoblast differentiation, we found that the highest expression occurred in proliferating MC3T3-E1 cells, when the lowest amount of Runx2 binding was observed (Figure 4C). The striking decrease in Ezh2 mRNA levels with increased Runx2 binding suggests that the transcription of Ezh2 is potentially regulated by Runx2.

12To characterize the Runx2 binding profiles in genes that are transcriptionally regulated by Runx2, we performed gene expression profiling at day 9 in MC3T3-E1 cells treated with control (Scr) short hairpin RNA (shRNA) or a Runx2-specific shRNA (shRunx2) that knocks down Runx2 expression. Runx2 protein levels decreased by 80% in cells treated with shRunx2 (Additional file 9). This knockdown in turn inhibited expression of differentiation marker genes and osteoblastogenesis as demonstrated by decreased Alp staining (Additional file 9). We found 159 genes whose expression was responsive to Runx2 knockdown, with |log2(IshRunx2/IScr)| > 1.5, and a false discovery rate (FDR) <0.05, where I is the measured normalized probe-set intensity. These genes included 115 up- and 44 down-regulated genes (Figure 5A). The 15 genes most responsive to Runx2 knockdown (Table 1) included the well-defined Runx2 targets Ocn/Bglap2, Bsp/Ibsp, and Mmp13, which are known to be activated by Runx2. Notably, the genes most upregulated by shRunx2 have not been characterized as Runx2 targets (Table 1) except for Usp18, which encodes a protein involved in the ubiquitin degradation pathway.

13We characterized the relationship between Runx2 occupancy and genes affected by shRunx2 knockdown, compared to non-responsive genes. The number of peaks, genomic distribution of peaks, and fold change of peak signals were compared among the gene groups at days 0 and 9 (Figure 5B-D). Genes that were downregulated by shRunx2 at both days 0 and 9 had more Runx2 peaks when compared to the control non-responsive genes (Figure 5B). It should be noted, however, that shRunx2 downregulated genes are longer than the upregulated genes (Additional file 10). This finding was also observed when we included the Runx2 binding profiles at day 28 for Runx2-regulated genes (Additional file 10). In contrast, the shRunx2 upregulated genes that appeared on day 9 had fewer Runx2 peaks compared to control genes (Figure 5B). There were also significant differences in overall peak distribution between shRunx2 responsive genes and control (Figure 5C). Genes that were downregulated tended to exhibit more intronic and intergenic enrichment of Runx2 peaks, while shRunx2 upregulated genes were strongly enriched in intergenic but reduced in intronic binding. We further examined the distribution of Runx2 peaks in up- and down-regulated genes as a function of changes in Runx2 binding during differentiation from days 0 to 9 (Figure 5D). The fold change in day 9/day 0 peak signals showed increased Runx2 binding predominantly at upstream and promoter regions for the shRunx2 downregulated genes; in shRunx2 upregulated genes, Runx2 binding diminished in the introns and showed no significant change in promoter regions. Therefore, shRunx2 downregulated and upregulated genes exhibited distinct preferences for Runx2 binding in genomic loci as reflected by peak distributions in Figure 5B,C, and in relation to differentiation (Figure 5D).

14To complement the above analysis examining the genome-wide distribution of Runx2 responsive peaks, we used the PeaksToGenes program to determine the enrichment of day 9 Runx2 signals at defined intervals within and surrounding the gene bodies (Figure 5E). For shRunx2 downregulated genes, Runx2 binding had the strongest enrichment surrounding TSSs, which includes proximal promoter, 5′ UTR, exons and introns all within five contiguous deciles (Figure 5E; Additional file 11). This analysis demonstrated that, of genes affected by shRunx2, 66.1% of upregulated genes have a low level of Runx2 binding (cluster V in Figure S5B in Additional file 11) but most downregulated genes (69.1%) have higher level intensities of Runx2 binding (Figure S5A,B in Additional file 11). This finding shows that Runx2 responsive genes at day 9 (shRunx2/Scr) are primarily regulated by Runx2 surrounding TSSs.

15As another computational analysis, EMBER (Expectation Maximization of Binding and Expression pRofiles) was used to relate measured changes in gene expression to the spectrum of Runx2 occupancy observed during osteoblast differentiation (Figure 3B,C; Additional file 12). In analogy with discovering a sequence motif from a collection of functionally related DNA sequences, EMBER optimizes an ‘expression pattern’ from a collection of genes related by patterns of transcription factor binding and uses this motif to determine which genes are regulatory targets of the transcription factor (details of EMBER are summarized in Additional file 3). Using this approach, we discovered and compared expression patterns from different sets of Runx2 binding regions (Figure 3). The Runx2 peaks were partitioned into 42 subsets of Runx2 binding regions (7 clusters with 6 genomic location categories) and it was observed that all groups of Runx2 binding show a correlative relationship to gene expression during osteoblast differentiation (Additional file 12). We noted, however, that intronic Runx2 binding regions - for all time-dependent patterns of Runx2 occupancy - are less informative than other groups, indicating that intronic Runx2 binding may be less useful than binding at other class elements as a predictor of gene regulation.

16Finally, downregulated and upregulated genes were compared on the basis of their evolutionary conservation (Figure 5F). Functional genomic elements are often characterized by conservation, which has been used to guide the prediction of transcription factor binding sites. This finding has been shown to distinguish functionally verified from un-verified binding sites in a large-scale study of transcription factor binding site function on human promoters. For each ChIP-Seq peak, a Runx motif was used to identify the single most likely Runx2 binding site and the mean phyloP score (for conservation among 30 vertebrate species) across the binding sites was computed. For each gene, phyloP scores were averaged among all ChIP-Seq peaks that were associated with individual genes to give an average measure of Runx2 conservation. We found that genes downregulated by shRunx2 had Runx2 binding sequences that were more conserved than those in upregulated or non-responsive genes (Figure 5F), suggesting that Runx2 regulation may be more evolutionarily conserved in genes that are activated by Runx2 during osteoblastogenesis.

17Taken together, the complementary methods described above suggest that Runx2 employs different mechanisms to regulate gene expression: 1) shRunx2 downregulated genes show increased Runx2 binding at promoter and far upstream regions; and 2) based on EMBER, intronic binding does not imply Runx2-mediated gene regulation to the same degree as Runx2 binding at promoter, exon and upstream regions. Thus, during the normal course of osteoblast differentiation, Runx2-activated genes (shRunx2 downregulated) are regulated through both promoter and non-promoter regions; and regulation of Runx2 repressed genes (shRunx2 upregulated) also occurs in promoter and other genomic regions, but with less Runx2 binding.

18Among potential Runx2 targets, we identified Tnfrsf19 (tumor necrosis factor receptor superfamily, member 19), which is involved in bone formation as a Wnt-responsive regulator of mesenchymal stem cell commitment to osteoblastic lineage. Tnfrsf19 exhibited enrichment of Runx2 binding to the promoter as well as intronic regions (Additional files 6 and 13). During differentiation of osteoblasts, Tnfrsf19 mRNA levels increased dramatically more than 100-fold from day 0 to day 9 (Figure S7A in Additional file 13). Consistent with Affymetrix data, depletion of Runx2 resulted in significant decreases in Tnfrsf19 expression, indicating direct Runx2 regulation (Figure S7B in Additional file 13). We found that Runx2 occupancy was increased at intronic and promoter regions of the Tnsrsf19 locus during osteoblast differentiation (Figure S7C in Additional file 13).

19Adamts4 and Crabp2 were selected for further analyses based on their responsiveness to shRunx2 (with fold change cutoff of 1.3; Additional file 14) and their Runx2 binding patterns predominantly in non-promoter regions. To test the functionality of Runx2 binding at these putative regulatory regions, we cloned non-promoter Runx2 binding regions and measured transcriptional activity by luciferase reporter assay. Adamts4 is expressed in osteoblasts and osteocytes and encodes an enzyme that degrades aggrecan. Runx2 exhibits multiple peaks across the Adamts4 locus, with increased occupancy during osteoblast differentiation (Figure 6A; Additional file 6). Adamts4 expression during osteoblast differentiation was increased at day 9 and remained steady to day 28 (Figure 6B). Knockdown of Runx2 suppressed the expression of Adamts4 (Figure 6C) by 40%, supporting Adamts4 as a direct target of Runx2 during osteoblastogenesis. We characterized the functional activity of two prominent Runx2 binding regions, peak A in intron 1 and peak B in the last exon (Figure 6A). The peak A region increased luciferase activity in MC3T3-E1 cells by over 20-fold. In contrast, peak B functioned as a suppressor of luciferase activity (Figure 6D). These results indicate that the two Runx2 sites can function as a positive and negative regulator of Adamts4; however, the weaker, negative regulation by peak B on the luciferase reporter may be due to the lack of a native chromatin context. The large increase in luciferase activity observed from peak A is consistent with a previous study that demonstrated Runx2 upregulates ADAMTS4 in human SW1353 chondrosarcoma cells. Here we established in osteoblasts that Runx2-mediated upregulation of Adamts4 can occur at non-promoter regulatory elements, and is not restricted to the proximal promoter region as previously shown.

20Crabp2 is a cytoplasmic retinoic acid binding protein previously reported to be upregulated during osteoblastogenesis. We observed that Runx2 constitutively occupies the Crabp2 locus in the first intron, while binding increases upstream and downstream of the Crabp2 gene body during differentiation (Figure 6E; Additional file 6). Crabp2 was upregulated during MC3T3-E1 cell differentiation and knockdown of Runx2 decreased Crabp2 mRNA level (Figure 6F). In proliferating cells, peak regions C, D and E demonstrated minimal luciferase reporter activity. However, in differentiating cells, all peaks exhibited a significant increase of luciferase activity, with peak region D showing a more than 30-fold activation (Figure 6G). The knockdown of Runx2 reduced luciferase activity (Figure 6H), further supporting the function of these regions in mediating Runx2 regulation of Crabp2.

21These findings of novel genes bound and regulated by Runx2 through different types of genomic elements support an emerging concept that non-promoter elements can regulate gene transcription. Our results also indicate that Runx2 mediates complex fine-tuning of gene expression in osteoblasts by both activating and repressing regulatory elements that are located in non-promoter regions.

22Through comprehensive genomic analysis of Runx2 by ChIP-Seq, we describe widespread Runx2 binding throughout the genome of differentiating osteoblasts. In addition to Runx2 interaction at promoters, we find Runx2 binding in non-promoter regions regulating novel targets that are silenced or expressed at different stages of osteoblast differentiation. Runx2 peaks cluster into temporal and functional categories associated with genes in a broad range of cellular programs, including bone development, negative regulation of proliferation, active matrix formation and mechanisms for mineral deposition that reflect the progression of osteogenesis. Our data have identified new Runx2-regulated genes (Tnfrsf19, Adamts4, Crabp2, and Ezh2) that have established roles in bone formation, and more importantly, have extended the understanding of Runx2-mediated gene regulation to a broader range of cellular functions during osteoblast differentiation.

23Time-dependent Runx2 binding patterns underlie the dynamic gene regulation by Runx2 during osteoblastogenesis. We identified a large number of peaks, consistent with the increasing protein levels of Runx2 from the early osteoprogenitor to the mature osteoblast/osteocyte. The results from binary clustering, together with subsequent GO term analyses by GREAT, identified many categories associated with skeletal development and bone homeostasis. These findings support our initial hypothesis that distinct binding patterns of Runx2 at different stages of osteoblastogenesis have novel functional implications.

24Our clustering analysis partitioned Runx2 peaks into two main categories: less variable steady-state binding (cluster 1, ubiquitous peaks) and more dynamic binding groups (clusters 2 to 7, stage-specific clusters). Steady-state binding of Runx2 persists from proliferative pre-osteoblasts to differentiating osteoblasts with binding signals plateauing during matrix deposition and mineralization stages. Peaks in this category represent genes related to housekeeping processes such as protein folding, negative regulation of mitotic cell cycle, and mRNA catabolism and processes not previously related to Runx2 (Additional file 4). The GO terms associated with stage-specific clusters include negative regulators of other cell lineages (that is, fat and smooth muscle cells) as well as positive regulators of osteogenesis. Thus, dynamic Runx2 binding primes, enhances and stabilizes the osteoblast phenotype as well as suppresses non-osteoblast lineages. Many Runx2 bound genes in cluster 4 (days 9 & 28, differentiation) have been demonstrated to contribute to in vivo bone formation. Thus, the genomic profiling of Runx2 binding in our in vitro model system is consistent with the known properties of Runx2 in bone formation. More importantly, our profiling reveals pathways previously unknown to be controlled by Runx2 underlying biological mechanisms of general cellular processes.

25Only a small proportion of sequence-specific transcription factors, such as Myc, have narrow distributions of binding around proximal promoters of genes. In contrast, non-promoter binding is recognized for other transcription factors, including STAT1, RUNX1, ERα, CTCF, and HNF4α. In a previous study, overexpression of Runx2 in prostate cancer cells revealed extensive non-promoter binding. In our study, endogenous Runx2 binding across the genome was characterized during osteoblastogenesis. We found that over 70% of Runx2 occupancy was localized to non-promoter regions (intergenic, intron, exon, TTS, and upstream regions that constitute the bulk of the genome) during the differentiation of osteoblasts. From days 0 to 9 there was a two-fold increase in the number of non-promoter peaks, indicating a functional association with the differentiation process. Runx2-dependent regulation through non-promoter peaks around the Adamts4 and Crabp2 genes provided direct evidence that non-promoter binding of Runx2 controls gene expression.

26Although we have demonstrated the importance of non-promoter Runx2 binding events, Runx2 peaks at promoter regions have critical regulatory roles as well. In our data, Runx2 occupancy has the highest enrichment at promoter regions when compared with other genomic locations. Examples of this regulatory mode can be seen in some well-characterized Runx2 targets such as Bsp and Ocn, in line with established evidence that Runx2 can regulate these genes in a promoter-dependent manner. The genes downregulated upon Runx2 silencing also displayed a clear enrichment of Runx2 signal at promoter regions (Figure 5B), exemplified in Tnfrsf19.

27Gene expression in different biological settings is influenced by higher order three-dimensional chromatin complexes that involve looping of promoter and non-promoter elements, blurring the distinction of defined regulatory regions. It is plausible that some Runx2 peaks in promoter and non-promoter regions may serve as nucleation sites for modifications of chromatin structures necessary for gene expression. Furthermore, transcription factors can interact with RNA polymerase II, CTCF, and other factors via higher-order chromatin conformations. During osteoblastogenesis, Runx2 exhibited a distribution pattern among genomic elements similar to that of CTCF, suggesting that, like CTCF, Runx2 may have a functional role that extends beyond direct regulation of transcription. Consistent with this idea, Runx2 is able to form discernible foci associated with the nuclear matrix: a nuclear framework for organizing higher order chromatin structures. Runx2 truncated of the NMTS (nuclear matrix-targeting signal) domain results in diminished nuclear matrix association and disrupted expression of Runx2 target genes. It is also noteworthy that the genes upregulated upon Runx2 knockdown have preferential Runx2 binding in intergenic regions, indicating that distal elements may have a regulatory role through long-range interactions. In addition, Runx2 binding at distal regions may contribute to chromatin remodeling through interacting with chromatin-modifying enzymes, as has been well documented at regulatory elements in osteoblasts and other cell models.

28Runx2 displays complex binding patterns similar to other lineage-specific transcription factors, such as PPAR-γ, MyoD, and GATA3. By systematic annotation of Runx2 peaks, multiple integrative analyses of gene expression combined with Runx2 binding profiles and direct experimental validation of individual targets, we defined Runx2 binding with biological outcomes during osteoblast differentiation. These analyses revealed that, for a small set of genes, the enrichment and binding patterns of Runx2 were indicators of gene expression. Genes that have decreased expression in the absence of Runx2 (by shRunx2 treatment) have a greater number of Runx2 peaks and greater fold change of peak signal at promoter and non-promoter regions when compared to non-responsive and gene length-matched controls. In contrast, genes upregulated by Runx2 knockdown tend to have fewer Runx2 peaks and smaller relative fold changes in peak signals when compared to controls. One notable finding that arose from our analysis was that genes that were downregulated in the absence of Runx2 had both a greater evolutionary conservation of Runx2 binding sites and tended to be longer than non-responsive and upregulated genes. It is unclear why this particular subset of genes (that is, genes normally activated by Runx2) would retain these features throughout evolution; however, this is an interesting point for future investigations. Although we demonstrated that Runx2 binding influences gene expression, a proportion of Runx2 peaks were found to have no direct function in transcriptional control of genes. This may be due to efficient but not complete knockdown of Runx2 by viral-mediated shRNA. Alternatively, similar to other transcription factors, many binding regions were found to be nonfunctional in transactivating luciferase reporters. This finding suggests that binding of transcription factors may have a distinct function other than direct control of gene expression, consistent with previously described non-transcriptional functions and genome-organizing capabilities of Runx2.

29Our findings provide a new level of understanding of the Runx2-mediated transcription program as defined by genome-wide Runx2 binding essential for osteoblastogenesis. Our data support that Runx2 functions at promoter and non-promoter regions at both previously known and novel targets. The impact of our study examining the global occupancy of endogenous Runx2 in differentiating osteoblasts sets a framework for novel mechanisms underlying bone biology.

30The calvaria-derived preosteoblast cell line MC3T3-E1 (Subclone 4) was obtained from ATCC (Manassas, VA, USA) and maintained in ascorbic acid-free alpha-MEM (Hyclone, Novato, CA, USA) supplemented with 10% fetal bovine serum (FBS; Hyclone), 2 mM L-glutamine, 100 U/ml of penicillin and 100 μg/ml streptomycin (Pen/Strep, Invitrogen, Carlsbad, CA,). To induce osteogenic differentiation, complete alpha-MEM was supplemented with 280 μM ascorbic acid and 10 mM beta-glycerophosphate (Sigma Aldrich, St. Louis, MO, USA). Cells were maintained at 37°C in a humidified 5% CO2 environment and media replaced every 2 to 3 days for the duration of all experiments.

31Total RNA was isolated using Trizol reagent (Invitrogen) according to the manufacturers’ specifications. Total cellular RNA treated by DNaseI (Zymo, Irvine, CA, USA) was primed with random hexamers and reverse transcribed into cDNA using Superscript First-strand cDNA Synthesis Kit (Invitrogen) according to the manufacturer’s instructions. Gene expression was determined by quantitative real time PCR (qPCR) using iQ™ SYBR Green PCR Master Mix (BioRad, Hercules, CA, USA) in an ABI Prism 7300 thermocycler (Applied Biosystems, Foster City, CA, USA). For each gene, the expression level was normalized to that of Hprt1 using 2-ΔΔCT method. Experiments were performed in triplicate and results are presented as mean values ± SEM. Primers for qPCR reactions were designed by FoxPrimer, and are available in Additional file 15.

32Lentiviruses carrying Runx2-shRNA and previously described control Scramble-shRNA were used to infect MC3T3-E1 cells. Infected cells were subsequently detected by green fluorescent protein, sorted and grown to 90% confluency followed by osteogenic differentiation for 9 days. Knockdown experiments were performed in three biological replicates.

33Microarray samples were handled following the manufacturers’ recommended protocols (Affymetrix, Santa Clara, CA, USA). Briefly, RNA isolated from MC3T3-E1 cells (day 0 and 9 scramble shRNA, and day 9 Runx2 shRNA) were reversely transcribed into cDNA using WT Expression Kit (Ambion, Austin, TX, USA), labeled and fragmented with GeneChip WT Terminal Labeling and Controls Kit. Labeled cDNAs were then hybridized to GeneChip Mouse Gene 1.0 ST Array rev.4 using a GeneChip Hybridization, Wash, and Stain Kit. Hybridization signals were obtained by GeneChip Scanner (Affymetrix). Microarray data were analyzed using Bioconductor (version 2.11) packages affy and limma in R (version 2.15.1). Briefly, after performing RMA (robust multichip average) normalization of microarray expression levels and filtering, differential expression was detected using a Bayesian moderated t-test. The Benjamini-Hochberg FDR was applied to correct for multiple testing. The Affymetrix expression profiles were annotated to RefSeq genes. These expression data have been deposited in the Gene Expresion Omnibus (GEO) database under accession number GSE53982.

34Nuclei extracts were prepared from MC3T3-E1 cells using a protocol modified from Dignam et al.. Primary antibodies and dilutions were: mouse anti-Runx2 monoclonal (Clone 8G5, MBL International, Woburn, MA, USA; 1:1,000); rabbit anti-H3 3H1 monoclonal (Cell Signaling, Danvers, MA, USA; 1:2,000). HRP-conjugated secondary antibodies and dilutions were: goat anti-mouse IgG (Santa Cruz, Dallas, TX, USA; 1:3,000); goat anti-rabbit IgG (Santa Cruz; 1:3,000). Detection of HRP was performed using a Western Lightning Plus Kit (Perkin Elmer, Waltham, MA, USA) followed by exposure onto Biomax Light film (Kodak, Rochester, NY, USA).

35At days 0, 9, and 28 of differentiation, approximately 1 × 108 MC3T3-E1 cells were washed with PBS (phosphate-buffered saline) and then fixed on a plate with 1% formaldehyde for 8 minutes to crosslink DNA-protein complexes. The fixed cells were washed with ice-cold PBS, harvested, and pelleted. Nuclei extraction was performed using a protocol modified from Dignam et al.. Isolated nuclei were sonicated using a Misonix S-4000 ultrasonic sonicator to obtain sheared chromatin ranging from 0.2 kb to 0.6 kb. Sheared chromatin was used for immunoprecipitation with Runx2 antibody (M-70, Santa Cruz) or immunoglobulin G (IgG) (12-370, Millipore, Billerica, MA, USA) followed by purification using Protein-G Dynabeads (Invitrogen). Precipitated chromatin was washed with solutions of increasing salt concentration, and eluted and subsequently uncrosslinked at 65°C. DNA was recovered by phenol/chloroform/isoamyl alcohol extraction followed by ethanol precipitation. Libraries of purified DNA were generated using Illumina SR adapters (Illumina, San Diego, CA, USA) following the manufacturer’s manual. DNA libraries were selected for inserted fragments of 200 ± 50 bp, and single-end 36 base reads were generated on an Illumina Genome Analyzer II at the UMASS Deep-Seq core facility. Base calls and sequence reads were generated by Illumina CASAVA software (version 1.6; Illumina). Two independent biological repeats of Runx2 ChIP-Seq libraries were prepared for each time point, and two Input libraries were prepared with sonicated DNA from day 9 MC3T3-E1 cells.

36Single-end 36 base sequences from Runx2 ChIP-Seq and input libraries were mapped to the mouse genome (assembly mm9) using Bowtie (version 0.12.8). Runx2 peaks and read counts were determined by MACS (version 1.4.1) using default settings. Runx2 peaks that were significant at the P < 10-10 level were retained for subsequent analysis. For each of the three time points in our data, read counts were normalized to 10 million reads. The UCSC genome browser was used to visualize Runx2 peaks.

37Runx2 binding regions were classified on the basis of genomic location categories and annotated to known RefSeq genes. Runx2 peaks were grouped into seven clusters based on the presence or absence of a peak. Runx2 peaks in each cluster were analyzed for GO terms using GREAT (version 2.0.2), using default association rules between ChIP-Seq peaks and annotated genes. The conservation of Runx2 motifs associated with the genes responsive to shRunx2 were determined using phyloP. Statistical significance was determined using Kolmogorov-Smirnov test.

38PeaksToGenes (Additional file 3) was used to test Runx2 binding in relation to the genes from microarray analysis. PeaksToGenes defines genomic intervals relative to all RefSeq genes, and in each window uses a non-parametric Wilcoxon rank sum test to calculate the probability and binding frequency. The comparisons were individually made between each responsive group and the non-responsive group to Runx2 shRNA. Runx2 binding and Runx2-mediated transcription control were also evaluated by EMBER. Analogous to discovering a sequence motif from functionally related DNA sequences, EMBER optimizes an expression pattern from a collection of genes’ expression data related by profiles of transcription factor binding and uses this information to determine which genes are potential regulatory targets of the transcription factor. For a detailed description of the PeaksToGenes and EMBER analysis, please refer to Additional file 3.

39Functional genomic elements can be characterized by evolutionary conservation and have been shown to distinguish functionally verified from unverified transcription factor binding sites on human promoters. For each ChIP-Seq peak (combined from day 0, 9, and 28 datasets), a Runx motif was used to identify the single most likely Runx2 binding site and the mean phyloP score (for conservation among 30 vertebrate species) across the binding sites was computed. In order to compare conservation among genes that fell into different regulatory groups (shRunx2 downregulated, upregulated and non-responsive genes) based on our expression microarray measurements, phyloP scores were averaged among all ChIP-Seq peaks that were associated with individual genes to give an average measure of Runx2 conservation for each gene. To minimize the effect of gene length on Runx2 binding analysis, we compared the conservation of shRunx2-downregulated and -upregulated genes to randomly sampled, length-matched non-responsive genes. Statistical significance was determined by pair-wise conservation comparisons using Kolmogorov-Smirnov test.

40The raw sequences and peak-related files in BED and WIG formats representing processed data have been deposited in the GEO database under accession number GSE54013.

41Selected Runx2 peak regions were cloned with MluI combined with BglII or XhoI into a pGL2-SV40-Luc reporter (Promega, Fitchburg, WI, USA). Primers used in cloning are listed in Additional file 16. A pGL2-SV40-Luc reporter with minimum SV40 promoter was used as mock control. Reporter plasmids with Runx2 peak regions and pGL2-SV40-Luc empty vector were co-transfected with pcDNA3.1-Runx2-WT or pcDNA3.1-EV into MC3T3 cells. Transfected cells were then differentiated for 7 days before luciferase activities were determined by Dual-Glo Luciferase Assay Kit (Promega). A second set of luciferase assays was done in differentiating MC3T3-E1 cells stably infected with a doxycycline-inducible pLKO-puro-Tet-on-Rx2shRNA lentiviral vector. This vector was constructed by re-cloning a previously validated Runx2 shRNA sequence to Tet-pLKO-puro plasmid (catalog number 21915, Addgene, Cambridge, MA, USA). Luciferase activities from MC3T3-E1 cells treated with or without 2.5 μg/ml doxycycline were examined at day 7 after differentiation (Figure 6).

42bp: base pair; ChIP: chromatin immunoprecipitation; FDR: false discovery rate; GEO: Gene Expresion Omnibus; GO: Gene Ontology; qPCR: quantitative real time polymerase chain reaction; shRNA: short hairpin RNA; TSS: transcription start site; TTS: transcription termination site; UTR: untranslated region, SEM, standard error of mean.

43The authors declare that they have no competing interests.

44HW and JARG conceived and designed the experiments. HW performed the experiments. TWW, JRD, and HW performed bioinformatics data analysis. HW, JARG, TWW, and JRD interpreted the data with insights from PWLT, JBL, JLS, AVW, and GSS. HW, TWW, JARG, JRD, PWLT, JBL, and JLS wrote the paper with help from AVW and GSS. All authors read and approved the final manuscript.