|Home | About | Journals | Submit | Contact Us | Français|
Identifying gene regulatory elements and their target genes in human cells remains a significant challenge. Despite increasing evidence of physical interactions between distant regulatory elements and gene promoters in mammalian cells, many studies consider only promoter-proximal regulatory regions. We identify putative cis-regulatory modules (CRMs) in human skeletal muscle differentiation by combining myogenic TF binding data before and after differentiation with histone modification data in myoblasts. CRMs that are distant (>20 kb) from muscle gene promoters are common and are more likely than proximal promoter regions to show differentiation-specific changes in myogenic TF binding. We find that two of these distant CRMs, known to activate transcription in differentiating myoblasts, interact physically with gene promoters (PDLIM3 and ACTA1) during differentiation. Our results highlight the importance of considering distal CRMs in investigations of mammalian gene regulation and support the hypothesis that distant CRM-promoter looping contacts are a general mechanism of gene regulation.
The identification of genomic sequence regions that regulate genes in a condition-specific manner is essential to understanding how the same genome sequence can give rise to the diversity of cell types and functions observed in an organism. In an organism with a small genome, such as yeast, the majority of gene regulation can be explained by transcription factor (TF) binding and chromatin modifications within approximately 600 bp to 1 kb of DNA sequence upstream of the regulated gene [1, 2]. In metazoans, numerous prior studies in a range of organisms from sea urchin to mammals have identified cis regulatory modules (CRMs), consisting of clusters of TF binding sites, located next to (or within the introns of) the genes whose expression they regulate . However, as compared to the yeast genome, metazoan genomes, in particular mammalian genomes, have a much higher proportion of noncoding sequence, and recent research has highlighted the importance of more distant CRMs in gene regulation within these genomes [4–6]. Much detailed work on distantly located transcriptional enhancers in Drosophila has shown the importance of such distant regulatory elements and the ways that they can be directed to their target genes by insulator boundaries and promoter targeting sequences . In mammalian systems, the rules governing enhancer-promoter links are less evident, but certain examples of distant regulatory elements have been studied in detail. For example, the locus control region (LCR) of the murine β-globin locus forms a GATA-1-dependent looping interaction with actively transcribed globin gene promoters located approximately 40–60 kb away in erythroid cells [8, 9]. Such results suggest that a complete understanding of gene regulation will require searching for CRMs distant from target genes and further studies of how these CRMs are directed to and regulate their target genes.
The differentiation of human skeletal myoblasts into mature muscle fibers requires the coordinated regulation of many genes and is essential for the development and maintenance of proper muscle function. This differentiation process can be easily induced and monitored in cell culture, making it a tractable model system for investigating the regulatory mechanisms underlying dynamic, tissue-specific gene expression in human. Prior studies have measured changes in gene expression and TF binding during differentiation in primary human skeletal muscle cells or in the similar C2C12 mouse skeletal muscle cell line [10–14]. These studies have shown that a master regulatory TF, MyoD, initiates the cascade of gene regulation that leads to fusion of undifferentiated myoblast cells into multinucleated, elongated myotubes with developed contractile elements [10, 15]. Other TFs in the myogenic regulatory factor (MRF) basic helix-loop-helix (bHLH) family -- myogenin (MyoG), Myf5, and MRF4 -- some of which are activated directly by MyoD, work together with MyoD and other TFs such as Serum Response Factor (SRF) and the Myocyte Enhancer Factor 2 (Mef2) family to regulate the expression of genes involved in the continuation and completion of this differentiation process [10, 15, 16]. Using the sets of differentially expressed genes, conservation of sequence across species, and the DNA binding site motifs of myogenic TFs, CRMs responsible for coordinating changes in expression during myogenic differentiation have been predicted computationally [12, 17, 18]. However, unlike the β-globin locus in which some of the classic distant regulatory interactions have been characterized previously, little is known about the role of distant regulatory modules in human skeletal muscle differentiation.
Previous efforts to characterize the transcriptional regulatory network in skeletal muscle differentiation have measured binding of myogenic TFs using chromatin immunoprecipitation followed by microarray hybridization or sequencing (ChIP-chip or ChIP-Seq, respectively) in murine C2C12 cell lines [10, 11, 14] and have tested the regulatory function of some predicted TF-bound elements with reporter assays . While the results of these studies have revealed previously unknown regulatory connections between TFs and differentiation processes, they tended to focus on TF binding near promoters (approximately 1–4 kb upstream of transcriptional start sites (TSSs))[10, 11, 18]. In contrast, methods to predict CRMs involved in myogenic differentiation have identified candidate regulatory modules by searching sequences up to 10 kb  or 50 kb  away from TSSs. A handful of the more distantly located predicted CRMs (approximately 20–30 kb upstream or downstream of TSSs) were found to activate expression in reporter assays specifically during myogenic differentiation and to be bound by myogenic TFs . These results suggest that gene regulatory myogenic TF binding occurs at locations distant from the proximal promoter regions that have been the focus of most prior studies.
In this study, we identified genomic regions with a potential regulatory role in myogenic differentiation by combining experimental measurements of myogenic TF binding before and after muscle differentiation with previously published locations of enhancer-associated chromatin modifications in myoblasts. We found such putative CRMs located both proximal to and distant from muscle genes, and observed that distant CRMs showed changes in TF binding during differentiation more often than did CRMs located in proximal promoter regions. Investigating the potential regulatory roles of such distant candidate CRMs, we found that a few are located near known microRNAs, which we propose might be involved in myogenic differentiation. To determine whether distant CRMs with no microRNA or protein-coding gene promoters nearby might form long-range looping interactions with muscle gene promoters, we examined two distant CRMs previously found to drive gene expression specifically during myogenic differentiation as case examples for further inspection. We found that these CRMs form differentiation-specific physical interactions with their closest target genes. Our results provide further support that a complete understanding of transcriptional regulation in mammals will require consideration of CRMs located distant from their target genes in the genome sequence.
To search for candidate cis regulatory regions utilized in human skeletal muscle differentiation, we characterized the binding of the known myogenic TFs MyoG, MyoD, and SRF to regions surrounding muscle genes before and after differentiation. ChIP-chip analysis of TF binding in 100 kb regions surrounding each of 104 muscle genes resulted in 1,781 significant (Bonferroni-corrected; α = 0.05), biologically replicated 1-kb regions of DNA bound by either MyoG or SRF immediately before (0 hours) differentiation and 3,702 1-kb regions bound by either MyoG or SRF 48 hours after differentiation.
Prior studies have shown that considering clusters of TF binding sites and certain chromatin modifications can aid in the identification of functional regulatory regions, narrowing the set of candidate CRMs from the large set of regions bound by single factors to a smaller, higher-confidence set [19–21]. Therefore, we investigated which combinations of TF binding and chromatin features might be predictive of CRMs involved in myogenic differentiation. We considered the overlap of MyoG and SRF binding from our ChIP-chip data and the changes we measured in this binding from 0 h to 48 h of differentiation, as well as supplementary single-replicate ChIP-chip data that we collected for MyoD binding at 0 h and 48 h and for the enhancer associated protein p300  binding at 48 h after differentiation. In addition to our own experimental data, we considered previously published data on histone modifications in human myoblasts (corresponding approximately to the 0 h of differentiation timepoint), focusing on marks that have been shown to be associated with enhancers both near and far from gene TSSs (H3K4me1 and H3K27Ac) . As an initial guide for future predictions of putative CRMs, combinations of the above features were evaluated for their ability to distinguish between 9 positive control CRMs and 9 negative control genomic regions. These control regions were identified previously as either active or inactive, respectively, in regulating gene expression during myogenic differentiation, as evidenced by their being occupied by MyoD, MyoG, and SRF by ChIP-qPCR and by their ability to activate reporter gene expression during differentiation  (Table 1).
We observed MyoG binding at 48 h after differentiation for all 9 positive controls and only 2 of the negative controls, making this a strong indicator of CRMs active during myogenic differentiation. The binding of SRF, MyoD, and p300 at 48 h post-differentiation co-localized with MyoG binding at subsets of the known CRMs, but including the ChIP data for these factors did not help to further distinguish the positive control CRMs from the negative control regions (Table 1). The binding of MyoD or SRF alone at 48 h post-differentiation was less specific than MyoG binding for distinguishing these sets of regions: we observed SRF binding at 7 positive and 7 negative controls and MyoD binding at 7 positive and 3 negative controls.
Combining TF binding information with data on enhancer-associated histone marks H3K4me1 and H3K27Ac in human myoblasts  further helped to distinguish the positive and negative controls. We found that the CRMs bound by MyoG at 48 h post differentiation and that had either H3K4me1 or H3K27Ac modifications in myoblasts (i.e., immediately before differentiation) includes all of the positive controls and none of the negative controls. MyoG binding at regions with enhancer-associated chromatin modifications not only identified positive control CRMs located in proximal promoter regions, but also distinguished previously validated distant (greater than 20 kb away from TSS) CRMs from negative control genomic regions (Figure 1A and B). Since this combination of enhancer marks and MyoG binding best discriminates the positive control CRMs from the negative controls, in subsequent analyses we considered all regions with these features as a set of putative muscle differentiation CRMs.
The promoter-associated histone mark H3K4me3  also provided good distinction between the known CRMs and the negative controls, but the levels of this mark were only marginal for validated muscle differentiation enhancers more distant from promoters (Figure 1A), in agreement with previous observations that H3K4me3 is generally not observed at enhancers located more than 2 kb away from TSS . Therefore, we did not include this feature in subsequent analyses of CRMs since our analyses broadly encompassed both proximal and distant regulatory regions.
MyoG binding during differentiation and enhancer-associated histone marks in myoblasts not only distinguished positive control CRMs from negative control genomic regions, but also overlapped significantly with computational predictions of CRMs identified by searching for clusters of conserved MRF DNA binding site motif matches with the algorithm PhylCRM . This MRF motif (the motif recognized by MyoG and MyoD) was more enriched among a set of genes upregulated during myogenic differentiation than any other combination of MRF, Mef2, SRF, or Tead motifs in a previous study . Of the 1,975 regions enriched for MyoG binding and enhancer histone marks, 1,128 contain a CRM predicted by PhylCRM (Supplementary Figure S1). This significant overlap (p<1×10−14 ; Hypergeometric test) between computationally predicted CRMs and experimentally identified putative CRMs further demonstrates that PhylCRM can be used successfully to predict likely CRMs and that the observed MyoG binding events are likely mediated by direct interactions with MyoG binding sites within the predicted CRMs. Neither correlation is perfect, however: 43% of regions with experimentally observed MyoG binding and chromatin modifications were not PhylCRM-predicted CRMs, and 67% of PhylCRM-predicted CRMs do not show these in vivo indicators of regulatory function. TF binding at regions not predicted by PhylCRM could be indirect association through other TFs, rather than by the myogenic factor binding motif. Likewise, PhylCRM predicted regions that are not bound could be CRMs in other conditions. Given the limitations of both experimental TF binding measurements and computational predictions, it is important to combine such different sources of evidence to identify putative CRMs.
As a larger set of putative elements than the few positive and negative controls considered earlier, PhylCRM predictions allow to further evaluate the patterns of TF binding at putative regulatory regions. Interestingly, the binding of some TFs, such as SRF, that were not as specific as MyoG for positive control CRMs, do show significant overlap with PhylCRM-predicted CRMs. As measured by the sum of sensitivity and specificity, SRF binding at 48 h is nearly as effective as MyoG binding at 48 h in discriminating PhylCRM-predicted CRMs from genomic regions that did not contain PhylCRM hits (Supplementary Table S1). This suggests that SRF binding is an informative feature of many muscle differentiation CRMs, even though it is not highly specific in our small set of validated positive and negative controls. The degree of overlap between computational predictions and different combinations of TF binding and/or chromatin marks across different timepoints is shown in Supplementary Table S1.
Many of the genes whose predicted CRMs were included on our custom-designed ChIP-chip microarray were chosen because they were upregulated over a timecourse of human skeletal muscle differentiation . To evaluate whether changes in TF binding might be a hallmark of cis regulatory regions of muscle genes and might effect changes in their gene expression levels, we compared TF binding at 0 and +48 hours relative to the induction of myoblast differentiation by serum removal. Changes in TF occupancy occurred at many genomic regions during differentiation; for example, 1,033 of the 3,007 regions bound by MyoG at 48 h were not bound by any of the examined myogenic TFs (MyoG, MyoD, or SRF) before differentiation (0 h). However, in most cases including information about these changes in TF binding did not help to distinguish the positive control CRMs from the negative control regions since the patterns of TF binding changes during differentiation vary across the known regulatory regions. For example, the ACTA1 promoter is bound by MyoG and SRF at both 0 and 48 h while the COX6A2 promoter is unbound before differentiation and then bound by MyoG, MyoD, and SRF at 48 h. The validated distant enhancer region between the PDLIM3 and SORBS2 genes is bound by MyoD and SRF at 0 h and then by an expanded set of factors (MyoG, MyoD, and SRF) with stronger peak intensities at 48 h (Figure 1A). Only in the case of SRF binding did the ChIP timecourse data help to distinguish the positive control CRMs from negative control regions: of the 7 known CRMs bound by SRF at 48 h, 6 were also bound at 0 h, while in contrast only 1 of the 7 negative control regions bound at 48 h was also bound at 0 h.
In order to determine the value of any of the measured ChIP datasets in predicting the locations of TF binding events at 48 h during differentiation, we constructed a dependency network using the WinMine toolkit . A dependency network, like a Bayesian network, finds predictive relationships between variables using conditional probability distributions . The derived network indicates which variables (in this case, TF binding events at a given genomic region) are well predicted by other variables. The dependency network for our data suggests that MyoD binding at 0 h is more important for predicting the binding patterns of MyoG, SRF, MyoD, and p300 at 48 h after induction of differentiation than are either MyoG or SRF binding at 0 h or H3K4me1 marks in myoblasts (Supplementary Figure S2). These findings are consistent with prior literature on the pioneering role of MyoD in establishing cis regulatory sites in myogenic differentiation .
Most regions previously known to have a role in myogenic differentiation (including most of the positive controls considered in the previous analyses) are located within proximal promoter regions of protein-coding genes (within 1–2 kb of TSS), but some of the putative regulatory elements we identified are more distant from any protein coding genes (> 20 kb away from TSS; Figure 1A). To further investigate the role of distant genomic regions in the regulation of myogenic differentiation, we examined the distribution of distances between genes and putative regulatory elements enriched for myogenic TF binding and enhancer-associated histone modifications. There is an enrichment of these putative regulatory regions close to gene TSS (within 1 kb) as compared with the locations of all possible 1-kb regions on our custom ChIP-chip array (p<1×10−14, Hypergeometric test; Figure 2). Despite this skew toward proximal promoter regions, about 80% of regions bound by MyoG at 48 h with H3K4me1 or H3K27Ac modifications fall outside of traditional promoters (1 kb from TSS) and, strikingly, 25% of such putative CRMs are more than 20 kb away from any TSS. Interestingly, genomic regions occupied by different TFs and with different binding dynamics during differentiation differ in their distances from gene TSS. Genomic regions that either lose or gain MyoG binding from 0 h to 48 h tend to be located further from TSSs than those that have constant MyoG binding at 0 and 48 h (p = 1.39×10−4 for constant binding versus gain of binding at 48 h; Kolmogorov-Smirnov (KS) test; Figure 2A). Increased distance from gene promoters is also observed for regions that change in both MyoG and SRF binding between 0 and 48 h of differentiation as compared with those that are bound by both MyoG and SRF at both timepoints. In addition, regions bound by all 3 of these TFs (MyoG, MyoD, SRF) tend to lie closer to genes than regions with only one or two TFs bound (p < 2.2×10−16 for MyoG 48 h vs. MyoG, SRF, and MyoD at 48 h; KS test; Figure 2B).
For 69% of the regions with enriched MyoG binding at 48 h and H3K4me1 or H3K27Ac marks in myoblasts, the closest gene is not known to be involved in myogenic differentiation and is not differentially expressed during this process. In fact, 90% of such putative CRMs are located outside of 1 kb muscle gene promoter regions and 56% are more than 20 kb away from a muscle gene TSS (Figure 2C). Furthermore, for 30% of these MyoG-bound, H3K4me1 or H3K27Ac modified regions, associating with the closest muscle gene along the linear chromosome requires skipping over intervening non-muscle genes.
We specifically included 100 genomic regions on our custom ChIP-chip array design that contained 287 high-scoring PhylCRM predictions of myogenic CRMs with clusters of evolutionarily conserved myogenic TF binding site matches, but were located far (>50 kb) from any annotated muscle gene TSS. These distant predicted CRMs exhibited MyoG binding and enhancer associated histone marks at a somewhat lower rate (26%) than PhylCRM–predicted CRMs <50 kb away from muscle gene TSS (35%). However, permutation testing showed that this level of overlap was greater than expected (p = 0.02) for randomly selected genomic regions represented on the array and much greater than expected (p < 1×10−4) for genomic regions with the same median distance to the nearest TSS of any protein-coding gene (21.8 kb). The significant overlap between MyoG binding events, H3K4me1 or H3K27Ac histone marks, and the PhylCRM computational CRM predictions indicates that these putative distant CRMs may indeed be biologically functional CRMs.
To explore the possible functions of putative myogenic CRMs that are distant from muscle protein-coding genes, we searched for microRNA genes annotated in miRBase  that are located adjacent to these CRMs. We identified 14 different putative CRMs that are located within 25 kb of 10 different microRNA genes and that are closer to these microRNAs than to any known muscle protein-coding genes (Table 2). In 5 cases, the putative CRMs are closer to one of these microRNAs than to any protein-coding gene (highlighted in green in Table 2). Some of these microRNAs that may be regulated by the myogenic TFs bound at these putative CRMs already have an annotated function in muscle differentiation (highlighted in yellow in the table). For example, miR-133 and miR-1 influence skeletal muscle proliferation and differentiation, with miR-133 specifically downregulating the myogenic TF SRF, and miR-1 directly downregulating HDAC-4 and thus indirectly modulating Mef2 gene expression . Additionally, miR-499 is differentially expressed during murine skeletal muscle differentiation  and is known to be involved in cardiac myocyte differentiation  and other processes in skeletal muscle . The 7 remaining microRNAs (miR-9-2, miR-550a-2, miR-550b-2, miR-675, miR-3925, miR-4298, and miR-4322) we identified based on their proximity to putative myogenic CRMs do not yet have known regulatory functions in muscle differentiation and would be interesting to investigate further for their potential myogenic regulatory roles in future studies.
We used Chromosome Conformation Capture (3C) [32, 33] experiments to investigate whether the putative myogenic CRMs that we found located distant from muscle genes form physical interactions with their adjacent, putative target genes during myogenic differentiation. We selected two CRMs that were bound by MyoG at 48 h after induction of differentiation, exhibited H3K4me1 and H3K27Ac enhancer-associated histone marks, were found to drive expression of a reporter gene during differentiation , and were located at least 20 kb away from the TSS of the closest known protein-coding gene. The “ACTA1 CRM” (Figure 3A) is located 23 kb downstream of the TSS of ACTA1, a gene which encodes the skeletal muscle alpha actin protein, an essential part of the contractile element in muscle. The “PDLIM3/SORBS2 CRM” (Figure 3B) is located 32 kb upstream of PDLIM3, which encodes a cytoskeletal organizing protein localized to the Z line of the sarcomere and which may also regulate the myogenic differentiation transcriptional network by affecting the nuclear localization of SRF . This CRM is also located 244 kb downstream of SORBS2, another gene differentially expressed during myogenic differentiation.
To assay the physical interactions between these two CRMs and their adjacent genes, we performed 3C experiments on human muscle cells 48 hours before, immediately (“0 hours”) before, and 48 hours after induction of differentiation by serum removal. We identified “differentiation-specific interactions” by searching for any peaks of interaction in the 0 or 48 h post-differentiation cells, outside of the neighboring +/− 15 kb that may interact highly by random collisions with the fixed bait, that were higher than the 48 h pre-differentiation interaction level (see Methods). The results of these experiments show evidence for a differentiation-specific interaction between the “ACTA1 CRM” and the proximal promoter region spanning from 1.6 kb upstream to 1.3 kb downstream of the ACTA1 TSS (Figure 3A). We confirmed this looping interaction in cells at 48 h post-differentiation using two different restriction enzymes, BglII and AflII (Supplementary Figure S3). We also found differentiation-specific interaction between the “PDLIM3/SORBS2 CRM” and the region spanning from 3 kb upstream to 4 kb downstream of the PDLIM3 TSS (Figure 3B). This interaction was observed in primary myoblasts obtained from two different individuals (one male (Figure 3B), and one female (Supplementary Figure S4)).
The physical interactions between both of these previously identified CRMs  and their corresponding target gene promoters are absent 48 h before differentiation, but then become evident at 0 h, when the cells are confluent and elongating and differentiation is induced by serum removal. These CRM-promoter interaction peaks become more pronounced at 48 h after induction of differentiation (Figure 3). These results suggest that a differentiation-induced chromosome conformation is established at these gene loci as the cells are approaching confluence and may be initiating the differentiation process even before stimulation of differentiation by serum removal. These physical interaction data parallel the previously published observation that the expression of the ACTA1 and PDLIM3 genes begins to increase from −48 h to 0 h  and our observations described above that predicted CRMs are often already bound by myogenic TFs at 0 h. Interestingly, the PDLIM3/SORBS2 CRM was observed to interact only with the PDLIM3 promoter region 23 kb downstream of the CRM but not with either of the two alternative promoters of SORBS2 located 240 kb and 385 kb upstream of the PDLIM3/SORBS2 CRM (Supplementary Figure S5). Thus, while the PDLIM3/SORBS2 CRM acts on its target promoter over a distance, the CRM selectively interacts with the closer promoter rather than the farther promoter. These observations suggest a model in which distant enhancers participate in activating gene expression by looping to contact the promoters of the genes closest to them. Alternatively, the promoter selectivity could be due to promoter sequence elements that differ between genes to specify the appropriate interaction or cell-type specific insulator architecture . Interestingly, while the PDLIM3/SORBS2 CRM did not interact with the SORBS2 promoters, it did show evidence of differentiation-specific interaction with an intronic region of one of the SORBS2 isoforms where a peak of H3K4me1 enrichment is observed (Supplementary Figure S4). This intronic region may be an additional regulatory region, and the observed contacts between these two noncoding regions raise the interesting possibility of combinatorial regulation of genes by multiple CRMs analogous to combinatorial TF input into transcriptional enhancers.
In this study we found a variety of evidence that CRMs involved in myogenic differentiation can be located distant from the proximal promoter regions of muscle genes. When 50 kb of sequence upstream and downstream of muscle genes are considered, many genomic regions with the same TF binding patterns and histone modifications as those observed in known myogenic CRMs are located more than 20 kb away from muscle genes (56% of the putative set), and 90% of them fall outside of 1-kb proximal promoters.
We found that distant CRMs were more likely than proximal CRMs to show changes in MyoG and SRF binding after induction of human muscle differentiation. This finding is in agreement with a prior study that characterized the genome-wide binding of MyoD during murine myogenic differentiation with ChIP-Seq . In that study, Cao et al. noted that 90% of promoter proximal MyoD peaks were bound by MyoD both before and after differentiation and that most of the “differentially bound” regions tended to be located further from gene TSS. “Differentially bound” peaks of MyoD binding were found in that study to be more likely to regulate differentiation-specific gene expression changes , suggesting that the distant differentially bound MyoD peaks in that study and the distant putative CRMs we observed here to be differentially bound by MyoG and SRF are likely to be differentiation-specific regulatory regions. These results are also in agreement with the previous general observation that the chromatin state at promoters is often shared across cell types while distant putative regulatory regions often have highly cell-type specific histone modifications . Beyond trends in binding dynamics with distance from gene TSS, we also observed that regions bound by three myogenic TFs tend to be located closer to promoters than regions bound by one or two TFs. This may relate to the phenomenon of TF binding “hotspots” (genomic regions where nearly every assayed TF is observed to be bound) observed in fly and worm [37, 38]. In this scenario, the regions bound within proximal promoters may be “hotspots” for TF binding, while the distant regions are not hotspots but rather are bound selectively across timepoints and cell types.
The need to search for CRMs outside of proximal promoter regions increases the difficulty of a computational search by increasing the amount of DNA sequence to be examined. However, we previously found that the PhylCRM algorithm, which searches for clusters of conserved TF binding sites, was able to predict CRMs distant from muscle genes; these predicted CRMs were subsequently validated experimentally by ChIP-qPCR and by reporter assays . Here, we found that two of these predicted and validated CRMs, the ACTA1 CRM and the PDLIM3/SORBS2 CRMs, both located more than 20 kb from muscle genes, interact with adjacent gene promoters as assayed by 3C. Additionally, PhylCRM predicted a CRM 340 kb away from the nearest muscle gene TSS near a microRNA with a known role in muscle differentiation; we found this region to be bound by the myogenic TFs MyoG, MyoD, and SRF 48 h after the initiation of differentiation. Thus, computational predictions combined with experimental measurements of cell type specific histone modifications and TF binding should provide a strong starting point for identifying active transcriptional enhancers both within and far from promoters in future studies of other cell types and processes. The DNA binding motif for the myogenic factor Mef2, in combination with MyoG and SRF motifs, was found previously to help in identifying CRMs using PhylCRM . The Mef2 antibodies available at the time of the ChIP-chip experiments presented here did not produce high-quality DNA binding profiles, but future experiments characterizing the binding of Mef2 during differentiation might provide further evidence for CRM locations.
After distant CRMs are identified, a remaining challenge is to determine which genes they regulate and the mechanisms by which they affect gene expression. The differentiation-specific interactions observed in this study by 3C for the CRMs located 20–35 kb away from ACTA1 and PDLIM3, combined with the myogenic TF binding observed with ChIP-chip and the upregulation of expression of these genes during differentiation, suggest a model in which these distal enhancers assist in activating gene expression by looping to contact the promoters of the closest genes. These are the first such distal interactions between a CRM and promoter that have been identified in mammalian muscle differentiation. Combined with the mounting evidence of CRM-promoter interactions in other cell types and systems, these data suggest that long-distance physical interaction is a general mechanism of enhancer action rather than a mechanism specific to a particular system. That these enhancer-promoter interactions begin to form at 0 h of differentiation suggests that a differentiation-specific set of looping interactions may be established before the resulting changes in gene expression and cell phenotype begin. This is consistent with findings in other systems: in mouse erythroid progenitors, the β globin LCR region exhibits interactions that are established before changes in gene expression occur , and interactions between enhancers and promoters implicated in hormone-responsive gene expression in mouse adenocarcinoma cells are already present at lower levels before hormone treatment .
When identifying target genes for CRMs outside proximal promoter regions, some previous studies have used the simplifying assumption that CRMs are likely to interact with gene promoters lying within a domain bounded by CTCF sites, on the basis that CTCF often acts as an insulator across which boundaries interactions are less likely to occur [14, 41]. This assumption matches experimental results in the case of the PDLIM3/SORBS2 CRM, which interacts with the PDLIM3 promoter in the same CTCF domain (as measured by a previous study in human skeletal muscle myoblasts ), but does not interact with the SORBS2 promoters that are separated from the CRM by several CTCF sites. However, the interaction we observed between the ACTA1 CRM and the ACTA1 promoter crosses a CTCF site, as does the interaction between the PDLIM3/SORBS2 CRM and the predicted CRM located 27 kb upstream of a SORBS2 promoter. Thus, our results suggest that the links between CRMs and their target genes cannot always be predicted accurately by CTCF sites.
How enhancer-promoter interactions are established and why they are advantageous to myogenic differentiation are as yet to be determined. It is possible that the myogenic TFs observed to bind to these interacting regions are also responsible for mediating these interactions. However, a full understanding of the interaction mechanism will require future work to search for sequence motifs recognized by other factors potentially involved in establishing chromosomal domains (CTCF, for example ), to experimentally purify factors associated with the interacting genomic regions , and to disrupt implicated DNA binding sites to test whether they are required for the physical interactions.
Although this study primarily focused on the regulatory links between distant CRMs and target genes, we observed a differentiation-specific interaction between one distant CRM (the PDLIM3/SORBS2 CRM) and another region exhibiting the H3K4me1 enhancer-associated histone modification 265 kb away in the intron of SORBS2. This interaction raises the interesting possibility of combinatorial regulation of genes by multiple CRMs, analogous to combinatorial TF input into transcriptional enhancers. Interestingly, while the PDLIM3/SORBS2 CRM interaction with the PDLIM3 promoter seems to be already present at the initiation of differentiation (0 h), the interaction with this other putative CRM does not appear until 48 h into differentiation. This may indicate a distinct role of a CRM-CRM interaction in the timecourse of myogenic differentiation. Multiple distant (>50 kb away) enhancer regions have been previously shown to interact to select the appropriate promoter for expression of the major histocompatibility complex regulator CIITA . Future work will be needed to determine the prevalence of such interactions involving multiple CRMs and what role they might play in a process such as myogenic differentiation.
This study identifies putative CRMs for myogenic differentiation distant from muscle gene promoters and demonstrates that some of these CRMs physically interact with their target genes. These results highlight the importance of looking beyond the proximal promoter to understand the transcriptional regulation of genes involved in differentiation. These distant CRM-promoter interactions may relate to the global changes in genome organization observed as cells undergo differentiation . In the future, combining computational enhancer predictions, experimental tests of distant CRMs, and global measurements of chromosome conformation  will help to further elucidate the mechanisms of gene regulation during this differentiation process.
Adult primary human skeletal myoblasts (Lonza) were grown in SkGM-2 medium (Lonza). Myogenic differentiation was stimulated by switching the culture medium to DMEM-F12 with 2% horse serum (Sigma) when the cells reached about 70% confluence. All time points referred to in this study are with respect to the time of switching to differentiation medium. The majority of the results presented here were obtained using adult male skeletal myoblasts; female skeletal myoblasts were used in one 3C experiment to confirm the reproducibility, and lack of gender-specificity, of the differentiation-specific physical interaction of the PDLIM3/SORBS2 CRM.
We designed a custom ChIP-chip array using the Agilent 1×244K oligonucleotide array format (Agilent Technologies Inc., AMADID # 015037) to achieve 25-bp resolution tiling of 100 kb surrounding each of 104 selected human genes. The genes were chosen because they were either upregulated over the timecourse of myogenic differentiation , or annotated by Gene Ontology  as being a part of the muscle sarcomere. For each of these genes, the array covers −40 kb to +10 kb relative to transcription start and −10 kb to +40 kb relative to transcription stop. Additionally, the array sequence covers: 5 kb around each of 9 negative control regions previously found to be not bound by myogenic TFs ; 5 kb around each of 35 potentially regulatory noncoding SNPs (i.e., polymorphisms predicted to disrupt myogenic TF binding sites located within a computationally predicted CRM); and 5 kb around each of 100 predicted CRMs which are distant (i.e., greater than 50 kb away from the TSS) from annotated muscle-specific or differentially expressed genes [12, 46].
ChIP experiments were performed using a modification of the ChIP protocol published described previously . Briefly, cells at various timepoints relative to differentiation were crosslinked with 1% formaldehyde for 20 min and then quenched by glycine addition. Cells were harvested by scraping and then lysed with detergents and mechanical stress. After cell lysis, crosslinked chromatin was sheared by sonication to an average length of ~600 bp, and then immunoprecipitated with relevant TF-specific antibodies (MyoD: sc-760; Myogenin: sc-576; SRF: sc-335; p300: sc-585; all from Santa Cruz Biotechnology, Inc.). A separate, mock ‘no antibody’ (i.e., protein A beads without antibody) sample was processed as a negative control. The supernatant from the mock IP was saved as ‘Input’ DNA. The purified DNA from each ChIP was amplified and labeled with Cy5 (IPed DNA) or Cy3 (Input DNA) according to a published protocol . Arrays were hybridized, washed, and scanned at 5-micron resolution according to recommendations from Agilent Technologies, Inc.
Agilent Feature Extraction software (version 188.8.131.52; Agilent Technologies Inc., Santa Clara, CA) was used to calculate log10 ratios of IP DNA/Input DNA (Cy5/Cy3) for each microarray spot in each experiment and to remove spots flagged as outliers because of debris on the slide or excessively high background. These log ratios were averaged across technical replicate experiments. ChIPOTle software  was then used to identify peaks of statistically significant TF binding enrichment for the averaged technical replicate datasets (while keeping each biological replicate separate) using a sliding window size (200 bp) that accounts for the density of tiling and the average length of the immunoprecipitated DNA fragments. A threshold on the calculated p values for each peak was set so that the resulting FWER (α) was 0.05 by the Bonferroni correction. The final set of significantly bound regions was then determined for each TF and each timepoint by finding genomic regions with binding peak overlap between biological replicates (peaks with overlapping genomic coordinates or centers separated by less than 800 bp were counted as overlapping) and without peaks in the negative control ‘no antibody’ experiment.
Distances between significant ChIP-chip binding peaks and gene TSSs were calculated using RefSeq gene annotations, excluding microRNA genes. For the purpose of the analyses described in this study, we define “muscle genes” as genes that meet at least one of the following criteria: a) previously found to be differentially upregulated (FDR 5%) over the timecourse of human myogenic differentiation (clusters 0–5 in the Warner et al. publication) ; b) annotated with one or more of the following Gene Ontology (GO) terms: myofibril, sarcomere, contractile fiber, muscle cell differentiation, muscle contraction, or structural constituent of muscle (GO database release 2011-01–06) . The locations of microRNA genes were obtained from miRBase Release 16 . All genomic coordinates are from Hg18.
The WinMine toolkit  was used to create a dependency network describing the value of ChIP binding data for certain TFs at 0 or 48 h in predicting other TF binding at 48 h. The model was learned on a 70% training set of the data using the default kappa value of 0.01 and resulted in 63% probability of correctly predicting the values in the remaining (30%) test set of data.
The PhylCRM algorithm was run with default parameters as described previously , searching for clusters of conserved MRF binding sites using the MYF motif (representing the specificities of MyoG and MyoD) from a prior study . We chose to use MRF alone for the PhylCRM predictions because this motif was more enriched among a set of genes upregulated during myogenic differentiation than any other combination of MRF, Mef2, SRF, or Tead motifs in a previous study . The phylogenetic tree used to evaluate binding site conservation in PhylCRM included human, chimp, macaque, mouse, rat, dog, cow, and opossum. From evaluation of the PhylCRM window scores for some previously experimentally tested CRMs, we set the threshold PhylCRM window score for calling a cluster of conserved MRF sites a predicted CRM at 1.5.
3C experiments were performed following previously described protocols . Briefly, approximately 5×107 human muscle cells were crosslinked in a final concentration of 1% formaldehyde and harvested at each of two timepoints: 48 hours before and 48 hours after differentiation by serum deprivation. Chromatin extracted from these cells was digested with either an AflII or BglII restriction enzyme (New England Biolabs) and then incubated with T4 DNA ligase (New England Biolabs), such that regions of chromosome contact were ligated together. After reversing the formaldehyde crosslinks, physical interactions between genomic regions were detected with PCR primers specific for each ligated interaction product. Primer pairs with melting temperatures of 56–60°C, within 2°C of each other, and of 35–50% GC content, were designed to uniquely amplify an approximately 100–300 bp region straddling the location of potential ligation between two restriction fragments. PCR-amplified products were detected by quantitative real-time PCR (qPCR) on an iCycler (BioRad) with SybrGreen Supermix for iQ (BioRad) and standard protocols. The quantification by qPCR was validated by visualization on an agarose gel and quantification using QuantityOne gel image quantification software (BioRad). The qPCR quantification results were utilized for the final analysis.
To control for differences in primer efficiency, DNA fragments generated from bacterial artificial chromosomes (BACs) spanning the genomic regions of interest were digested and randomly ligated (ACTA1: RP11–1111E20, PDLIM3/SORBS2: CTD-2559A19, CTD-2194A4, and RP11–78H20). Quantified PCR products from each 3C reaction were then normalized by the quantified results from the same PCR amplifications performed on this BAC control library. For those regions for which an ANOVA test of interaction frequencies within a given genomic region (excluding fragments within 15 kb of the fixed primer, which are likely to show high levels of interactions from random collisions ) indicated the presence of a significant peak (p<0.05) within the data, we tested the observed peaks for statistically significant differences in interaction level between timepoints using a Student’s t-test.
The raw and processed ChIP-chip microarray data are available in GEO series accession number GSE29865 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE29865).
The authors thank Job Dekker and Ye Zhan for providing training in the 3C technique, Jason Warner for assistance with array design and training in ChIP-chip experiments, Cherelle Walls for assistance with cell culture and ChIP-chip experiments, and Anthony Philippakis for assistance with PhylCRM and array design.
This work was supported by NIH/NHGRI grant # R21 HG005149 (M.L.B.). R.P.M. was supported in part by a National Science Foundation Graduate Research Fellowship.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.