The RUNX2 Cistrome in Osteoblasts

Background: RUNX2 is indispensable for mature osteoblast differentiation and function. Results: During osteogenic differentiation, RUNX2 along with C/EBPβ bind throughout the genome near genes that are essential for the phenotype. Conclusion: Differentiation leads to a restricted RUNX2 cistrome that correlates with changes in gene expression. Significance: Broad genomic localization of RUNX2 during osteoblastogenesis highlights its role as a dynamic lineage-determining factor. RUNX2 is a transcription factor that is first expressed in early osteoblast-lineage cells and represents a primary determinant of osteoblastogenesis. While numerous target genes are regulated by RUNX2, little is known of sites on the genome occupied by RUNX2 or of the gene networks that are controlled by these sites. To explore this, we conducted a genome-wide analysis of the RUNX2 cistrome in both pre-osteoblastic MC3T3-E1 cells (POB) and their mature osteoblast progeny (OB), characterized the two cistromes and assessed their relationship to changes in gene expression. We found that although RUNX2 was widely bound to the genome in POB cells, this binding profile was reduced upon differentiation to OBs. Numerous sites were lost upon differentiation, new sites were also gained; many sites remained common to both cell states. Additional features were identified as well including location relative to potential target genes, abundance with respect to single genes, the frequent presence of a consensus TGTGGT RUNX2 binding motif, co-occupancy by C/EBPβ and the presence of a typical epigenetic histone enhancer signature. This signature was changed quantitatively following differentiation. While RUNX2 binding sites were associated extensively with adjacent genes, the distal nature of the majority of these sites prevented assessment of whether they represented direct targets of RUNX2 action. Changes in gene expression, however, revealed an abundance of genes that contained RUNX2 binding sites and were regulated in concert. These studies establish a basis for further analysis of the role of RUNX2 activity and its function during osteoblast lineage maturation.

numbers of hypertrophic chondrocytes, and an inability to form bone (20 -22). Interestingly, loss of a single RUNX2 allele in mice results in hypoplastic clavicles and delayed closure of the fontanelles, and in cleidocranial dysplasia in humans (14,23). It is also apparent that in addition to SOX9, RUNX2, and OSX, osteoblast differentiation and function are also regulated by the temporal and overlapping expression of a myriad of additional transcription factors that include MSX2, TWIST1, SATB2, NFATc1, ATF4, STAT1, C/EBP isoforms, and multiple nuclear receptors as well (2). Many of these factors interact directly with RUNX2 and/or OSX either to facilitate or suppress transcription factor events that are important to osteoblast differentiation and/or function. Nevertheless, while many of the players are known, the overall temporal activities of each and their interactions with each other are not well understood.
RUNX2 represents one of three members of the RUNT domain family of transcription factors (RUNX-1, -2, and -3) that play central roles in cell fate determination, directing hematopoiesis, skeletogenesis and neurogenesis among others (17). Interestingly, RUNX1 also plays a role in osteoclast formation (24). Aberrant expression of these proteins can also lead to a variety of cancers (25,26). RUNX2 is targeted to the nucleus via a C-terminal domain (27,28). RUNX family members are also DNA-binding proteins that interact directly via a common CBF␤ heterodimer partner (17). The heterodimer's activities at target gene loci are controlled, however, via diverse mechanisms that modulate both RUNX2 abundance and its ability to bind to DNA and that mediate extensive post-translational modifications that include phosphorylation, acetylation, and ubiquitinylation (29,30). These modifications influence RUNX2 turnover, coregulator recruitment and relative transcriptional competence. This complexity in regulation highlights the relevance of the direct interactions that occur between RUNX2 and the many additional classes of factors that include cellular kinases, DNA-binding (C/EBP␤, SMADS, ␤-catenin) and non DNA-binding (HDACs) transcription factors, and regulatory factors that control chromatin structure and function as well. Temporal control of RUNX2 activity ensures the timely expression of gene networks that orchestrate bone cell development and function.
A number of genes are known to be directly regulated by RUNX2 (2). These targets include osteocalcin (Bglap2), osteopontin (Spp1), Fgf18, Rankl, Vdr, and Runx2 itself. The typical DNA response sequence located at these gene loci that acts in cis is comprised of 5Ј-T/A/C G T/A/C GG T/G-3Ј (31) and in a gene-specific context functions both to induce or suppress the expression of genes. A clear understanding of the network of genes that are directly regulated by RUNX2 has not been delineated however, due primarily to this factor's diverse biological role in osteoblast lineage cells as well as its complex regulation. A recent ChIP-chip analysis which coupled the binding of RUNX2 at a genome-wide collection of promoter regions in osteosarcoma cells to gene expression has defined a number of potential genetic sites of action of RUNX2 (32). In addition, an unbiased genome-wide ChIP-seq analysis using RUNX2-negative prostate cells engineered to express RUNX2 in a doxycycline-conditional manner suggests a limited RUNX2 cistrome that is enriched for the cognate RUNX2 binding site motif (33).
Despite these studies, an in-depth analysis of the RUNX2 cistrome in osteoblast lineage cells has not been described.
In the current study, we used the mouse pre-osteoblastic cell line MC3T3-E1 (POB) and it's in vitro differentiated osteoblast progeny (OB) to investigate the properties of the RUNX2 cistrome in osteoblast lineage cells and to explore the relationship between RUNX2 binding and gene expression. We show through ChIP-seq analysis that while the RUNX2 cistrome is extensive in POBs, it undergoes significant restriction upon differentiation to OBs. Many of the sites identified confirm those found in previous studies through single site analyses. Interestingly, C/EBP␤ was found to be a frequent genome-wide co-occupant at many sites of RUNX2 binding. Sites occupied by either RUNX2 or C/EBP␤ or both exhibited enriched levels of methylated histone H3K4 (H3K4me1) and acetylated histones H3 and H4, known epigenetic signatures that highlight authentic enhancers. Finally, we show that components of the RUNX2 cistromes correlate directly with both genes whose expression levels are unchanged as well as genes whose expression were altered during in vitro differentiation. Changes in RUNX2 binding are noted at genes whose activities are both increased and decreased during the process of differentiation.
Cell Culture and Differentiation-Mouse MC3T3-E1 cells (early passage line from Sudo, H. et al. (34)) were cultured in minimum Eagle's medium alpha modification supplemented with 10% heat-inactivated fetal bovine serum from Hyclone (Logan, UT) and 1% penicillin-streptomycin from Invitrogen to obtain POBs. Cells were grown to confluency and then changed into differentiation media (10 mM ␤-glycerophosphate and 50 g/ml ascorbic acid) and cultured for an additional 15 days, replenishing the media every 2-3 days, to prepare OBs. At the conclusion of differentiation, cells were stained with Alizarin Red and Von Kossa for mineral deposition and bone nodule formation as well as alkaline phosphatase (ALP) and an ALP enzymatic activity assay. For Alizarin Red staining, the cells were washed in PBS, fixed with 10% formaldehyde solution for 20 min at room temp, stained with Alizarin Red solution for 20 min (Alizarin Red S (Sigma-Aldrich) in ddH 2 O, pH 4.2), rinsed with tap water and dried. For Von Kossa staining, cells were washed with PBS, 5% silver nitrate solution was added for 30 min at room temp, rinsed with ddH 2 O, exposed to carbonateformaldehyde solution for 5 min, rinsed with tap water and dried. For ALP staining, cells were washed in PBS, then stained with Alkaline Phosphatase kit (Sigma, 86R-1KT) according to manufacturer's protocol. For ALP enzymatic activity assay, cells were washed in PBS twice, lysed for 30 min shaking at room temperature in Lysis Buffer (20 mM Tris, 0.5 mM MgCl 2 , 0.1 mM ZnCl 2 , 0.1% Triton X-100), scraped from plates and transferred to microfuge tubes for assay or storage at Ϫ80°C. The activity assay was performed according to the Alkaline Phosphatase Kit (Sigma 104-100). Samples were compared against a standard curve of diluted p-nitrophenol solution (Sigma N7660) in 0.02 N NaOH. Activity is expressed as International Units (IU) per liter (L), normalized to total lysate protein as assayed by Bradford. A Student's t test was completed for statistical significance compared with D0 (POB).
Gene Expression Analysis-POBs or OBs were first prepared from MC3T3-E1 cells as described above. RNA was isolated using the TRI-Reagent protocol (MRC, Cincinnati, OH) and amplification methods according Roche NimbleGen protocols. Briefly, 10 g of total RNA was DNase treated and reverse transcribed with the Superscript double stranded cDNA synthesis system (Invitrogen, Carlsbad, CA). Double stranded cDNA (dscDNA) was verified for integrity by an Agilent 2100 Bioanalyzer and qPCR analysis. dscDNA was then labeled with Cy3 as described by Roche NimbleGen. Labeled samples were hybridized to mouse (mm8) 385k microarrays (Roche NimbleGen) in triplicate. Extracted data were processed with the Arraystar 5 package from DNASTAR (Madison, WI) using a moderated t test and multiple testing correction, FDR (Benjimini-Hochberg). All data passing 95% confidence limits greater than 2-fold were included in up-or down-regulated data analysis. Data are displayed as log2 values (log2) as well as fold change values in supplemental Table S1. For all other qRT-PCR and validation, cells were differentiated for 15 or 30 days, 1 g of isolated total RNA was DNase treated, reverse transcribed using the High Capacity cDNA Kit (ABI, Foster City, CA), and then diluted to 100 l with RNase/DNase free water. qPCR was performed using primers specific to a select set of differentially expressed genes by Taqman analyses available in supplemental Table S1.
Chromatin Immunoprecipitation (ChIP) followed by Sequencing (ChIP-seq)-Chromatin immunoprecipitation was performed as described previously (35). Antibodies are listed in Reagents above. Briefly, samples were subjected to immunoprecipitation using either a control IgG antibody or the indicated experimental antibody (VDR, RXR, C/EBP␤, RUNX2, H3K5Ac, H3K9Ac, H3K4me1, H3K4me3, or H3K36me3). For differentiated cells, there was substantial mineralized matrix. To remove the matrix from the assay, the cell/matrix mix was subjected to 3 ϫ 15s pulses with a Polytron Homogenizer (Power Gen 125, Fisher Scientific) while in NCP 1 buffer (10 mM Hepes, pH 6.5, 10 mM EDTA, 0.5 mM EGTA, 0.25% Triton X-100). The resulting homogenate was centrifuged (200 ϫ g) through a 10% sucrose buffer at a 1:1 ratio (10 mM Tris-HCl, 15 mM NaCl, 60 mM KCl, 1 mM EDTA, 10% sucrose). The resulting cell pellet was re-suspended in NCP 2 buffer (10 mM Hepes, pH 6.5, 1 mM EDTA, 0.5 mM EGTA, 200 mM NaCl) and the remainder of the protocol was followed (35). The isolated DNA (or Input DNA acquired prior to precipitation) was then validated by quantitative real time PCR (qPCR) and further prepared for ChIP-seq analysis. ChIP-seq libraries were prepared using the NEBNext DNA sample prep kit (NEB, E6000L) with the Bioo NEXTflex ChIP-seq Barcodes (Bioo Scientific, Austin, TX, 514122) according to manufacturer's protocols. During the NEBNext prep, the Illumina adapters were replaced with the Bioo Scientific Barcoded adapters according to Bioo protocols. ChIP-DNA ligated libraries were purified with Agencourt AMPureXP Magnetic Beads (Beckman-Coulter, A63881). A pre-size selection PCR procedure was performed using Phusion polymerase, NEXTflex Primer Mix, and purified ligation product for 4 PCR cycles according to the Bioo Protocol. Libraries were size-selected using Invitrogen E-gels to 400 -500 bp. Samples were then amplified by PCR for 14 cycles using Phusion polymerase, NEXTflex Primer mix and the size-selected DNA as per the Bioo protocol, followed by an Agencourt bead clean up. Libraries were validated for integrity using the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). Clusters were formed and sequenced on the Illumina HiSeq2000 sequencers by the University of Wisconsin, Madison DNA Sequencing Facility in the University of Wisconsin, Madison Biotechnology Center. DNA clusters were generated using a cBot Single Read Cluster Generation kit (ver. 3) on an Illumina cBot (Illumina) according to the manufacturer's instructions to obtain an average of 1.5 ϫ 10 8 clusters for each flowcell lane. All 50mer sequencing runs were performed on an Illumina HiSeq2000 using the Illumina Sequencing kit (ver. 3). Fluores-cent images were analyzed using the CASAVA 1.8.2 (Illumina) to obtain FASTQ formatted sequence data. ChIP-seq samples were analyzed and prepared according to the ENCODE consortia best practices guidelines (36). Twelve barcoded libraries were run simultaneously per lane over 4 lanes. The raw FASTQ for each sample was concatenated from the 4 lanes prior to mapping to create a single sample. Each ChIP sample was repeated in biological replicate (Pearson Correlation Coefficient, PCC Ͼ 0.85) in the same manner (minimum of 2 replicates). Sequences were mapped to the mouse genome (mm9) using BOWTIE (best -m 1) to yield unique alignments (37). All sequencing tag statistics (total and mapped reads) are listed in supplemental Table S1.
Bioinformatic and Statistical Analyses-Samples were further processed by QuEST and HOMER. Peaks were accepted if they passed criteria for both methods. QuEST 2.4 was run using the recommend settings for transcription factor (TF) like binding with the following exceptions: kde_bandwith ϭ 30, region_ size ϭ 600, ChIP threshold ϭ 35, enrichment fold ϭ 3, rescue fold ϭ 3 (38). HOMER analysis was applied using the default settings for peak finding (39). A False Discovery Rate (FDR) cut off of 0.001 (0.1%) was used for all peaks. The tag density for each factor was normalized to 10 7 tags and displayed using the UCSC genome browser. Positive peak regions were interrogated further using Homer as well as the Genomatix Software Suite (40). Motif analysis (de novo and known) was performed using the HOMER software and Genomatix. Peak overlaps were processed with HOMER and Galaxy (41). Peak comparisons between replicates and between differentiation time points were processed with the EdgeR and DESeq statistical packages in R where p values were drawn from the tag densities under each peak (42)(43)(44). Samples were further sorted by p value and FDR. All data are deposited in the Gene Expression Omnibus (GSE41955).

RUNX2 Cistrome Analysis in POBs and in Vitro
Differentiated OBs-The MC3T3-E1 cells (POBs) used in this study were initially validated for osteoblast differentiation following a 15 day treatment with ascorbic acid and ␤-glycerophosphate (OBs) using Von Kossa and Alizarin Red staining procedures that verify appropriate differentiation, matrix secretion, nodule formation, and calcification indicative of mature osteoblasts. Cells were also examined for the presence of alkaline phosphatase and alkaline phosphatase enzymatic activity as well. As can be seen in Fig. 1, each of these assays confirms that this cell line undergoes osteoblast differentiation in response to these typical osteogenic-inducing agents. To explore the RUNX2 cistrome in these osteoblast-lineage cells, we then conducted a quantitative genome-wide comparison of RUNX2 binding sites in POBs and day 15 OBs. Accordingly, confluent POBs or differentiated OBs were subjected to replicate ChIP-seq analysis as described under "Experimental Procedures" using an antibody for RUNX2. The normalized list of genomic binding sites (derived by multiple methods, FDR Ͻ 0.001) for RUNX2 in both POB and OB cells were then used for all downstream bioinformatic analysis. As seen in Fig. 2A, 12,674 high confidence RUNX2 binding sites were observed throughout the POB genome. A more restricted 6,272 sites were observed in OBs, suggesting a striking reduction in the RUNX2 cistrome following differentiation. As 4,601 of these sites were found to overlap the RUNX2 cistrome in undifferentiated POBs, 8,058 and 1,641 sites appear to be unique to POBs and OBs, respectively. Interestingly, while Fig. 2a documents the loss of a large number of unique RUNX2 sites in POBs (8,058) and the appearance of a smaller number in OBs (1,641), absolute levels of RUNX2 binding based exclusively upon tag density changes (p Ͻ 0.05) at many of the observed sites were also variable as well (Fig. 2B). This latter assay does not take into account a statistical cut off alone, but rather the read densities as well. Thus, although the level of RUNX2 binding remained statistically unchanged at 7,875 sites, relative levels of RUNX2 were statistically higher (2-fold or greater, p Ͻ 0.05) at 5,380 sites in POBs than OBs and 1,045 sites in OBs than POBs. Thus, both absolute as well as relative changes in RUNX2 binding were observed by ChIP-seq analysis. The reduction in the RUNX2 cistrome is likely due to differentiation, as Western blot analyses, as documented in Fig.  2C suggest that the overall levels of RUNX2 protein expression following differentiation are maintained if not slightly elevated (changes not significant).
The RUNX2 cistromes in these two differentiation states also exhibited a number of common features as well. As seen in Fig.  2D, an analysis of their location across the genome, as correlated by GREAT (Genomic Regions Enrichment of Annota- tions Tool) association of genomic regions to nearby genes (45), revealed that while RUNX2 binding sites were indeed located near annotated gene promoters (ϳ18% within 5kb of TSS), they were preferentially found distal to gene promoters within both introns and intergenic regions. To further characterize sites of genomic RUNX2 binding in the two cell types, we next performed a de novo motif finding analysis for over represented binding sites present in these data sets. As can be seen in Fig. 2E, a classic RUNX2/AML/RUNX1 binding motif was found in these regions 61% of the time (compared with 10% in 50,000 GC matched background regions). Importantly, many of these RUNX2-bound sites contained multiple motifs for RUNX2, suggesting the possibility that the overall level of RUNX2 mass at many sites was influenced by the absolute number of RUNX2 binding motifs within the peak. Interestingly, potential motifs for other transcription factors including TEAD, IRF4, C/EBP isoforms, and AP1 were also observed (subset shown in Fig. 2E). Finally, Nearest Neighbor Analysis (NNA) (using GREAT) revealed a distance-dependent relationship between each site and the nearest annotated gene(s), which suggests that the average number of RUNX2 binding sites per associated gene is greater than one (data not shown). Critically, however, this metric is not definitive for the existence of a specific regulatory relationship between sites of RUNX2 binding and the genes located nearby. As seen in Fig. 2F, GO term analysis of genes expressed in POBs, OBs or both and associated with RUNX2 binding sites suggests that all sets of genes are associated with signal transduction components essential to osteoblast differentiation and transcription (including Sp7, Vdr, and Sox9 known to be involved in skeletal development) and gene networks involved in osteoblast function such as matrix organization, cell adhesion, and skeletal development/ossification. They include fibroblast growth factor 18 (Fgf18), alkaline phosphatase (Alpl), osteopontin (Spp1), collagen 1a1 (Col1a1) and the Runx2 gene itself, all of which are previously known targets of RUNX2 (46 -51). The ChIP-seq data tracks for several of these as well as novel genes are illustrated in Fig. 2G, where we document changes in relative RUNX2 binding using overlaid data track analysis as a function of the osteoblast differentiation state. Spp1, for example, represents a classic RUNX2 target and indeed, as previously characterized, RUNX2 was found near its promoter (52)(53)(54). Several additional RUNX2 binding sites were also present at the Spp1 locus, however, highlighting the concept that RUNX2 target genes are potentially regulated via multiple sites that are frequently located distal to target gene promoters. Importantly, the overlaid ChIP-seq profiles revealed a significant decrease in the level of RUNX2 as a function of differentiation at several sites. As documented in the data tracks in Fig. 2G, both Fgf18 and Col1a1 represent additional examples. In both cases, while the ChIP-seq analyses confirmed previous data suggesting functional RUNX2 binding sites near the promoters for Fgf18 (between Ϫ0.3 and Ϫ0.8 kb) (47) and Col1a1 (Ϫ1.376 kb) (48 -50), they also revealed the presence of a number of additional upstream sites as well. As with Spp1, several of these sites also showed a striking change in RUNX2 binding as a function of differentiation. ChIP-qPCR is shown for validation of ChIP-seq in Fig. 2H for Spp1 and Fgf18 regions. Finally, RUNX2 also has been observed to regulate many additional genes including Alpl (51) and Sp7 (55) that are known to be relevant to osteoblast structure and function.
Co-occupancy of C/EBP␤ at the RUNX2 Cistrome-RUNX2 activity is known to be regulated at various levels by a multitude of transcription factors (2). Among these is C/EBP␤, a DNAbinding protein that is highly active in osteoblast-lineage cells during differentiation (56). Indeed, motif finding analysis seen in Fig. 2E suggests that C/EBP␤ may play a role in RUNX2 activity. To explore this further, we conducted a similar ChIPseq analysis of the C/EBP␤ cistrome in both POBs and OBs and assessed the level of overlap between the two proteins on a genome-wide scale. The analysis in Fig. 3A revealed that like RUNX2, C/EBP␤ also bound broadly to the POB genome (15,622 sites) and that its cistrome was restricted following differentiation to OBs (10,847 sites). The C/EBP␤ protein levels were unchanged in differentiation states by Western blot analysis (data not shown). In addition, while an extensive number of sites were found common to both differentiation states, as also seen with RUNX2, each cell type retained a unique subset of C/EBP␤ binding sites as well. Similarly, GREAT analysis revealed that these sites were located 1) both promoter-proximal as well as promoter-distal (Fig. 3B, data for both POB and OB overlap in distances from TSS), 2) in a distant-dependent fashion near annotated genes (Fig. 3A) and 3) that the average number of binding sites per associated gene was also greater than one (data not shown). De novo motif finding analysis, as documented in Fig. 3C, also confirmed that the most over represented sequence found at 56 -67% of these sites reflected a classic C/EBP␤ binding motif; the abundance of RUNX2 and other transcription factor binding motifs was very high as well. Perhaps most importantly, as seen in Fig. 3D, a comparison of   Arrows link the C/EBP␤ binding site subsets (POB specific, OB specific, and POB/OB overlap) to annotated neighboring genes using GREAT. B, GREAT analysis of C/EBP␤ sites as a function of distance from TSS in kilobases (kb) for both POB and OB. C, de novo sequence over-representation analysis of C/EBP␤ peaks by HOMER in POBs and OBs. Abundance shown as percentage (black) compared with 50,000 GC content matched sequences (red). D, Venn diagram overlap of RUNX2 and C/EBP␤ peaks for both POB and OBs. E, ChIP-seq tag density (ChIP-Fragment Depth per bp per peak) for all RUNX2 peaks (left) and all C/EBP␤ peaks (right). RUNX2 peaks for POB (yellow) and OB (dark yellow), C/EBP␤ peaks for POB (green) and OB (dark green). F, the top 4 scoring Gene Ontology (GO) terms are presented for Molecular Function (MF), Biologic Process (BP), or Cellular Component (CC) for the regions in which RUNX2 and C/EBP␤ overlap. G, ChIP-seq tag density tracks for RUNX2 and C/EBP␤ (POB, yellow; OB, blue; overlap, green). Further details as in Fig. 2. the location of both RUNX2 and C/EBP␤ binding sites in each of the two osteoblast-lineage states revealed that while many sites of C/EBP␤ occupancy were independent of RUNX2, the two proteins co-occupied sites 31% of the time in POBs and 25% of the time in OBs. This relationship can also be seen via RUNX2 and C/EBP␤ ChIP-seq tag density at these regions; as seen in Fig. 3E, there was a very close association between the two proteins. Although none of these findings prove a direct relationship between C/EBP␤/RUNX2 co-occupied binding sites and the genes that are located most proximally, GO term analysis of both potential target gene sets, as listed in Fig. 3F, revealed genes associated with terms such as adhesion, growth factor binding, matrix organization, motility, transcription factors and phosphatase activity. These findings support the idea that RUNX2 binding sites reside within regulatory regions that are potentially modular in nature and that these regions may serve as enhancers for the genes to which they are linked. Fig.  3G documents the presence of C/EBP␤ at sites within several gene loci, including Spp1, Sp7, and Tln several of which also contain RUNX2.

RUNX2 and C/EBP␤ Cistromes Are Marked by Mono-and Tri-methylated H3K4 and Acetylated H3K9 and H4K5 Indicative of Enhancers-A wide variety of selective modifications at
both DNA and histones defines a complex epigenetic layer that provides both structural/architectural and functional elements to the genome (57,58). Importantly, while many of these modifications highlight essential features of genes on a genomewide scale, a subset are now known to define the regulatory or enhancer regions of genes. Given RUNX2's role as a primary determinant of osteoblastogenesis, we anticipated that sites of both RUNX2 and C/EBP␤ binding would also be enriched for enhancer histone marks and that these and additional marks might well define the regulatory epigenome of osteoblast-lineage cells. Thus, we used ChIP-seq analysis to map sites of enrichment of histone methylation reflecting transcription units (H3K36me3), gene promoters (H3K4me3) and enhancers (H3K4me1 and H3K4me3) and histone acetylation enrichment indicative of gene activity (H3K9ac and H4K5ac) in both POB and OB cells, and contrasted the relative levels of these marks during differentiation (58). CEAS (Cis-regulatory Element Annotation System (59)) analysis depicted in Fig. 4A, which assesses ChIP-seq tag density collectively across a normalized gene body, confirmed appropriate enrichment of these individual histone marks as has been previously described (58). As anticipated, while transcription factor binding was characterized by well-defined peaks, histone marks reflecting either overall gene activity or the presence of an enhancer were comprised of rather extended peaks that frequently spanned binding sites for multiple transcription factors (see data tracks). Fig.  4B documents the read densities for each histone modification centered at overlapping genomic RUNX2 and C/EBP␤ sites independent of the state of differentiation. Interestingly, while enhancer marks were broadly elevated across enhancers, the level of these marks were frequently depressed or absent specifically at sites of RUNX2 and/or C/EBP␤ occupancy (footprints) (60,61), indicating that RUNX2 binding (or a complex) reduces the level of histone modification, perhaps due to nucleosome displacement. Despite this, it is clear from these epigenetic data as documented in Fig. 4C that sites of binding for these two transcription factors are indeed located within the regulatory regions of genes. Importantly, distinct enhancer histone marks could also be observed at sites that did not contain RUNX2 and/or C/EBP␤, suggesting the presence of additional regulatory enhancers at these genes that were not defined by the presence of RUNX2 or C/EBP␤ and likely mediate gene regulation through additional transcription factors. Finally, we examined whether the level of these epigenetic marks was changed at the genes that were nominally associated with differential RUNX2 binding in POBs and OBs. No correlation was observed between the levels of histone enhancer (H3K4me1) or activity (H3K9ac and H4K5) marks and individual gene subsets defined via the RUNX2 binding sites that were either common to both POBs and OBs or unique to either cell differentiation state (data not shown). As DNA bound RUNX2 can exert both stimulatory and suppressive activities, this result is probably not surprising.
The Relationships between the RUNX2 Cistrome and Gene Expression-Studies suggest that RUNX2 is an early determinant of osteoblast lineage development and function and that the ability of RUNX2 to bind to DNA and to regulate gene expression directs these biological activities (2,4). We therefore explored the possibility that RUNX2 binding activity might change at subsets of genes whose expression levels changed as a result of differentiation, and that changes in histone modifications might be detectable as well. To accomplish this, we first conducted a DNA microarray analysis of RNA isolated from POBs and OBs and processed the data as described under "Experimental Procedures." All data passing a 95% confidence limit greater than 2-fold (up or down) were included in the analysis and are summarized as log2 values in supplemental Table S1. As documented in the heat map shown in Fig. 5A, although the expression levels of the vast majority of the genes evaluated by microarray were unchanged (Ͻ2-fold), 721 genes displayed either an increase (417) or decrease (304) (Ͼ2-fold) in expression between the two differentiation states. Many of these changes were analyzed via direct RT-qPCR and are documented in Fig. 6. Recognizable osteogenic genes that are modulated include Ank, Dmp1, and Spp1. Interestingly, despite the steady increase in Alp staining and enzymatic activity documented during differentiation (Fig. 1), changes in the expression of the Alpl gene at the RNA level was not observed (Fig. 6) despite repeated analyses. The Bglap2 gene was similarly not up-regulated. The reason for the absence of Alpl and Bglap2 as well as Runx2 induction following differentiation is not clear, although it is possible that changes in the expression of these genes occurs prior to day 15. Despite this, however, the changes in gene expression that are observed together with the cellular analyses documented in Fig. 1 are supportive of the differentiation and maturation of these MC3T3-E1 cells from earlier precursors to matrix-secreting and mineralizing cultures over the time frame of these genomic analyses. We then contrasted the relative distribution of all RUNX2 peak-associated genes (obtained by GREAT analysis) identified individually in POBs, OBs, or both with genes whose expression levels were changed during differentiation. As can be seen in Fig. 5B, this relative distribution across the three categories was similar. In fact, 60 of the 108 genes that were up-regulated and 32 of the 47 genes that were down-regulated in OBs were also found within the other two categories. Indeed, as documented in Fig. 5C, only when a comparison was made between GREAT-associated genes whose overall basal level of expression was in the lowest 20% (genes that are generally "off") versus the highest 20% (genes that are "on" and highly expressed) was this distribution of genes significantly elevated toward the highly expressed cohort. Given these results, it was not surprising as seen in Fig.  5D that the GO terms for the differentially expressed genes enriched in each category were similar. This observation suggests that regardless of whether the expression levels of the  ). B, ChIP-seq tag density focused on overlapping RUNX2 and C/EBP␤ peaks for the POB and interaction with histone markers. C, ChIP-seq tag density profiles for RUNX2, C/EBP␤, H4K5Ac, H3K9Ac, H3K4me1, H3K4me3, and H3K36me3 (POB, yellow; OB, blue; overlap, green) at the Igfbp5 and Tln1 gene loci. Further details as in Fig. 2. genes are changed during differentiation or not, both categories represent potential targets of RUNX2 action.
As discussed earlier, gross histone modifications were not observed on a genome-wide scale when attempting to correlate RUNX2 binding. If these gene groups were further subdivided into those that are up and down-regulated Ͼ2-fold, however, measurable correlative changes can be observed as assessed by CEAS analysis in Fig. 5E. Accordingly, the level of H4K5ac is decreased at genes that are down-regulated and increased at genes that are up-regulated during differentiation; similar cor-  , moderated t test). Data are displayed as log2 intensities; scale is indicated at bottom where red represents high and blue low expression. B, RUNX2 GREAT associated genes (6,856, 4,367, and 2,063) are cross referenced to up-and down-regulated genes (60 up-regulated and 32 down-regulated genes are in common with all three peak states). C, GREAT-associated genes are compared with the top 20% (4,487 genes) and bottom 20% (3,483 genes) quintiles for gene expression (overlapping accession numbers were removed) in POB and OBs. D, gene ontology of each category from B shows similar enrichment for these categories. E, CEAS comparison of average gene profiles for H4K5Ac in POB (solid lines) and OBs (dashed lines). Up-regulated genes are shown in red; down-regulated genes are shown in blue; all genes are shown in black. F, ChIP-seq tag density profiles for RUNX2, C/EBP␤, H4K5Ac, H3K9Ac, H3K4me1, H3K4me3, and H3K36me3 (POB, yellow; OB, blue; overlap, green). Further details as in Fig. 2. relations were also evident with other histone marks. Although these changes were modest across this cohort of regulated genes, a number of target genes that are highly regulated show a striking increase or decrease in the overall level of H3K4me1, H3K9ac and/or H4K5ac upon differentiation when correlated with RUNX2 binding (Fig. 4C and Fig. 5F); genes whose expression remains unchanged did not show such epigenetic changes.

DISCUSSION
RUNX2 is a lineage determination factor that is expressed early in transitioning MSCs and essential for commitment to the osteoblast lineage (2,4,5,12,14,62). Its role is clearly multifactorial in that while RUNX2 orchestrates the initial formation of osteoblast precursors, it is silenced during chondrogenesis and then re-expressed together with RUNX3 to drive hypertrophic chondrocyte differentiation (4). Other members of the RUNX family are similarly involved in cell fate decisions through the regulation of cellular growth and differentiation, and the aberrant or miss-expression of each of the members is associated with unique cancer phenotypes (25,(63)(64)(65). Interestingly, the role of RUNX2 is primarily early in the osteoblast lineage. Accordingly, while genetic removal of RUNX2 results in a total absence of osteoblasts in mice in vivo, selective deletion of RUNX2 in mature osteoblasts does not appear to be deleterious for either continued differentiation of the lineage or for maintenance of the bone mineralization program (66). Clearly, other proteins are involved in these processes as well including OSX, which is both induced by RUNX2 and is capable of directing downstream maturation and function of the osteoblast lineage even in the absence of RUNX2 (21). Additional participating factors include cell autonomous regulatory agents involved genetically in the osteoblast program, those that interact to control RUNX2 activity directly and those whose expression levels or activities are regulated by external, ligand-regulated signal transduction pathways (29). Regardless of these details, however, RUNX2 functions as a lineage-determining factor together with synergistic partners to modulate the expression of both broad as well as cell-specific gene networks by contributing to the formation, maintenance and activity of the enhancers that regulate the expression of the individual gene components of these networks. The importance of RUNX2 is highlighted by the additional fact that this protein also plays a significant epigenetic role in osteoblastogenesis by maintaining its presence at key genes throughout mitosis (26). Many of these features are consistent with that of other lineage determining factors, each of which acts in concert with additional lineage-specific factors to orchestrate individual programs of differentiation.
Central to understanding RUNX2 actions further is a quantitative and qualitative genome-wide assessment of the cis-regulatory sites to which RUNX2 binds both in early osteoblast precursors as well as more mature osteoblasts and identification of the properties of these binding sites that are essential to the regulation of the genes associated with the osteoblast phenotype. As determined by ChIP-seq analysis, we showed that the RUNX2 cistrome was comprised of 12,674 high confidence binding sites that were widely dispersed across the POB genome; this cistrome contracted to 6,272 sites upon differentiation to the more mature cell phenotype. This decrease is likely consistent with the idea that RUNX2 represents an early phenotype regulator whose influence becomes increasingly restricted as differentiation proceeds. As total RUNX2 protein levels were largely unchanged following differentiation, this reduction seems likely to result from either a change in RUNX2 activity or a change in the differentiated cell's genome. Despite the reduction in the RUNX2 cistrome following differentiation, the features of the two were similar: cis-regulatory sites are located both promoter-proximal as well as promoter-distal, frequently positioned many kilobases from potential target genes and contain one or more consensus RUNX2 response elements (RREs) that are known to represent sites of direct interaction between RUNX2 and DNA. In addition, both the genome-wide analyses and visual inspection of numerous data tracks at known target genes indicated that more than one RUNX2 regulatory region was often present at potential target gene loci. The fact that multiple RUNX2 enhancers frequently regulate specific target genes leads to additional complexity in defining the role for RUNX2 at these particular genes. Importantly, the RUNX2 cistromes also overlapped extensively with those of C/EBP␤, an osteoblast lineage-active factor that is known to interact directly with RUNX2 to influence RUNX2 activity (56,67,68). C/EBP␤ plays distinct roles in the differentiation of osteoblasts, chondrocytes and adipocytes and, like RUNX2, its cistrome also undergoes restriction upon osteoblast differentiation. Thus, C/EBP␤ represents a prototype co-regulatory DNA binding factor for RUNX2. Given the abundance of proteins that interact with RUNX2 coupled with the frequent presence through our de novo analysis of additional DNA motifs found near sites of RUNX2 binding, it is likely that numerous other proteins form complexes at these sites as well. Finally, sites to which either RUNX2, C/EBP␤, or both bind were also marked by post translational histone modifications that represent either bone fide signatures of enhancers (H3K4me1 and 3) or indicators of gene activity (H4K5ac and H3K9ac) (58). Importantly, the levels of these marks frequently changed at many loci upon differentiation indicating possible changes in gene activity. These data provide significant insight into the sites of action of RUNX2, the significance of C/EBP␤ as a coregulatory DNA binding factor for RUNX2 activity, and of the identity of target genes they potentially regulate. Our study of the RUNX2 cistrome in osteoblasts appears unique as similar analyses have not yet been reported. Two recent studies of RUNX2 binding activity have been conducted, however, the first reported a ChIP-seq analysis of the RUNX2 cistrome in a RUNX2-null prostate cell line engineered to express RUNX2 in a doxycycline-sensitive manner (33) and the second reported a ChIP-chip analysis of RUNX2 binding at a genome-wide collection of gene promoter regions in the osteosarcoma cell line MG63 (32). The latter identified a collection of 2,265 RUNX2 binding sites that were correlated with adjacent genes, several of which were investigated further as contributing to cancer cell migration and adhesion. As RUNX2 sites located distal to gene promoters represent the norm rather than the exception, as illuminated unequivocally in the present study, the general lack of correlation between RUNX2 binding and target gene regulation evident in this earlier study is likely due to the lack of detection of distal RUNX2 sites. Interestingly, siRNA mediated down-regulation of RUNX2 was not entirely effective in depleting RUNX2 binding at many target genes (32), suggesting that this approach could be problematic. Approximately 1,603 RUNX2 binding sites were noted across the prostate cell line genome upon induction of RUNX2 (33). A parental prostate cell line control expressing endogenous levels of RUNX2 was not examined, however. Thus, it is unclear whether the number of sites measured for RUNX2 represented an appropriate quantitation of the RUNX2 cistrome in this cell line or simply a low estimation due to the technical nature of the model. As the role of RUNX2 in gene expression is complex, it is possible that while the controlled induction of RUNX2 appears to represent a novel approach for defining bone fide target genes in prostate cells, the acute appearance of RUNX2 in this system may underestimate this protein's role at many target genes. In view of the number of sites identified in the present study, the latter seems the more likely possibility. In both studies, however, direct linkage was made with several genes that were indeed regulated by RUNX2, including several involved in motility, migration, adhesion and metastasis.
While quantitation and characterization of the RUNX2 cistrome in POBs and OBs is informative, perhaps of greater importance is an assessment of how these sites and the changes that are observed during differentiation correlate with changes in gene expression patterns. An analysis of gene expression in POBs and OBs revealed a large overlap between the two transcriptomes and that expression levels of most of the genes were changed Ͻ2-fold as a consequence of differentiation. However, 721 changes in gene expression Ͼ2-fold (up or down) were noted. Regardless of regulation, many of the genes in both of these categories were identified as contributing to the osteoblast lineage phenotype and therefore represent potential targets of RUNX2 action. Demonstrating this fact directly, however, is more difficult. Accordingly, while GREAT analysis points to some 11,232 genes associated with the RUNX2 cistrome in POBs and 6,430 genes associated in OBs (and 4,601 RUNX2 sites and 4,367 associated genes there were found to be in common in both), the overall net distribution of RUNX2 sites at these gene subsets was not different from that seen for the RUNX2 sites associated with the 721 regulated genes. In fact, only when a comparison was made between genes with the lowest (20%) and highest (20%) levels of overall expression was an increase noted in the abundance of RUNX2 binding sites at genes that were expressed most highly. Are these RUNX2 binding sites involved in gene regulation in general? We conclude that many are indeed linked to the regulation of adjacent neighboring genes based upon the fact that the majority of the sites located within regions are heavily marked by epigenetic enhancer signatures and that levels of these histone marks are increased at genes that are up-regulated during differentiation. This conclusion requires qualifications, however. First, the identity of the target gene for an individual enhancer cannot be confirmed through proximity alone. Thus, as noted in this RUNX2 study as well as those for other transcriptional regula-tors, enhancers are frequently located many kilobases from target genes (69 -72) and in several cases have been found to regulate genes on other chromosomes (73,74). A second qualification is that while DNA binding is likely a pre-requisite for RUNX2 activity, it alone does not indicate the state of transcription activity for RUNX2 (2,29). Indeed, RUNX2 activity is regulated by its abundance, through post-translational modifications that include both phosphorylation and acetylation and affect both concentration and activity, and through its ability either to interact synergistically with adjacent DNA binding factors or to recruit co-regulatory factors. A final qualification is a temporal one; while many sites of RUNX2 binding may be fully operational, it is possible that other sites reflect historical interactions (memory) or interactions that predict future activity. All of these qualifications point to the importance of these data sets as a future resource for further mechanistic studies aimed at individual genes and gene networks.
Regardless of the complexities associated with linking RUNX2 binding sites to genes on a genome-wide scale, we note the presence of RUNX2 near many genes that are either known or likely targets of RUNX2 action by virtue of their contribution to the osteoblast phenotype. With respect to known targets of RUNX2 action, which include Spp1 (52)(53)(54), Vdr (72), Rankl (75,76), and Runx2 itself (46), it is interesting that although our ChIP-seq analyses generally confirmed sites known from these previous studies to bind RUNX2, they frequently revealed additional sites as well providing novel insight into the mechanism through which RUNX2 regulates these target genes. This increase in complexity contributes in no small way to the difficulty of assessing the nature of RUNX2 regulation. For example, while RUNX2 binding activity may change in concert with gene expression at some regulated genes, RUNX2 binding at other regulated genes may be selectively altered at only one of several adjacent sites. In addition, sites may be differentially occupied as well depending upon the state of maturity. Detailed dissection of RUNX2 action at these gene loci will be required to sort out the potentially unique role(s) of each regulatory region in RUNX2 action. Importantly, in addition to confirmation of known gene targets, genome-wide analysis of the transcriptome herein has also revealed the presence of RUNX2 at a plethora of new gene targets as well. Clues as to their role in osteoblast biology can be found in the GO term categories that emerge for the gene expression analyses. Thus, as stated earlier, the greatest value of these data sets may be as a resource for continued investigation of individual genes.
In summary, we have defined the RUNX2 cistromes in both early and late osteoblasts in vitro, and characterized important features of these cistromes relevant to the activity of the lineage-determining factor RUNX2. We have also correlated these cistromes to gene expression patterns identified in early and late osteoblasts and in particular to those genes whose expression levels undergo change during the differentiation process. Our results add to an understanding of the mechanism through which RUNX2 plays its regulatory role in the development and maintenance of the osteoblast phenotype.
Addendum-While this manuscript was in revision, a report by Wu and colleagues (77) was published that described the results of a similar analysis of the RUNX2 cistrome in MC3T3-E1 clone 4, a more rapidly differentiating version of the MC3T3-E1 cell line used in our studies. The results presented in this new report provide evidence for a much more extensive RUNX2 cistrome than that described herein, although the overall genomic properties of these sites are similar to those we have identified. It is difficult to determine why the RUNX2 cistrome in this report by Wu et al. is considerably larger than that we have observed, although it is not due to differences in the bioinformatic analyses that were performed in the two laboratories. Nevertheless, we note that our data were generated from nearly identical ChIP-seq replicates of the RUNX2 cistrome across the MC3T3-E1 genome at each state of differentiation. The replicate nature of the data sets generated by Wu et al. is unclear as only one set appears to be available through GEO. More importantly, however, our analyses also include the observation that a high percentage of the RUNX2 sites that we identified also contain co-localized C/EBP␤ and epigenetic histone signature marks of active enhancers (H3K4me1, H4K5ac, and H3K9ac). These features of the RUNX2 cistrome we have identified are generally absent at many of the additional RUNX2 sites identified in the report by Wu et al. and all are available at GEO. While future studies will be necessary to resolve the minor differences highlighted in these two reports, the general identification and characterization of the RUNX2 cistrome in differentiating MC3T3-E1 cells provides new insight into the DNA binding properties of this key osteoblast-lineage determining factor, and will likely enable more detailed studies on its function.