Transcriptional regulatory logic of the diurnal cycle in the mouse liver

To identify DNA regulatory elements controlling diurnal transcriptional rhythms in the mouse liver, we mapped DHSs every 4 h during a full LD cycle. Specifically, C57BL/6 male mice were kept in standard 12-h light–12-h dark cycles, and four animals were euthanized every 4 h for 1 d followed by liver dissection (Materials and methods). DNase I hypersensitivity libraries were produced, sequenced, and mapped to the mouse genome using standard methods (Materials and methods). To monitor transcription activity in the same conditions, we generated chromatin immunoprecipitation followed by DNA sequencing (ChIP-seq) samples for the histone modification H3K27ac (marking active regulatory elements [33]) and resequenced previous total Pol II ChiP-seq libraries [12] at increased coverage (Materials and methods and S1 Table). Circadian clock outputs result in the rhythmic transcriptional activation of hundreds of genes, notably through binding of CLOCK:BMAL1 heterodimers [16,17]. To validate our assays, we therefore examined the known circadian output gene Dbp (S1 Movie), maximally transcribed at ZT8 [30], to determine whether cutting frequency at DHSs exhibited diurnal variation. We detected several DHSs in the vicinity of Dbp, with high intensity and narrow signals surrounded by low noise levels (Fig 1A). These DHSs near Dbp coincided well with regions identified using classical DHS mapping [30] and overlapped with BMAL1 ChIP-seq regions [17] (S1 Fig). As exemplified by DHSs near the transcription start site (TSS) of Dbp, we observed that DHSs were located in regions with lower H3K27ac signals in between H3K27ac-enriched islands, which is suggestive of TF-induced nucleosome displacement [34–36] (Fig 1B). The DNase I hypersensitivity changed diurnally, notably at the TSS (Fig 1C), where the oscillations in DNase I hypersensitivity, Pol II, and H3K27ac peaked in sync at ZT10 (Fig 1B). Moreover, all DHSs within 15 kb of the Dbp TSS displayed oscillations of varying amplitudes (maximally about 7-fold), consistent with [30], but all with the same phase as the TSS (Fig 1D), suggesting regulatory relationships between these regions and gene transcription.

We next analyzed the Npas2 gene (S2 Movie), another known clock target [37]. Npas2 is a target of RORs and peaks in the late nighttime around ZT22 [38]. We detected several DHSs along the transcribed region of Npas2 (Fig 1E), including proximal (defined as 1–10kb from a TSS) and distal (defined as >10kb from a TSS) sites. The distal sites displayed high-amplitude oscillations of DNase I signals and H3K27ac (Fig 1F). Normalized signals at the Npas2 TSS also peaked at the expected phase, with maximal signal at ZT22 for all three marks studied (Fig 1G). Finally, all DHSs associated with Npas2 (those having Npas2 as their closest TSS), including numerous distal regions, likewise cycled with phases around ZT22, some with amplitudes higher than 10-fold (Fig 1H). The examples of the Dbp and Npas2 loci suggest that our genome-wide study detected DHSs with high resolution and that the temporal patterns of DNase I cuts reflected diurnal activities of these elements.

To comprehensively map putative regulatory elements across the genome, we merged our DNase I hypersensitivity time points and performed peak finding (Materials and methods). This revealed 62,418 DHS sites, covering around 2% of the mappable genome (considering a width of 600 bp for each DHS site), which is comparable to previous studies across mouse tissues [39] (all sites and associated signals in S2 Table). Because we aimed at associating DHSs with nearby genes to infer regulatory relationships, we first discarded transcripts from Ensembl annotations that were not expressed in our samples. For this, we used histone modifications, Pol II profiles, and now also DNase I signals at transcription start and end sites of annotated transcripts to train a supervised learning method (support vector machine) that distinguishes expressed (active) from nonexpressed genes, similar to our previous work [12] (Materials and methods). To infer putative regulatory relationships, we then annotated each DHS to the nearest active TSS. Distances between DHSs and TSSs followed a bimodal distribution, with a first mode around 100 bp from the TSSs and a second 10 kb from the TSS (S2A Fig). Consistent with previous reports [40,41], one-third of our DHSs were found within 1 kb of a TSS, while almost half were located more than 10 kb from a TSS (S2B Fig). At TSSs, the genomic distributions of DNase I cuts, Pol II, and H3K27ac signals (centered on TSSs) were consistent with accessibility of DNA being determined by nucleosome displacement and Pol II complex assembly (S2C Fig) [42]. At distal DHSs, profiles of H3K27ac showed a dip in the peak center, consistent with occupation by TFs and nucleosome displacement (S2D Fig), while the weaker Pol II signals could reflect distal assembly of the transcriptional complex [43] or interactions between enhancer regions and the TSS through DNA looping [44,45].

To determine whether DHSs reflected DNA-bound transcription regulators, we searched for short windows protected from cleavage (or footprints) [46] within a ±300 bp window around the center of each DHS. This identified previously reported footprints, as illustrated for the well-characterized promoter of the Albumin (Alb) gene [47] (S2E Fig). In the promoter region of Rev-erbα (Nr1d1), the detected footprints coincided with E-boxes and high BMAL1 ChIP-seq signals (Fig 2A). Overall, the majority (70%) of DHSs within 1 kb of a TSS contained at least one footprint, while this proportion dropped to one-half for proximal (defined as DHSs within 1–10kb of a TSS) or distal (>10 kb of a TSS) DHSs (Fig 2B). Since transcribed DNA is known to be DNase I sensitive [48], the DHSs without footprints might reflect transcription. To test this, we analyzed the number of footprints in DHSs outside of promoter regions and further marked with H3K36me3, a mark coinciding with transcribed gene bodies [12,49]. Indeed, DNase I hypersensitive regions without footprints were frequently (90%) linked with highly transcribed genes (Fig 2C). Thus, DHSs at TSS seemed to contain more footprints than distal DHSs, and transcription elongation explains why some DNase I hypersensitive regions did not exhibit a footprint.

We next studied whether DNase I hypersensitivity, Pol II density, and H3K27ac quantified at the identified DHSs displayed diurnal rhythms using harmonic regression (Materials and methods). The number of cyclic regions identified at different significance thresholds clearly indicated that Pol II and H3K27ac oscillated at a larger number of DHSs compared to the DNase I signal itself, both for proximal and more distal DHSs (Fig 3A). To select rhythmically active regions, we assessed the combined rhythms of the three marks at each DHS as previously described using Fisher's combined test [12,50], which yielded 4,606 DHSs (7.3%, false discovery rate [FDR] < 0.05). For all three signals, the amplitude of the oscillations was larger at distal DHSs (the median peak-to-trough amplitude was 2-fold for DNase I and H3K27ac and higher for Pol II) compared to TSSs, and Pol II had larger amplitudes than either DNase I or H3K27ac (Fig 3B). Moreover, the peak times of the oscillations in DNase I signals were, except for some small deviations, similarly distributed as peak times in gene transcription and H3K27ac [6,11,12], with a weak evening peak around ZT10 and a marked late-night peak around ZT22 (Fig 3C). We next considered the relationships of peak times in the DNase I, Pol II, and H3K27ac rhythms. It is known that many chromatin marks exhibit diurnal rhythms that are tied to transcription [6,11,12,16], and similarly, enhancer RNAs (eRNAs) were shown to be transcribed in sync with their cognate transcripts [28]. We observed that DNase I cuts, Pol II, and H3K27ac displayed synchronous oscillations at DHSs (Fig 3D). Such relationships were maintained after removing DNase I–sensitive regions situated in the transcribed region of active genes (S3 Fig), indicating that this phenomenon was not a mere reflection of transcription [51]. To test whether the signals measured at DHSs near TSSs were temporally correlated with those at proximal or distal DHSs, we examined pairs of oscillating DHSs (FDR < 0.1, Fisher’s combined test), of which one was located near a TSS (<1 kb) and the other in an intergenic region positioned at least 2 kb and at most 20 kb from any TSS. While no pair reached statistical significance for DNase I signals (at the level of FDR < 0.1), probably reflecting that DNase I signals are noisier than the two other marks, we found 1,611 pairs oscillating for H3K27ac and 630 for Pol II. The two peak times were highly correlated, with differences within 1 h (Fig 3E), suggestive of enhancer–TSS interactions [40].

To understand how the circadian clock and the feeding–fasting cycle control diurnal gene expression in the liver, we studied mRNA expression (from [24]) and Pol II loading at TSSs in WT and Bmal1^-/- mice subject to the same night-restricted feeding regimen (S4 Fig). First, we observed that a similar number of genes oscillated in the WT and Bmal1^-/- genotypes; however, these genes had an overlap of about 30% for Pol II and 50% for mRNA (when we selected oscillations with p < 0.05). This indicates that genes with a diurnal expression differ between WT and Bmal1^-/- mice (S4A Fig). While such comparisons depend on cutoffs, stratifying by peak-to-trough amplitudes clearly showed that high-amplitude rhythms are more abundant in WT as compared to Bmal1^-/- mice (S4B Fig) and that this was more pronounced for mRNA than for Pol II loading at TSSs. For example, we found 12 genes with more than 10-fold mRNA amplitudes in WT and only 3 in Bmal1^-/- mice. Genes with Pol II or mRNA rhythms in both genotypes showed highly correlated phases, with a tendency for a slight average delay (~1 h in Pol II and less in mRNA) in the absence of a circadian clock (S4C and S4D Fig). Functional annotation using KEGG and Reactome pathways and comparison between mRNA rhythms in WT and Bmal1^-/- animals revealed that genes annotated for circadian rhythm, bile secretion, steroid metabolism, ribosome biogenesis, and SREBP signaling were enriched in the WT condition at the expected times of the day (S5A Fig). In Bmal1^-/- mice, circadian rhythms were obviously absent, but interestingly, most of the other functions were still oscillating, notably SREBP and ChREBP signaling. In general, the lipid related processes even had higher amplitudes in the Bmal1^-/- animals (S3 Table, S5A and S5B Fig). These findings are consistent with our observation that genes bound by BMAL1 or CLOCK cycle with much greater amplitudes in the WT compared to Bmal1^-/- mice. Conversely, genes bound by nutrient-responsive and systemic TFs (CREB, SREBP, HSF1) exhibited either much smaller differences (CREB) between these two conditions or no difference (SREBP, HSF1) (S6A Fig).

To identify transcriptional regulators underlying rhythmic transcription by the circadian clock and feeding–fasting cycles, we used a computational approach that combines temporal Pol II loading at TSSs in WT and Bmal1^-/- mice with annotated TF-binding sites in accessible chromatin regions as defined by our DHSs. Using DHSs and a collection of about 1,900 position-weight matrices for TF–DNA affinities (Materials and methods), we identified DNA sequence motifs that explain rhythmic Pol II patterns in WT and Bmal1^-/- mice. Briefly, we modified previously described linear models [17,52,53] to identify transcriptional activities (strictly speaking, DNA motifs) represented by phase (time of maximal activity) and amplitude (Materials and methods). In this model, motif activities are linearly combined (as in the phase vector model [54]) according to the presence of corresponding DNA motifs within nearby DHSs. This enabled us to take into account, in addition to the proximal promoter, a collection of putative regulatory regions that may control the expression of a given gene (Fig 4A). Specifically, we considered motifs in DHSs located within a certain window from active promoters, and first estimated the optimal window size according to the quality of the fit. We found that the inclusion of DHSs up to 50 kbp improved the fits in both genotypes (Fig 4B), suggesting that enhancers (represented by distal DHSs) contribute to circadian gene transcription. In WT mice (Fig 4C, S4 Table), our modeling confirmed that known circadian TFs showed the strongest activities, as reflected by the emergence of ROR response elements (RREs) [18,55] with predicted peak activity at ZT22, D-Box elements at ZT12 [56], and E-boxes around ZT8, as previously described [57]. Other motifs that had previously been associated with diurnal transcription in the liver were also identified. These included FOX motifs around ZT20 and ZT5 [58–62], the CREB motifs at ZT7 [63–67], GR motifs around ZT10 [68], SREBP motifs at ZT19 [24,69,70], HSF1 at ZT16 [21,22,26], and ETS TFs in the morning [28]. Further analysis regarding cooccurrences of these motifs at the TSSs overlapping a DHS with cycling DNase I, H3K27ac, and Pol II revealed, among several interesting associations (S6B Fig), that E-box motifs were positively associated with CREB and D-Box motifs but negatively associated with SREBP, FOX, and RRE motifs (S6B Fig). To further substantiate the E-box–CREB association, we selected DHSs overlapping BMAL1 ChIP-seq sites and considered the distribution of CREB motifs as well as the CREB ChIP-seq signal around the BMAL1 peak centers (S6C Fig). CREB motifs were indeed slightly enriched, and a CREB ChIP-seq signal was clearly present on those sites (S6D Fig).

In Bmal1^-/- mice (Fig 4D, S4 Table), activities of E-box, RRE, and D-Box motifs were not detected or greatly reduced, as expected in the absence of a functional circadian oscillator. On the other hand, TFs linked with metabolic functions, notably those associated with feeding rhythms (e.g., FOX, CREB, SREBP) were identified among the strongest contributors in the absence of a functional clock. Similarly, TFs whose activity depend on systemic signals (e.g., GR and HSF1) were also found with peak activity times that were similar in the WT and Bmal1^-/- mice. Interestingly, CREB was found among the most delayed TF activities, with a predicted delay of 6 h (S4 Table). To test this prediction, we measured nuclear levels of CREB and pCREB (activated CREB, as measured by phosphorylation on Ser 133) using western blots of nuclear extract from four independent livers every 2 h in WT and Bmal1^-/- mice (Fig 4E and S7 Fig). This showed an oscillatory pattern of CREB activity in WT mice (as previously reported in [71]), which was still oscillating in the Bmal1^-/- livers, suggesting that CREB activity is regulated by food-related signaling in clock-impaired mice subjected to a night-restricted feeding regimen. However, while we observed an average phase delay of approximately 2 h in Bmal1^-/- mice, this trend was not significant (p = 0.5, Chow test), presumably owing to interindividual variability in the feeding patterns. Indeed, such interindividual variability is common in metabolism-related, time-dependent regulatory processes, as reported for the rhythmic activation of the TORC1 and AMPK pathways [72].

We next examined how BMAL1 binding might influence cycling of activity-related signals (DNase I, H3K27ac, Pol II) at DHSs, as well as DNA accessibility more generally. To address this, we first used our previously published BMAL1 ChIP-seq data [17], which had revealed a highly nonhomogenous peak-strength distribution, showing a handful of very strong peaks associated mostly with circadian and rhythmically transcribed genes (~20) and a majority of weak sites linked to a broad range of genes without obvious rhythmic expression [17]. Consistent with this peak-strength distribution, the percentage of cycling DHSs among BMAL1-overlapped DHSs strongly increased with binding strength (S8 Fig). While we found that DNase I signals oscillate at all the strongest BMAL1 sites (this concerns only about 20 sites), the percentage of BMAL1 sites with cycling DHSs decreased to about 43% for the entire set (S8 Fig). Since tandem E-boxes were found in a large majority of the strongest BMAL1 sites [17], we hypothesized that tandem E-boxes might also promote differential DNase I signals, which we observe below (Fig 5).

Since BMAL1 was proposed to act as a pioneer-like DNA-binding factor [36], we next investigated the effect of Bmal1 knockout on DNase I signal and thus performed DHS mapping at ZT6 in our Bmal1^-/- mice [73], near the maximal DNA-binding activity of BMAL1 in WT mice. DNase I hypersensitivity at strong BMAL1 sites (detected in ChIP-seq), such as in the Rev-erbα locus, was markedly decreased (by about 5-fold) in Bmal1^-/- mice, whereas control (unbound) regions like the Gsk3 promoter showed little difference (Fig 5A). Overall, we observed a clear shift in DNase I hypersensitivity in Bmal1^-/- mice: namely, regions bound by BMAL1 in the wild type [17] showed, on average, about 2-fold reduced DNase I cuts in Bmal1^-/- animals (Fig 5B). The effect was markedly stronger when selecting BMAL1 sites more stringently—notably in BMAL1 sites with target gene expression within the ZT2–ZT10 interval—and for such sites that also contained tandem E-box sites (Fig 5C). These findings are consistent with the proposed pioneering function of the BMAL1-CLOCK complex [36]. In comparison, DNase I signals at those sites in the wild-type at the minimum of BMAL1 activity (ZT18) only showed reduction compared to ZT6 for the more stringent selections (Fig 5C), which could reflect residual DNA-binding activity of BMAL1 at the trough time [17]. To verify that the differential DNase I signals were specific to BMAL1 sites, we performed the same analysis at sites bound by the E-box–binding protein USF1 [74], which did not show similar differences between WT and Bmal1^-/- animals (Fig 5D and 5E).

Owing to the 3-D structures of protein–DNA interactions, genomic patterns of DNase I cleavage around transcription factor binding sites display factor-specific footprints [32,75–78]. We previously showed that BMAL1 binds DNA rhythmically and that strong BMAL1 binding was frequently associated with tandem E-boxes [79] separated by 6 or 7 nucleotides, which were bound by one or two CLOCK:BMAL1 dimers [17]. Here, we analyzed DNAse I footprints at BMAL1-binding sites as a function of time. Starting from BMAL1 ChIP-seq sites, we modified a “mixture model” for DNase I cuts [80] to determine the optimal boundaries of the footprints at each time point as well as the probability that the factor is bound to DNA (calculated here as the probability that the DNase I showed a footprint) for every site (Materials and methods, S1 File). We then analyzed footprints at BMAL1-binding sites containing tandem E-boxes separated by 6 bp (E1E2-sp6). At ZT6, close to the maximal DNA-binding activity of BMAL1, both E-boxes in the E1E2-sp6 motif appeared to be protected from digestion. In contrast, at ZT18, only the 5′ E-box displayed a footprint consistent with occupation by a transcription factor (Fig 6A, full time course in S9 Fig). Moreover, the footprint at ZT18 was undistinguishable from that in the Bmal1^-/- mice, suggesting that other transcription factors bind BMAL1 sites when BMAL1 activity is low. The estimated proportion of E1E2-sp6 motifs showing a footprint indicative of two CLOCK:BMAL1 dimers varied across time points, with a maximum of 65% at ZT10 and minimum of 20% in the Bmal1^-/- animals (Fig 6B). Also, the binding dynamics of BMAL1 at E1-E2-sp7 (tandem E-boxes separated by 7 bp) was largely similar to that for E1-E2-sp6, though E1-E2-sp7 had both E-boxes predominantly protected only at ZT6, suggesting spacer-specific binding dynamics (S10 Fig). In contrast, the footprints at BMAL1-binding sites with single E-boxes did not show significant changes in time or in the Bmal1^-/- mice (S11 Fig), again suggesting that other bHLH transcription factors can also bind at BMAL1 sites. In fact, footprints at DNA regions bound by the bHLH transcription factor USF1 in ChIP-seq [74] were largely similar to that of BMAL1 sites with single E-boxes, though the fraction of sites with clear footprints was reduced for USF1 compared to BMAL1 (S12 Fig).

To better understand the time-dependent footprint at BMAL1 sites and to gain insight into how the CLOCK:BMAL1 heterodimer occupies its tandem E-box–containing target sites, we used recently established 3-D protein structures of single CLOCK:BMAL1 complexes combined with molecular modeling (Materials and methods). Our models strongly support formation of CLOCK:BMAL1 heterodimers in a heterotetramer configuration at peak activity of these factors and residual binding of the dimer or other transcription factors during low-activity times. Two 3-D models of the heterotetramer configuration were constructed. In the first model, the spacing between the two E-boxes was 6 bp (sp6) (Fig 6C, S3 Movie, S9 Fig, S2 File), and in the second model the spacing was 7 bp (sp7) (S4 Movie, S10 Fig). For the model of the single CLOCK:BMAL1 complex, we used the crystal structure of the heterodimeric CLOCK:BMAL1 (Protein Data Bank [PDB] id: 4F3L) [81], into which we built the missing parts of the flexible loops. To link the single CLOCK:BMAL1 model to the E-box, we employed the complex crystal structure of CLOCK:BMAL1 basic helix–loop–helix domains bonded on the E-box (CACGTG) (PDB id: 4H10) [82]. We then superimposed the two single CLOCK:BMAL1 E-box models with the sp6 DNA and the sp7 DNA, forming the respective symmetric heterotetramer models. We found that the 6-bp spacing between the two E-boxes was optimal to establish favorable interactions between the two CLOCK:BMAL1 heterodimers (Fig 6D), involving mainly residues (e.g., K335, Y338, Q352, E380, and E384) located in the PAS-B domain of the CLOCK in a dynamic H-bond network [83]. Similarly, the 7-bp spacing seemed also able to favor a heterotetramer conformation, producing only a minor twist of 10° in the three interval base pairs. However, a conformation with base-pair spacing less than 6 or more than 7 would make complex formation difficult because of conformational constraints. Thus, the modeling results are consistent with two CLOCK:BMAL1 heterodimers binding to two E-boxes separated by 6 or 7 base pairs, and the DNase I footprints with characteristic and dynamically changing shapes suggest exchanges of different transcription factor complexes on the DNA during the diurnal cycle.

Finally, we examined temporal footprints at DHSs bound by other rhythmically active TFs: REV-ERB, HSF1, SREBP, and CREB (S13 Fig). Interestingly, and unlike what we observed for CLOCK:BMAL1, the shapes of the footprints for those factors did not change with time and were unaffected in the absence of BMAL1. However, the fraction of sites showing footprints coincided well with the maximal transcriptional activity of the different factors. For example, footprints centered on REV-ERBα–bound RREs showed the largest proportion of footprints at ZT22, which coincides with the trough activity of the REV-ERB repressors. The low percentage of bound (as detected in ChIP) RREs with footprints called by the model was low (<20%), which could reflect that nuclear receptors tend to have a low residence time and therefore display a lower DNase I cleavage-protection pattern [77]. For HSF1, the number of footprints was maximal at ZT18—approximately 4 h later than the previously reported peak activity [84]—and for the feeding-induced SREBP this number peaked during the night, as expected [24]. Lastly, high-confidence CREB-binding sites [85] showed clearly marked and invariable-width footprints throughout the 24 h in both WT and Bmal1^-/- mice, consistent with the finding that CREB activity is regulated posttranslationally on the DNA [63–66].

We mapped genome-wide DNase I hypersensitivity with 4-h time resolution in adult mouse livers. While most of the signal is likely derived from hepatocytes that constitute 80% of the liver mass, we cannot exclude that some signals, including rhythmic ones, originate from the remaining nonparenchymal cell types. Our analysis provided a comprehensive view on the dynamics of chromatin accessibility controlled by the circadian clock and feeding–fasting cycles. Overall, the identified hypersensitive regions clustered in about 60,000 DHSs. The latter tended to span several hundreds of base pairs and collectively covered 2% of the mappable genome. One-third of these regions were located near gene promoters and the remaining two-thirds were more distal from TSSs, which is consistent with what has been previously observed in mammalian cells [40]. On a genome-wide scale, 98,000 footprints were detected in about 60% of these accessible regions. Importantly, our data provided global insights into the temporal variations in DNase I hypersensitivity on the timescale of several hours to a day. Indeed, while it was previously shown that high-amplitude circadian genes such as Dbp showed nearby hypersensitive regions [30], it was not known how widespread these rhythms are across the genome. Here, we showed that thousands of DHSs exhibit rhythmic signals with peak-to-trough amplitudes that are comparable to those of Pol II signals. Accessible chromatin, as measured with DNase I hypersensitivity, is typically associated with transcriptionally active states and often reflects the presence of proteins bound to regulatory DNA elements [31,32,40]. However, we showed that DNase I–sensitive regions within gene bodies, notably in the case of highly transcribed genes, may reflect transcription elongation. As a consequence, they do not necessarily display DNA footprints such as the ones discerned in regulatory elements.

We then compared the temporal patterns of DNase I hypersensitivity with other frequently used transcriptional activity marks, in particular H3K27ac and Pol II. While DNase I signals, H3K27ac, and Pol II densities all showed abundant rhythmicity, H3K27ac abundance cycled at the largest number of DHSs, in particular at distally located sites. For both DNase I hypersensitivity and H3K27ac, the peak-to-trough amplitudes appeared higher in distal elements as compared to TSSs. Such dynamic accessibility might reflect increased protein-binding dynamics at enhancers, suggesting their potential role in controlling diurnal gene expression. This would parallel mechanisms underlying cell-type specificity, in which the modulation of histone marks and accessibility of chromatin at enhancers are among the major features associated with regulatory mechanisms [86]. The hypothesis that distal DHSs might represent enhancers for diurnal transcription was further supported by our observation that rhythms in pairs of putative enhancers and nearby TSSs showed a tight temporal correlation. In contrast to the observed delay between H3K4me3 enrichment and Pol II density, reported previously [12], no significant delays were observed between accessibility as measured by DNase I hypersensitivity and H3K27ac enrichment. This likely reflects that turnover of histone acetylation is faster than that of histone methylation [87]. We then used these temporal datasets to explore the involvement of putative enhancer regions in the cyclic recruitment of Pol II at the TSSs and subsequent transcription of the respective target genes. Our findings were consistent with a previous study on eRNA, which showed that eRNAs cluster in specific circadian phases and are correlated with Pol II occupancy and histone acetylation [6,28]. In addition, eRNA levels are correlated with the expression of nearby genes [28].

Our genome-wide study of Pol II loading and mRNA expression in WT and Bmal1^-/- mice kept under LD cycles and night-restricted feeding revealed that the number of genes exhibiting diurnal fluctuations did not drastically change in these behaviorally arrhythmic animals. This was a rather unexpected finding; however, we also found that the number of high-amplitude oscillations in mRNA accumulation was much reduced in Bmal1^-/- mice. In addition, we observed greater differences for mRNAs than for Pol II loading. Given our previous research [12] demonstrating little, if any, significant regulation of circadian gene transcription at the elongation level, these differences suggest the involvement of posttranscriptional mechanisms in controlling circadian gene expression, which is well documented [11,13,14]. All in all, our observations suggest that feeding cycles can entrain a significant set of low-amplitude transcriptional oscillations, while the circadian clock drives high-amplitude rhythms of a relatively limited number of transcripts (S4 Fig).

In this study, we accumulated compelling evidence for the contribution of distal regulatory elements in circadian transcription regulation. In fact, we observed that about 47% of DHS are located at more than 10 kb from the closest active TSS. Using penalized regression models, we predicted a collection of transcription factor binding motifs that best explain diurnal variation in transcriptional activity in both WT and Bmal1^-/- mice. Moreover, while the analysis of promoter sequences recently yielded insights into promoter architecture that favor rhythmic transcription [88], the inclusion of distal DHSs up to 50 kb improved the variance explained by our penalized linear model in WT and Bmal1^-/-. The obtained set of transcription factors that exhibited high-activity amplitudes in WT mice was similar to the one derived from a screen that used differential display of DNA-binding proteins [21]. On the other hand, comparison with Bmal1^-/- mice indicated that transcription regulators related to feeding–fasting cycles and rhythmic systemic signals were active in both genotypes, as would be expected. Among those, forkhead domain factors (FOX) have been implicated in cell cycle regulation and oxidative stress and are negatively regulated by insulin signaling [59]. Notably, FOXO1 and FOXO6, like the core clock [4], regulate the expression of key enzymes implicated in gluconeogenesis [60,89], collectively pointing towards FOX transcription factors as effectors of metabolic rhythms in liver.

We also found CREB to be among the most delayed transcription factor activities inferred by the generalized linear model in a Bmal1^-/- mouse liver. CREB is implicated in the nutrient response cycle and it regulates hepatic gluconeogenesis [62,63,68,71]. We were able to replicate the pattern of CREB activity in WT mice [71] and we showed that CREB activity is still oscillating in Bmal1^-/- mice, indicating that CREB is regulated by food-related signaling. Consistently, CREB activity during fasting was shown to be modulated by CRY1 and CRY2, which are rhythmically expressed in the liver [71]. Similarly, TFs that are responsive to systemic signals—such as HSF1 driving rhythmic transcription of heat shock proteins around ZT18 [21,22] or the GR sensitive to glucocorticoid hormones (GCs) released near the day–night transition [25,90–93]—were identified both in WT and Bmal1^-/- mice. Our identification of GR activity is consistent with the previous observation that hundreds of circadian transcripts, distinct from clock-controlled circadian genes, are under glucocorticoid control [94]. In summary, our analyses revealed an important number of genes that in the absence of the clock have higher cycling amplitude or even cycle de novo. Many of these genes have a link with metabolism, which could indicate a role for the clock in buffering nutrient-induced rhythmicity of metabolic gene expression.

CLOCK:BMAL1 complexes recruit coactivators (p300 and CBP) or corepressors (CRYs, PERs). At tandem E-box sites, these ChIP signals for BMAL1, CLOCK, and NPAS2 from [11,17] were clearly more enriched at CT8 compared to CT20 and higher at tandem compared to single E-boxes (S14 Fig). On the other hand, corepressors (CRYs and PERs) were clearly more abundant on DNA at CT20 (the effect for CRY1 was weaker but this was discussed already in [11]). Moreover, at CT20, the repressors showed a weaker signal at tandem E-boxes compared to single E-boxes. Finally, the cofactors p300 and CBP exhibited an opposite behavior: while p300 behaves similarly to BMAL1, CBP behaved like PER or CRY, suggesting a differential regulatory function.

Comparing DNase I signals between WT and Bmal1^-/- samples at ZT6 revealed that BMAL1-binding sites showed on average a decrease in DHS signals (Fig 5), which may be consistent with a pioneering function for the CLOCK:BMAL1 core clock transcription factor [36]. Moreover, our analysis of DNase I signals at nucleotide resolution revealed interesting dynamics in the shape of the footprint, which was reminiscent of our earlier proposition that strong and functional CLOCK:BMAL1 recognition elements—as with those found near the majority of core circadian clock genes—involved the binding of a dimer of CLOCK:BMAL1 heterodimers [17,74,79,83]. Here, we found that CLOCK:BMAL1 binding leaves a wide footprint spanning a tandem E-box element at the maximum activity and that this footprint shrinks to encompass only a single E-box at the minimum activity time point, resembling that footprint detected in Bmal1^-/- mice. This indicates that other E-box–binding transcription factors expressed in the liver, such as USF1, can occupy these E-box sites. These transcription factors may thereby function as placeholders to render these sites quickly accessible for CLOCK:BMAL1 heterodimers at the onset of the next circadian cycle. Indeed, USF1 has been shown to act as a nonallelic suppressor in certain mouse strains carrying a semidominant mutation of CLOCK [74]. Structural modeling of the TF–DNA complexes based on the CLOCK:BMAL1 crystal structures supported the establishment of a heterotetramer configuration at peak activity [17,74,79,83]. Other transcription factor footprints, (S13 Fig) such as for HSF, SREBP, RRE, and CREB, did not show a changing footprint shape. In the case of RRE, this might reflect a dynamic exchange of activating and repressing TFs (namely RORs and REVERBs) on their cognate binding sites [95].

All animal care and handling was performed according to the Canton of Geneva (Ueli Schibler, authorization no. 1010/3950/0) and Canton of Vaud (Nouria Hernandez, authorization no. VD2546.2 and Fred Gachon, authorization no VD 2720) laws for animal protection.

C57/BL6 male and Bmal1^-/- mice [73] 12–14-wk-old (at time of euthanasia) were housed in a 12-h light/12-h dark (LD) regimen. Starting 7 d before euthanasia, the mice were entrained to a 12 h/12 h LD regimen with ad libitum access to water. For the Pol II ChIP-seq, H3K27ac ChIP-seq, and microarray experiments, food was accessible only between ZT12 and ZT24 (restricted feeding). Mice used for DNase I-seq were entrained to a 12 h/12 h LD regimen with ad libitum access to both water and food. At each ZT2, ZT06, ZT10, ZT14, ZT18, ZT22, and ZT26, three to five mice were anesthetized with isoflurane and decapitated. The livers were perfused with 2 ml of PBS through the spleen and immediately collected. A small piece of liver tissue (approx. 100 mg) was snap-frozen in liquid nitrogen and kept at −80°C for RNA extraction. The remaining liver tissue was immediately homogenized in PBS containing 1% formaldehyde for chromatin preparation.

Mouse liver nuclei were prepared as described in [96]. Freshly prepared nuclei were suspended in ice-cold Ψ-buffer (11 mM KPO₄ at pH 7.4, 108 mM KCl, 22 mM NaCl, 5mM MgCl, 1 mM CaCl₂, 1 mM DTT) and pelleted. 5 x 10⁶ nuclei were suspended in 200 μl of Ψ-buffer supplemented with 0.2% of NP40 and 1 u/ml of DNase I (DPFF Worthington Biochemical Corporation). DNase I digestion was performed for 6 min at room temperature and the reaction was stopped by adding 200 μl of lysis buffer (50mM Tris-HCl at pH 8, 20 mM EDTA, 1% SDS, 200 μg/ml proteinase K). Protease digestion was performed overnight at 55°C. RNaseA (100 μg/ml) was then added and samples were incubated at 37°C for 1 h. DNA was then extracted twice with phenol-chloroform and precipitated with isopropanol in the presence of 0.5 M NaCl. DNAs were dissolved in 5 mM Tris-HCl pH 8. DNA from four animals were pooled, and 75 μg of DNA was loaded on 11 ml 10%–50% sucrose gradient in STE buffer (1M NaCl, 20 mM Tris-HCl at pH 8, 5 mM EDTA) and centrifuged at 30,000 rpm for 16 h at 20°C (SW 40 Ti rotor, Beckman Coulter Inc). The sucrose gradients were then fractionated, and DNA was precipitated by two volumes of ethanol in the presence of 5 μg of glycogen. Fractions containing DNA sized around 300 bp were pooled and used for Illumina library preparation.

For Bmal1^-/- animals, perfused livers were processed for chromatin preparation as described in [16]. The chromatin samples from the five mice were then pooled, frozen in liquid nitrogen, and stored at –80°C. For the ChIP experiments, the following antibodies were used: anti-RPB2 (Santa Cruz Biotechnology, sc-673-18). To determine the optimal amounts of each antibody, we performed pilot ChIP assays and determined the enrichment for a set of promoters by real-time qPCR according to [16]. A total of 1 ml of each chromatin suspension (containing about 60 μg of DNA) was incubated with 10 μg of anti-RPB2 in buffer A (20 mM Tris/HCl [pH 7.5], 150 mM NaCl, 2 mM EDTA) overnight at 4°C on a rotating wheel. 10 μl of protein A bead suspension (25% slurry in buffer A), preblocked with 10 μg/ml of salmon sperm DNA and BSA at 4°C overnight, was then added and the incubation was continued for 1 h at room temperature on a rotating wheel. The beads were then washed with dialysis buffer and ChIP wash buffer as described in [97]. Protein–DNA complexes were eluted from the beads, de-crosslinked, and treated with RNase A (and subsequently with proteinase K) as described in [16]. The DNA concentration was determined by fluorometry on the Qubit system (Invitrogen). A total of 10–12 ng DNA was used for the preparation of the library. Libraries for ultra-high-throughput sequencing were prepared with the ChIP-Seq DNA sample kit (Illumina) as recommended by the manufacturer.

For WT and Bmal1^-/- animals, H3K27ac ChIPs were performed according to the method described by [98] with a few modifications. The 100-μl chromatin aliquots were used for each IP and diluted with 900 μl of RIPA buffer (1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS in PBS at pH 7.4) and added to Dynal magnetic beads conjugated with sheep antimouse IgG dynabeads (Invitrogen, cat no: 110–31) pretreated with 3 μl of polyclonal antibody for H3K27ac (active motif, cat no: 39135) for immunoprecipitation of specific complexes. The samples were incubated overnight at 4°C on a rotator, then magnetic beads were washed seven times with lithium chloride wash buffer (100 mM Tris at pH 7.5, 500 mM LiCl, 1% NP-40, and 1% sodiumdeoxycholate) and once with 1X TE buffer (10mM Tris-HCl at pH 7.5, 0.1 mM Na₂EDTA). The chromatin complex was eluted using elution buffer (1% SDS, 0.1 M NaHCO₃) for 1 h at 65°C using Eppendorf thermomixer. The chromatin was then de-crosslinked overnight at 65°C and ChIP DNA purified using Qiagen PCR purification kit and eluted in 50 μl of elution buffer. For qPCR reaction, 1.5 μl of 1/10 diluted ChIP DNA was used. Libraries for ultra-high-throughput sequencing were prepared with the ChIP-Seq DNA sample kit (Illumina) as recommended by the manufacturer.

ChIP-seq of HSF1 was performed according to the method described by [17]. The HSF1 polyclonal antibody was from Stressgen (Enzo Life Sciences, ADI-SPA-901). For each IP, 5 μl of HSF1 antibody was used with 250 μl of precleared chromatin. A ChIP library was prepared using four independent ChIP experiments at ZT14, and one lane was sequenced to obtain about 20 million uniquely mapped reads.

Hepatic nuclear proteins were prepared as described in [27] using the NaCl-Urea-NP40 (NUN) procedure. 10 μg of the nuclear protein extracts was fractionated on an SDS-PAGE and transferred to a PVDF membrane for western blot analysis. Antibodies against CREB (Chemicon # AB3006) and phospho-CREB (pSer133) (Chemicon #AB3442) were used at 1:1,000 dilutions. Membranes were stained with naphtol blue black in order to quantify the protein loading.

At each time point, DNA sequenced reads were mapped to the mouse genome (Mus musculus NCBIM37 genome assembly [mm9; July 2007]) using bowtie through the HTS station portal (http://htsstation.epfl.ch↗) [99]. Duplicate reads were kept to avoid saturation because of high coverage of Hi-seq libraries and DNase I specificities. Quality controls, including the percentage of reads within enriched regions, indicated high overall enrichment at all time points, as about 50% of DNase-seq reads mapped to 1.3% of the genome considered to be accessible (S1 Table). Peak calling was done using ChIP-peak [100] (http://ccg.vital-it.ch/chipseq/chip_peak.php↗) on DNase I signals merged from all ZT time points with the parameters: cutoff = 100, vicinity = 400, window size = 600, threshold = 1,000. After peak calling, DNase I, Pol II, and H3K27ac signals were quantified at each time point within a window of ±300 bp around every peak center (±1 kb for H3K27ac). The values thus obtained were quantile normalized between time points for each mark. These quantile normalizations were then applied (by interpolation) to the genomic profiles in Figs 1 and 5. The wig files (S1 Fig) provided in GEO are normalized only to the library sizes.

Using ChIP-seq data for Pol II, H3K4me3, H3K36me3, and H3K27ac from [12] in the WT condition, a support vector machine classifier (SVM) was used to detect active transcripts among all Ensembl-annotated transcripts (version NCBIM37). We selected regions of interest to be ±300 bp around the TSS for Pol II and H3k4me3 and also ±300 bp around the TES for DNase I, and the last 600 bp of each transcript for the gene body mark H3K36me3. Read counts on the same strand as the transcript annotation were counted per 10 bp and quantile normalized across time. For training, a set of active and inactive transcripts were extracted, consisting of the top 10% and bottom 10% respectively, as determined by Pol II RPKM along each transcript. An SVM was trained on these active versus inactive transcripts and subsequently applied to all transcripts at each time point. Cross-validation indicated that the SVM had satisfactory false positive and false negative results for very high or very low Pol II signals (98% of test transcripts were correctly classified as either active or inactive at ZT10). Transcripts shorter than 600 bp were set to “active” if they had higher Pol II RPKM than the lower quartile of active transcripts. Transcripts were considered active when they were classified as active at minimally one time point. The active transcripts list was used to associate DHS with the closest active TSS. The annotation result provided 13,457 unique active genes linked with at least one DHS.

Rhythmicity analysis was done as previously described [12] using harmonic regression. Throughout (DHS signals, ChIP signals, mRNA expression), log₂ normalized signals were used in the harmonic regressions. The Fisher combined probability test [50] for Pol II, H3K27ac, and DNase I signals was computed to select rhythmic DHSs. This uses a chi-squared distribution with 2k (k = 3 marks) degrees of freedom. The resulting p-value was used to estimate FDRs via the linear step-up method. mRNA microarray in WT and in Bmal1^-/- mice from [12] were reanalyzed using harmonic regression.

Published datasets of ChIP-seq of CREB [63], USF1 [74], and REV-ERBα [101] in the mouse liver were quantified in our DHSs. ChIP-seq experiments such as SREBP [24] and BMAL1 [17] were included. Z-scores were computed for each ChIP-seq in each DHS. Z-scores greater than 2 were used for subsequent footprint analysis in Fig 6.

Footprints in DHSs were detected using Wellington (pyDNase library) [46] with parameters: -sh 20,36,5 -fdr 0.05 on all DNase samples concatenated. To analyze the shape of footprints, we extended a mixture model for DNase I cuts [80] to determine the optimal boundaries of the footprints at each time point as well as the probability that the factor is bound to DNA (calculated here as the probability that the DNase I showed a footprint) for every site (details in S1 File).

To identify rhythmic TF activities from temporal Pol II data, we adapted existing methods based on linear regression [102] to the circadian context [54]. Specifically, we estimated transcription factor motif activities A_f by fitting the following linear model:

where P_g denotes the 24-h component of the temporal Pol II profile (P_t) for gene g i.e., P_g = Σ_tP_te^iωt, with ω=2π24h−1. In practice, to perform linear regression with real numbers, we used real and imaginary parts as two-dimensional vectors. The matrix N_gf represents the susceptibility of gene g to the factor f and contains the motif content for factor f in all DHSs within a certain window of an active TSS (Fig 4A). To cover a large representation of TF motifs, we used FIMO [103] and scanned our DHSs using sets of position weight matrices (PWM) from JASPAR [104], TRANSFAC [105], SELEX [106], and WANG [107] (in total, ~1,900 matrices). We counted all motifs below a threshold of 10⁻⁴. The fitting was performed using the Elastic-net penalized linear regression model [102], which conveniently controls sparseness of the solution (in virtue of the L1 norm), grouping of redundant features (owing to the L2 term, this was important here, since we have a large and redundant set of matrices), and overfitting (using cross-validation). This method is available as an R package called GLMNET and uses the elastic-net penalized regression. Unless otherwise stated, we used an “alpha” (tunes the relative weights of the L1 and L2 penalties) value of 0.1. In Fig 4C and 4D, real and imaginary parts of the inferred activities A_f are plotted, showing both their amplitudes and peak activity times (phases).

For the single CLOCK:BMAL1, the crystal structure of the heterodimeric CLOCK:BMAL1 (PDB id: 4F3L) was used as an initial model [81]. In this structure, there are five flexible loops lacking density. The residues in positions 129–134 (length of 6 residues), 212–237 (26 residues), 257–275 (19 residues), and 291–309 (19 residues) were missing from BMAL1, and the residues 224–247 (24 residues) were missing from CLOCK. These missing parts were computed by Rosetta's loop modeling application (v3.5), an application that extensively remodels the backbone of the loops [108]. The loops were remodeled and refined by the CCD (Cyclic Coordinate Descent) algorithm [109]. The fragment files used by CCD were made by Robetta Server [110]. The CLOCK:BMAL1 structure, as a unique chain, was used as Rosetta input, and from the output we selected the lowest-energy loops for the single CLOCK:BMAL1 model. In order to bind the single CLOCK:BMAL1 model to the E-box, the complex crystal structure of CLOCK:BMAL1 basic helix–loop–helix domains bonded on the E-box (CACGTG) (PDB id: 4H10) was used [82]. This structure was superimposed to the single CLOCK:BMAL1 model with the UCSF Chimera visualization program (v1.5.3) [111]. In accordance with this superposition, the N-terminal helices of CLOCK and BMAL1 were replaced by the helices in the 4H10 structure from the protein data bank. The base-pair geometry of the DNA in the 4H10 structure was analyzed by the 3DNA software (v2.0) [112]. Two double-strand DNA models, spacing 6 (sp6) and spacing 7 (sp7), with sequence 5′-CACGTGAAAAAA(A)CACGTG-3′, were generated by 3DNA. The CACGTG parts were rebuilt based on the analysis of the DNA in the 4H10 structure. The spacer of 6 bp was built with the standard B-DNA backbone conformation for A–T pairs. For the final models, two CLOCK:BMAL1:E-box models were bound to the DNA models with a spacer of 6 bp (sp6) or 7 bp (sp7) by superimposing them with UCSF Chimera. In the sp6 model, we performed energy minimization for 12,500 steps with the NAMD simulation package v2.9. The model was parameterized by the AMBER force field (ff99bsc0) [113].

Wig files were generated using the bam2wig script [99] and were normalized by the number of mapped reads divided by 10⁷. DNase I signal is represented using the first position of the read alignments considered as the cutting position and without shifting strands. Pol II and H3K27ac are represented using the coverage by the whole read length after shifting forward (in the read orientation) by 80 bp and 90 bp for Pol II and H3K27ac, respectively. These wig files were then visualized on the UCSC genome browser (http://genome.ucsc.edu/↗).

Mauro Delorenzi (Swiss Institute of Bioinformatics, Bioinformatics Core Facility, Department of Oncology and Ludwig Center for Cancer Research, Faculty of Biology and Medicine, University of Lausanne), Béatrice Desvergne (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Nicolas Guex (Vital IT, Swiss Institute of Bioinformatics), Winship Herr (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Jacques Rougemont (Bioinformatics and Biostatistics Core Facility, School of Life Sciences, Ecole polytechnique Fédérale de Lausanne), Pascal Cousin (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Fabienne Lammers (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), François Mange (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Dominic Villeneuve (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Fabrice David (Swiss Institute of Bioinformatics, Bioinformatics and Biostatistics Core Facility, School of Life Sciences, Ecole polytechnique Fédérale de Lausanne), Roberto Fabbretti (Vital IT, Swiss Institute of Bioinformatics), Philippe Jacquet (Swiss Institute of Bioinformatics, Bioinformatics and Biostatistics Core Facility, School of Life Sciences, Ecole polytechnique Fédérale de Lausanne), Dmitry Kuznetsov (Vital IT, Swiss Institute of Bioinformatics), Marion Leleu (Swiss Institute of Bioinformatics, Bioinformatics and Biostatistics Core Facility, School of Life Sciences, Ecole polytechnique Fédérale de Lausanne), Robin Liechti (Vital IT, Swiss Institute of Bioinformatics), Olivier Martin (Vital IT, Swiss Institute of Bioinformatics), Aurélien Naldi (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Viviane Praz (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne, Swiss Institute of Bioinformatics), Leonor Rib (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne), Volker Vlegel (Vital IT, Swiss Institute of Bioinformatics), and Ioannis Xenarios (Center for Integrative Genomics, Faculty of Biology and Medicine, University of Lausanne, Swiss Institute of Bioinformatics, Vital IT).

High-seq Illumina sequencing data for the ChIP-seq (Pol II WT/Bmal1-/-, H3k27ac WT/Bmal1-/-, HSF1 [WT ZT14], and DNase I-seq [WT and ZT6 Bmal1-/-]) are available at GEO as the super series GSE60430.