Regulation of circadian clock transcriptional output by CLOCK:BMAL1

To characterize the mechanisms by which CLOCK:BMAL1 regulates the transcriptional activity of its target genes at the genome-wide level in the mouse liver, we first generated a list of high-confidence CLOCK:BMAL1 DNA binding sites by determining the overlap between CLOCK and BMAL1 ChIP-Seq peaks in the mouse liver [3]. This analysis resulted in a list of 3217 CLOCK:BMAL1 binding sites, of which 2458 peaks can be assigned to a direct target gene (i.e., a CLOCK:BMAL1 peak located between -10kb of a target gene transcription start site and +1kb of a target gene transcription termination site; see S1 Fig, S1 Table, and methods section for details). To determine the extent to which rhythmic CLOCK:BMAL1 DNA binding contributes to rhythmic transcription activation at the genome-wide level, we used a public mouse liver Nascent-Seq dataset that characterized the levels of nascent RNA expression over the course of a 24-hr day [28]. A Nascent-Seq dataset was preferred over RNA-Seq because nascent RNA expression directly reflects transcription activation, and is unaffected by the post-transcriptional regulations that contribute to rhythmic mRNA expression in the mouse liver [3, 28, 31]. We found that only a small fraction of CLOCK:BMAL1 target genes are rhythmically transcribed (~26%; S1 Fig). Noticeably, not all rhythmic target genes are transcribed during the day, i.e., coincidently with rhythmic CLOCK:BMAL1 DNA binding (ZT02-ZT12). Indeed, 38% of the rhythmic CLOCK:BMAL1 target genes exhibit a peak of transcription between ZT12 and ZT02, out-of-phase with the rhythmic DNA binding of CLOCK:BMAL1 (n = 124 CLOCK:BMAL1 peaks) (Fig 1A–1C; S1 Fig). Importantly, our analysis also reveals that the majority of CLOCK:BMAL1 direct target genes are either arrhythmically transcribed (AR; n = 654 CLOCK:BMAL1 peaks) or not expressed (NE; n = 291 CLOCK:BMAL1 peaks) (Fig 1A–1D; S1 Fig). To determine if this result may be due to comparing samples collected in constant darkness (ChIP-Seq) and in a light:dark (LD) cycle (Nascent-Seq), we also analyzed a mouse liver BMAL1 ChIP-Seq rhythm performed under LD condition [29]. BMAL1 binding phase and ChIP-Seq signal under LD condition both exhibit a remarkably high level of similarity to those under DD conditions, and this even for the AR or NE target genes (S2 Fig). This therefore suggests that the large number of arrhythmically transcribed or not expressed CLOCK:BMAL1 target genes is not a consequence of using datasets generated under different lighting conditions. Taken together, these results indicate that the mechanisms underlying CLOCK:BMAL1-mediated rhythmic transcription of core clock genes (i.e., Per1, Per2, Per3, Cry1, Cry2, Rev-erbα, Rev-erbβ and Dbp) are not prevalent at the genome-wide level. They also suggest that the rhythmic recruitment of CLOCK:BMAL1 at its target gene promoters and enhancers is not sufficient to activate transcription for the majority of its target genes.

To investigate the mechanisms underlying CLOCK:BMAL1 heterogeneous transcriptional output, we first examined if differences in the phase, intensity or location of CLOCK:BMAL1 DNA binding might explain the differences in transcription activation. The phase of CLOCK:BMAL1 DNA binding was found to be indistinguishable between all four transcriptional output categories, as both CLOCK and BMAL1 rhythmically bind to DNA with a peak between ZT3 and ZT9 for almost all target genes (Fig 1A, 1B and 1E). We then used CLOCK and BMAL1 ChIP-Seq signal as a readout to determine DNA binding intensity, and found that both CLOCK and BMAL1 ChIP-Seq signals are significantly higher at DNA binding sites targeting the in-phase transcriptional cyclers (Rinφ) when compared to peaks targeting the 3 other groups (out-of phase cyclers, arrhythmically expressed and non-expressed target genes) (Fig 1A, 1B and 1F; Kruskal-Wallis test, p < 0.05). Remarkably, the binding intensity of CLOCK and BMAL1 at non-expressed target genes (NE) is similar to the binding intensity observed at the out-of-phase transcriptional cyclers (Ro/φ) and arrhythmically transcribed (AR) target genes, suggesting that CLOCK:BMAL1 DNA binding alone does not directly activate transcription at most of its target genes (e.g., comparisons between Fig 1D and 1F). To verify that these results are not due to the cut-offs we used to partition CLOCK:BMAL1 transcriptional output, we performed similar analyses using direct correlations between BMAL1 or CLOCK ChIP-Seq signal and the phase of rhythmic transcription, as well as by partitioning rhythmic target genes in five groups of equal sizes. These analyses confirmed our results (S3 and S4 Figs). While rhythmically transcribed target genes peaking from ZT5 to ZT13 exhibit higher BMAL1 and CLOCK ChIP-Seq signal, no differences in DNA binding signal were observed between the rhythmically expressed targets peaking from ZT13 to ZT5 and the AR and NE groups (S3 Fig). In addition, we did not find any significant correlation between either CLOCK or BMAL1 ChIP-Seq signals and the phase of DNA binding or the phase of rhythmic transcription (S4 Fig). Because CLOCK:BMAL1 peaks targeting core clock genes are enriched in the Rinφ and Ro/φ groups and exhibit higher ChIP-Seq signal than clock-controlled (output) genes, we also compared CLOCK and BMAL1 ChIP-Seq signals between groups after removing peaks targeting the core clock genes (i.e., comparing clock-controlled genes only). Whereas BMAL1 ChIP-Seq signal intensity was still significantly higher at the Rinφ target genes compared to the three other groups, CLOCK DNA binding intensity was similar between all 4 groups (Fig 1F). Our data therefore indicate that while higher BMAL1 DNA binding signal may contribute to Rinφ transcription, the different transcriptional output of CLOCK:BMAL1 target genes cannot be explained solely by differences in CLOCK:BMAL1 DNA binding intensity.

We also examined if differences in the location of CLOCK:BMAL1 DNA binding sites are associated with differences in transcriptional output by mapping CLOCK:BMAL1 peaks to either the transcription start site (TSS), gene body or extended promoter (-10 kb to -1 kb from the TSS) of their target genes. While the AR and NE groups were found to be statistically different (chi square test; p < 0.05), we did not observe any differences between the rhythmic target groups (Rinφ and Ro/φ) and the arrhythmically or not expressed groups (Fig 1G). The vast majority of CLOCK:BMAL1 peaks were located within enhancers (i.e., gene body or extended promoter), and only ~10–19% of CLOCK:BMAL1 peaks were mapped to TSS. Finally, we examined if differences in the number of genes targeted by multiple CLOCK:BMAL1 peaks were associated with differences in transcriptional output. We found that in-phase transcriptional cyclers were more frequently targeted by multiple CLOCK:BMAL1 peaks, and that conversely, non-expressed target genes were less frequently targeted by multiple peaks (S5 Fig). However, the lack of differences between the Ro/φ and AR groups indicates that the presence of multiple ChIP-Seq peaks does not directly influence the rhythmicity of CLOCK:BMAL1 target genes.

Taken together, our analysis indicates that CLOCK:BMAL1 heterogeneous transcriptional output can not be simply attributed to differences in the phase, intensity or location of CLOCK and BMAL1 binding to the DNA. While stronger DNA binding intensity may contribute to rhythmic transcription during the light phase, additional mechanisms are likely to contribute to CLOCK:BMAL1 transcriptional output heterogeneity.

Circadian repression in mammals is initiated at the beginning of the night by the recruitment of the PER/CRY repressive complex and its associated histone deacetylases and methyltransferases to CLOCK:BMAL1 on DNA [23–25, 27, 32]. Because a differential recruitment of PERs and CRYs at CLOCK:BMAL1 DNA binding sites could lead to differences in CLOCK:BMAL1-mediated transcriptional output (e.g., decreased recruitment at arrhythmically transcribed target genes, delayed recruitment of out-of-phase transcriptional cyclers, etc.), we investigated the DNA binding profile of PER1, PER2, CRY1 and CRY2 at CLOCK:BMAL1 DNA binding sites for each of the four transcriptional output groups using publically available ChIP-Seq datasets [3].

Our analysis shows that PER1, PER2 and CRY2 are rhythmically recruited at CLOCK:BMAL1 DNA binding sites with little difference between the four transcriptional output groups (S6 Fig). Maximal DNA binding for PER1, PER2 and CRY2 occur at CT12-16 for all groups, and differences were mostly observed for CRY2, where higher ChIP-Seq signal was found for rhythmically expressed target genes (S6 Fig). On the other hand, analysis of CRY1 recruitment to CLOCK:BMAL1-bound enhancers revealed more pronounced differences between all four groups. CRY1 is a potent circadian repressor that is preferentially recruited at the beginning of the light phase just prior CLOCK:BMAL1 transcription activation (i.e., CT0-4), a mechanism proposed to poise CLOCK:BMAL1 for transcription activation [3]. We found that CRY1 recruitment at CT4 is significantly higher for rhythmically transcribed target genes (both Rinφ and Ro/φ) than for arrhythmically transcribed and non-expressed genes (S6 Fig). In addition, CRY1 recruitment was significantly decreased in non-expressed CLOCK:BMAL1 target genes than arrhythmic genes at CT4. These data thus suggest that CRY1 recruitment to CLOCK:BMAL1 DNA binding sites is, in addition to its well-characterized repressive effect, linked to rhythmic transcription activation. Consistent with this hypothesis are the higher levels for Ro/φ at CT12 compared to Rinφ (S6 Fig).

Based on the mechanisms mediating the delayed transcription of the CLOCK:BMAL1 target gene Cry1 [33], a model incorporating the nuclear receptors Rev-erb (repressor) and Ror (activator), and the D-box transcriptional factors E4bp4 (also called Nfil3; repressor), Dbp, Hlf and Tef (activators) has been proposed to explain the different phases of rhythmic gene expression in the mouse liver [33, 34]. In this model, co-binding of D-box transcription factors with CLOCK:BMAL1 is proposed to delay the phase of CLOCK:BMAL1 target genes from the morning to the afternoon (i.e., from ~ZT6 to ~ZT12), and additional binding of REV-ERBs and RORs would further delay the phase of transcription to the night (e.g., ~ZT18). To test if the binding of REV-ERBs and D-box transcription factors contribute to the delay of the out-of-phase CLOCK:BMAL1 target genes, we used publicly available ChIP-Seq datasets to determine REV-ERBα, REV-ERBβ [35], and E4BP4 [36] DNA binding intensity at CLOCK:BMAL1 enhancers. We find that REV-ERBα and REV-ERBβ DNA binding, which peaks at ZT10 for all target genes [35], is significantly higher at CLOCK:BMAL1 peaks targeting genes transcribed during the night [consistent with the model proposed based on Cry1 expression; 33, 34], and no differences were observed between Rinφ, AR and NE target genes (Fig 2A and 2B; Kruskal-Wallis test, p < 0.05). The binding of E4BP4, which is maximal at ZT22 [36], was also enriched at CLOCK:BMAL1 enhancers targeting the Ro/φ genes, but to a lesser extent than what was observed for the REV-ERBs (Fig 2C). In particular, no significant difference in enrichment was observed between the Rinφ and the Ro/φ groups, perhaps because co-binding of both CLOCK:BMAL1 and D-box transcription factors drives rhythmic transcription in the afternoon around ZT12, a time used for our cut-off to differentiate the in-phase from out-of-phase transcription cyclers. In summary, our analysis indicates that the binding of REV-ERBα and REV-ERBβ (and eventually E4BP4) at CLOCK:BMAL1 enhancers may, as suggested by others [33, 34], contribute to the delayed transcription of rhythmically expressed CLOCK:BMAL1 target genes.

Our inability to detect substantial differences in CLOCK:BMAL1 DNA binding that would explain the heterogeneity of CLOCK:BMAL1 transcriptional output suggests that mechanisms other than the recruitment of core clock proteins to target gene promoters control CLOCK:BMAL1-mediated transcription. The recent finding that CLOCK:BMAL1 promotes the removal of nucleosomes when bound to DNA may represent one of these mechanisms [20]. Indeed, by mediating the removal of nucleosomes, CLOCK:BMAL1 would enable other transcription factors to access CLOCK:BMAL1 enhancers (most transcription factors bind better to naked DNA than DNA wrapped around nucleosomes).

To test if CLOCK:BMAL1-mediated nucleosome removal can contribute to the heterogeneity of CLOCK:BMAL1 transcriptional output, we examined mouse liver nucleosome signal over the 24-hr day at CLOCK:BMAL1 DNA binding sites for each of the transcriptional output groups, using a public MNase-Seq dataset (micrococcal nuclease digestion of mouse liver chromatin at 6-time points and high-throughput sequencing of mononucleosomes [20]). Our analysis reveals that nucleosome signal is rhythmic at CLOCK:BMAL1 DNA binding sites for each of the transcriptional output categories, i.e., even at CLOCK:BMAL1 DNA binding sites targeting arrhythmically transcribed and non-expressed genes (Fig 3A–3D; S7 Fig). Importantly, the phase of the rhythms is similar for all groups and minimal nucleosome signal coincides with maximal CLOCK:BMAL1 DNA binding during the light phase. Closer inspection of the levels of nucleosome signal and rhythm amplitude reveals important differences between each of the four transcriptional output categories (Fig 3E). First, the amplitude of the rhythms is significantly decreased for arrhythmically transcribed target genes. While minimal levels of nucleosome signal during the day are similar between the AR and Rinφ groups, nucleosome signal remains low during the night (i.e., when CLOCK:BMAL1 is not bound to DNA) at CLOCK:BMAL1 peaks targeting AR genes (Fig 3E). This suggests that some transcription factors may be still bound to DNA during the night in the AR group (when CLOCK:BMAL1 is not bound to DNA), thereby preventing the reformation of nucleosomes. This may promote transcription at night and thus lead to arrhythmic transcription. Second, the overall nucleosome signal is significantly lower at CLOCK:BMAL1 peaks targeting Ro/φ genes than for Rinφ genes, without any significant effect on the amplitude of the rhythm (Fig 3E). In addition, the time of minimal nucleosome signal is delayed by 4 hours between Rinφ and Ro/φ: while it coincides with the time of maximal CLOCK:BMAL1 DNA binding for Rinφ genes (ZT06), minimal nucleosome signal is observed at ZT10 for Ro/φ genes. This delayed nucleosome signal for the out-of-phase transcriptional cyclers may be explained by the significant recruitment of REV-ERBα and REV-ERBβ (Fig 2A and 2B). Indeed, CLOCK:BMAL1 has been recently proposed to facilitate circadian repression by promoting the recruitment of REV-ERBα through chromatin decondensation [37]. Thus, the increased binding of REV-ERBs at CLOCK:BMAL1 enhancers at ZT10 may promote a further decrease in nucleosome signal. Furthermore, anti-phase binding of the RORs on ROREs during the night would prevent a full nucleosome re-compaction, thereby promoting lower levels of nucleosome signal at CLOCK:BMAL1 peaks targeting Ro/φ target genes. Finally, there are no significant differences of nucleosome signal between CLOCK:BMAL1 DNA binding sites targeting in-phase transcriptional cyclers than those targeting non-expressed target genes (Fig 3E). This intriguing result suggests that although CLOCK:BMAL1 is unable to promote transcription activation at NE target genes, its rhythmic DNA binding still mediates a rhythm in nucleosome signal. One possible explanation for this result is that CLOCK:BMAL1 decondenses the chromatin to facilitate the binding of other transcription factors, but those would not be recruited at NE target genes except under specific conditions (e.g. environmental stressors), thereby preventing activation of transcription under standard conditions.

Our data indicate that CLOCK:BMAL1 rhythmic DNA binding promotes the rhythmic removal of nucleosomes at all four transcriptional output categories. We then asked if CLOCK:BMAL1 can also promote the formation of transcriptionally active enhancers. To address this question, we used public datasets [3, 36, 38] to examine the rhythmic pattern of two independent markers of enhancer activity at CLOCK:BMAL1 DNA binding sites: the post-translational modification H3K27 acetylation (H3K27ac), which positively correlates with enhancer activity at almost all enhancers and TSS [39], and the expression levels of enhancer RNA (eRNA), which are relatively short non-coding RNA molecules (50–2000 nucleotides) transcribed at active enhancer regions [40].

While mouse liver H3K27ac ChIP-Seq signal is rhythmic and high during the light phase at CLOCK:BMAL1 DNA binding sites targeting in-phase transcriptional cyclers (consistent with CLOCK:BMAL1 directly facilitating the acetylation of H3K27; Fig 3F and S8 Fig), significant differences were observed at CLOCK:BMAL1 DNA binding sites targeting the other 3 transcriptional output categories. Rhythmic H3K27ac rhythm is delayed for the out-of-phase transcriptional cyclers, and the amplitude of H3K27ac rhythm is significantly dampened at CLOCK:BMAL1 DNA binding sites targeting arrhythmically transcribed genes (Fig 3F and S8A Fig). Remarkably, levels of H3K27ac are close to background levels at CLOCK:BMAL1 peaks targeting non-expressed genes. Given that CLOCK:BMAL1 rhythmically binds to relatively similar levels for all four transcriptional output categories, our analysis suggests that CLOCK:BMAL1 DNA binding does not directly contribute to the acetylation of H3K27.

To extend on this observation, we then examined another marker of enhancer transcriptional activity by assessing eRNA expression at CLOCK:BMAL1 DNA binding sites using a publicly available GRO-Seq dataset [36]. The analysis confirmed the results obtained with H3K27ac (Fig 3G). Rhythmic eRNA expression is only observed at CLOCK:BMAL1 enhancers targeting rhythmically transcribed genes, and eRNA expression at enhancers targeting non-expressed genes is dramatically decreased to levels close to background (Fig 3G). Importantly, these differences in eRNA expression between the four CLOCK:BMAL1 transcriptional output categories are further corroborated by similar variations in RNA Polymerase II (Pol II) ChIP-Seq signal at CLOCK:BMAL1 enhancers (S8C Fig). Altogether, our analysis therefore demonstrates that, contrary to what has been typically found for core clock genes, CLOCK:BMAL1 DNA binding is not sufficient to promote the activation of its enhancers. Instead, our results suggest that CLOCK:BMAL1 rhythmically opens the chromatin to facilitate the binding of other transcription factors at its enhancers, and that the nature of these transcription factors (e.g., activators, repressors) significantly contributes to CLOCK:BMAL1 transcriptional output.

To test our hypothesis that transcription factors bind at CLOCK:BMAL1 enhancers to contribute to their transcriptional activity and thereby impact on CLOCK:BMAL1-mediated transcription, we assessed if transcription factors were differentially recruited at CLOCK:BMAL1 DNA binding sites within each transcriptional output group. To this end, we performed a DNA binding motif analysis using HOMER Software Suite that we further validated using mouse liver transcription factor ChIP-Seq datasets.

As expected, the motif analysis revealed that CLOCK:BMAL1 DNA binding motif (e-box of the sequence CACGTG) is highly enriched at CLOCK:BMAL1 enhancers for all transcriptional output categories (Fig 4A). Surprisingly however, we found that motifs for liver-specific transcription factors (e.g., Cebp, Hnf1, Hnf4 and Hnf6) were also enriched for all four transcriptional output categories, and thus even at CLOCK:BMAL1 enhancers targeting non-expressed genes (Fig 4B, S9 Fig, and S2 Table). On the contrary, motifs for ubiquitous transcription factors (u-TFs; broadly expressed transcription factors with a transcriptional activity regulated by external factors) were almost always enriched for specific CLOCK:BMAL1 transcriptional output groups (Fig 4C, S9 Fig, and S2 Table). For example, CRE motif was enriched at all CLOCK:BMAL1 enhancers except those targeting out-of-phase transcriptional cyclers, and FXR motif was enriched at all CLOCK:BMAL1 enhancers except those targeting out-of-phase transcriptional cyclers. Noticeably, the motifs for NF-κB [which binds to DNA and becomes transcriptionally active upon infection and inflammation; 41, 42], and CTCF [which establishes discrete functional chromatin domains by promoting DNA looping; 43, 44, 45] were enriched at enhancers targeting non-expressed genes.

To assess the relevance of this difference of motif enrichments between tissue-specific (ts-TFs) and u-TFs, we determined the DNA binding pattern of several transcription factors at CLOCK:BMAL1 enhancers in the mouse liver using publicly available transcription factors ChIP-Seq datasets [46–51]. This in vivo analysis largely confirmed the computational motif analysis: most liver-specific TFs were found to bind at CLOCK:BMAL1 DNA binding sites independently of the transcriptional output, whereas u-TFs were more specifically enriched in specific CLOCK:BMAL1 transcriptional output groups (Fig 4D, E, S10 Fig, and S3 Table). For example, Hnf4a and Hnf1 are the only liver-specific TF to exhibit a differential binding between CLOCK:BMAL1 transcriptional output groups of the six TFs tested (Foxa1, Foxa2, Hnf1, Hnf4A, Hnf6 and Cepba). Conversely, all twelve u-TFs investigated exhibit DNA binding differences at CLOCK:BMAL1 enhancers between categories of transcriptional output (Fig 4D and 4E, S10 Fig). Although each u-TF bound to different subsets of CLOCK:BMAL1 enhancers, u-TF recruitment was generally higher in rhythmically expressed target genes and lower in non-expressed target genes compared to the arrhythmic CLOCK:BMAL1 target group. (Fig 4D and 4E, S10 Fig, and S3 Table). To further characterize the differences in TF DNA binding between CLOCK:BMAL1, u-TFs and ts-TFs, we computed a TF DNA binding variability index by calculating the standard deviation of the ChIP-Seq signal between the 4 CLOCK:BMAL1 transcriptional output groups (see methods for details). We found that the DNA binding variability at CLOCK:BMAL1 peaks is comparable between CLOCK:BMAL1 and ts-TFs, whereas there is significantly more variability for u-TFs than for CLOCK and BMAL1 when peaks targeting core clock genes are removed from the analysis (Fig 4G). While there are variability index differences among ts-TFs and u-TFs, this analysis further supports our finding that u-TF recruitment at CLOCK:BMAL1 peaks is globally more variable than for ts-TF (Fig 4G).

Altogether, our data indicate that the mechanisms by which CLOCK:BMAL1 regulates transcription of clock-controlled genes differ from the well-characterized CLOCK:BMAL1-mediated regulation of core clock gene expression. Specifically, our data show that although CLOCK:BMAL1 mediates rhythmic nucleosome removal at its enhancers, it is not sufficient to generate an active enhancer or drive rhythmic transcription. We thus propose a model whereby CLOCK:BMAL1 regulates transcription of clock-controlled genes by rhythmically opening chromatin to facilitate the binding of other transcription factors at its enhancers (Fig 5A). This possibility is supported by results showing that nucleosome signal is rhythmic at the DNA binding sites of several TFs when those sites are located close to a CLOCK:BMAL1 peak, and not rhythmic when CLOCK:BMAL1 binding is absent (S11 Fig). Consequently, the transcriptional activities of these transcription factors would dictate the transcriptional outcome of clock-controlled genes rather than CLOCK:BMAL1 (Fig 5A). For example, binding of positive transcription factors along with CLOCK:BMAL1 would activate enhancers and lead to transcription activation during the day, whereas binding of transcriptional repressors (e.g., REV-ERBα/β) would inhibit CLOCK:BMAL1 enhancer activity and thereby contribute to rhythmic transcription peaking during the night, in anti-phase with CLOCK:BMAL1 DNA binding (Fig 5A). If no transcription factors are recruited (e.g., inducible transcription factors), CLOCK:BMAL1 enhancers remain silent and target genes are not expressed or are arrhythmically expressed (Fig 5A). Arrhythmically expressed genes at CLOCK:BMAL1 enhancers may also have positive transcription factors bound at all times overriding the absence of CLOCK:BMAL1 DNA binding at night (see result section about rhythmic nucleosome signal and Fig 3E). Our results also suggest that u-TFs regulate CLOCK:BMAL1 transcriptional output more prevalently than ts-TFs. It may be that ts-TFs facilitate the binding of CLOCK:BMAL1 at tissue-specific enhancers rather than contributing to CLOCK:BMAL1 transcriptional output (see discussion). To validate this model experimentally, we investigated how i) Bmal1 knockout, and ii) changes in environmental conditions (that alter u-TFs transcriptional activities) affect CLOCK:BMAL1 transcriptional output.

If the activity of u-TFs contributes to CLOCK:BMAL1 regulation of clock-controlled gene transcription, then, a knockout of Bmal1 (which eliminates CLOCK:BMAL1-mediated transcription [52]) should differentially affect the expression of CLOCK:BMAL1 target genes. Specifically, target gene expression levels in Bmal1-/- mouse should be arrhythmic and low for the in-phase transcriptional cyclers (no recruitment of positive transcription factors at CLOCK:BMAL1 enhancers), while they should be arrhythmic and high for the out-phase transcription cyclers (no recruitment of repressors at CLOCK:BMAL1 enhancers). These effects should also be more obvious during the light phase when CLOCK:BMAL1 binds to DNA. In addition, Bmal1 knockout should have a reduced effect on arrhythmically and non-expressed target genes. These predictions were confirmed by analyzing a public dataset that characterized the genome-wide effect of Bmal1 knockout in the mouse liver using RNA-Seq of rRNA-depleted total RNA (Fig 5B and 5C) [53]. For both intronic and exonic RNA-Seq signal, the expression of Rinφ genes in Bmal1-/- mouse liver is at the trough level of wild-type mice, and at peak levels in Ro/φ genes. Moreover, Bmal1 knockout does not significantly affect the expression levels of arrhythmic and non-expressed CLOCK:BMAL1 target genes (Fig 5B and 5C).

Our model also predicts that the transcriptional output of CLOCK:BMAL1 target genes can be altered by environmental changes that affect u-TF DNA binding capacity. External signals that activate or repress the binding of u-TFs are predicted to impact CLOCK:BMAL1 cooperation with other transcription factors, and thereby change the transcriptional output of CLOCK:BMAL1 target genes. For example, signals that enable the recruitment of positive transcription factors at CLOCK:BMAL1 enhancers could increase the amplitude of rhythmic transcription and/or initiate the rhythmic expression of target genes that are arrhythmic under control conditions. Conversely, signals that inhibit the binding of transcription factors that normally cooperate with CLOCK:BMAL1 could blunt the rhythmic expression of some CLOCK:BMAL1 target genes. To test this hypothesis, we analyzed how fasting, which is known to affect the transcriptional activity of many u-TFs [54, 55], alters CLOCK:BMAL1 target gene expression in the mouse liver using a public dataset [56]. Strikingly, while the expression of Clock, Bmal1 and several direct rhythmic target genes (e.g., Phf17, Slc16a2) are unaffected by fasting, some other targets exhibit a significantly altered gene expression profile (Fig 5D–5G, S12 Fig for additional examples). For example, some rhythmic target genes become arrhythmically expressed under fasting (e.g., Sgk2, Flcn) while other targets exhibit an increased amplitude of expression (e.g., Gnat1, Gm129) (Fig 5E). Remarkably, some direct CLOCK:BMAL1 target genes that are arrhythmically or not expressed under ad libitum condition become rhythmically expressed under fasting condition (Fig 5F and 5G). Because not all CLOCK:BMAL1 target genes are equally affected by fasting, these results cannot simply be explained by a global change in CLOCK:BMAL1 transcriptional activity under fasting condition. One possibility is that fasting enhances or represses the transcriptional capabilities of several u-TFs that cooperate with CLOCK:BMAL1, thereby altering the transcriptional output of many direct CLOCK:BMAL1 target genes. Similar results were found by analyzing a public dataset investigating the effect of high-fat diet on rhythmic gene expression in the mouse liver (S13 Fig) [57].

Unless notified below, publically available datasets used in this paper were downloaded from the NCBI or EBI websites in either sra or fastq formats (see Table 1 for accession numbers). Files in sra format were converted to fastq files using the sratoolkit (version 2.3.5–2). Fastq files were mapped to the mouse genome (version mm10) using bowtie2 [96] or tophat2 [97]. For all datasets, we only considered reads that mapped uniquely to the mouse genome (i.e., one unique genomic location). Datasets were further filtered to remove duplicated reads using samtools (rmdup function) or a custom-made script. Additional information is provided for each dataset as Supplementary Materials and Methods.

Genomic locations of CLOCK and BMAL1 DNA binding sites in the mouse liver provided in the original paper (supplementary S2 Table) [3] were used to generate our list of high confidence CLOCK:BMAL1 DNA binding sites. Genomic locations were converted to the mm10 version of the mouse genome using UCSC genome browser liftOver tools, and processed as indicated in S1 Fig to generate our list of high confidence CLOCK:BMAL1 DNA binding sites. Overlap between CLOCK and BMAL1 ChIP-Seq peaks was determined using bedtools (intersectBed) and coordinates from BMAL1 ChIP-Seq datatsets were further kept to generate a list of 3217 CLOCK:BMAL1 peaks. We also used the original data provided by the authors in their S2 Table to assign CLOCK:BMAL1 peaks to their putative target genes (original analysis performed using HOMER tools). In particular, we defined a CLOCK:BMAL1 target gene as a gene with at least one CLOCK:BMAL1 peak located between -10kb of the transcription start site and +1kb from the transcription termination site. Using this criteria, 2458 CLOCK:BMAL1 peaks were assigned to a target gene, and the remaining 759 peaks were assigned as an intergenic CLOCK:BMAL1 DNA binding site.

The 2458 CLOCK:BMAL1 peaks assigned to a target gene were then parsed based on the transcription profile of their target genes using the Nascent-Seq datasets from Menet et al., 2012 [28]. We directly used the original Nascent-Seq expression values and the assessment of their rhythmic expression from the original paper without performing new analysis. Details on how genes were determined to be rhythmically transcribed are provided in Supplementary Materials and Methods. Using these data, CLOCK:BMAL1 peaks were parsed into 4 different categories of transcriptional output (see also S1 Fig):

The list of the 3217 CLOCK:BMAL1 peaks parsed into the different transcriptional output categories is provided in S1 Table. The phase of rhythmic CLOCK:BMAL1 DNA binding, ChIP-Seq signal, and genomic location of CLOCK:BMAL1 DNA binding sites were retrieved from the Koike et al., 2012 original paper supplementary S2 Table [3] and processed to generate the analysis presented in Fig 1.

ChIP-Seq, MNase-Seq and GRO-Seq signal was retrieved from bam files containing uniquely mapped reads (and duplicated reads removed) at CLOCK:BMAL1 enhancers using custom-made scripts [20]. Specifically, signal was retrieved at:

and normalized to the sequencing depth. Differences in the window size were calculated based on the width of the ChIP-Seq signal at CLOCK:BMAL1 DNA binding sites (e.g., H3K27ac ChIP-Seq signal is significantly wider than any transcription factor ChIP-Seq signal). Because we aimed at assessing the role of CLOCK:BMAL1 in removing a nucleosome at its DNA binding site, we chose a narrower window size of 150bp (see Fig 3A–3D). All analyses were performed at individual CLOCK:BMAL1 ChIP-Seq peaks, and this even for peaks targeting the same gene. Data presented in S4 Fig examined the role of multiple peaks targeting the same gene on BMAL1 and CLOCK ChIP-Seq signals. For all datasets, ChIP-Seq signal is displayed as the number of reads/bp per 100,000,000 reads.

Enhancers lying into CLOCK:BMAL1 target gene loci (-10kb from the transcription start site to +1kb from the transcription termination site) were identified using a public mouse liver DNAse-Seq dataset (see above) [98] and bedtools (intersectBed function). Enhancers were then parsed based on the presence or not of a CLOCK:BMAL1 ChIP-Seq peak (3155 out of the 3217 CLOCK:BMAL1 ChIP-Seq peaks are located into a DNaseI hypersensitive site). Because a majority of the 104,556 DHS peaks only displayed low levels of ChIP-Seq (transcription factors, Pol II, H3K27ac) and GRO-Seq signals [as shown in the ENCODE project, 99, 100], we filtered the number of DHS lying into a CLOCK:BMAL1 target gene by only considering those being into the top 40,000 DHS list (based on DNase-Seq signal), obtaining the following number of DHS peaks:

H3K27ac and Pol II ChIP-Seq signals, as well as GRO-Seq (eRNA) signal, were retrieved at those DHS sites (as well as those overlapping with a CLOCK:BMAL1 peak) using the DHS peak coordinate and normalized to 100,000,000 reads. Signal was then normalized to the coordinate length (in bp) to obtain the signal displayed as reads/bp per 100,000,000 reads. The coordinates used were, for the same reason as above for CLOCK:BMAL1 DNA binding sites:

Because our analysis revealed the existence of small but significant overall variations of H3K27ac and Pol II ChIP-Seq signal between time points (see S8 Fig), we further normalized the datasets by performing either a mean normalization (H3K27ac) or a ranking analysis (Pol II). For H3K27ac ChIP-Seq datasets [3, 38], averaged H3K27ac signal was calculated at the top 40,000 DHS peaks (the top 40,000 DHS peaks concentrate the majority of TFs DNA binding sites; peak center ± 1 kb; total of 104,556 total DHS peaks; dataset from Ling et al., 2010 [98]) for each time point. This averaged signal was then used to normalize the raw H3K27ac ChIP-Seq signal, by calculating for each time point the ratio between H3K27ac signal for each peak and this averaged signal (see S8 Fig). Pol II ChIP-Seq dataset [22] were normalized by performing a ranking normalization (method similar to a quantile normalization). To this end, Pol II ChIP-Seq signal was calculated at all 104,556 DHS peaks (peaks mapped in Ling et al., 2010 paper [98]), and sorted based on the ChIP-Seq values. The raw values for each DHS peak were then normalized using the sorted averaged ChIP-Seq signal at each of the 104,556 ranks for all time points.

Motif analysis was performed at CLOCK:BMAL1 enhancers (original peak coordinates) for each of the transcriptional categories using the findsMotifGenome.pl script from the HOMER suite. Parameters were as the following: -size given–len8. The resulting table was sorted by the q-value and a q-value less than 0.05 was considered significant. Percent enrichment (percent of target sequences with motifs / percent of background sequences with motif) was then calculated for motifs found to be significant in at least one of the CLOCK:BMAL1 transcriptional output category. Results of the motif analysis are provided as. S2 Table

To determine the variance of each TF DNA binding (CLOCK, BMAL1, ts-TFs and u-TFs) between the four CLOCK:BMAL1 transcriptional output categories, we computed a TF DNA binding “variability index” based on the analysis performed in. The variability index was calculated by summing up the standard deviation of the ChIP-Seq signal between the 4 transcriptional output groups, which was calculated for each decile (0.1 to 0.9) and normalized to the averaged signal for each decile (the standard deviation is higher for upper deciles because ChIP-Seq signals are higher). This index reflects differential DNA binding strength between groups, as similar binding between the 4 groups results in small standard deviation values for each decile, and thus a small variability index. Conversely, differences in DNA binding signal between groups result in larger standard deviation values and thus a larger variability index. S10 Fig

To determine if the results described in this paper are specific to CLOCK:BMAL1, we also performed an analysis on genes not targeted by CLOCK:BMAL1, but exhibiting similar profiles of expression to the 4 CLOCK:BMAL1 transcriptional output categories (Rinφ, Ro/φ, AR and NE). To this end, 125 genes were randomly selected for each of the 4 groups, using criteria similar to those used to define CLOCK:BMAL1 transcriptional output (see above). Levels of expression for each group were not significantly different to those of CLOCK:BMAL1 target genes (Kruskal-Wallis test). Nucleosome signal, H3K27ac ChIP-Seq signal, Pol II ChIP-Seq signal, eRNA expression, tissue specific and ubiquitous transcription factor ChIP-Seq signal were all calculated as described above for CLOCK:BMAL1 target genes. Statistical analysis was also performed similarly to CLOCK:BMAL1 transcriptional output.

Statistical analysis was done using JMP, Version Pro 12.0.1. SAS Institute Inc., Cary, NC, 1989–2007. Differences in sequencing signal, represented in the boxplot graphs, were analyzed for statistical enrichment using the nonparametric Kruskal-Wallis test. Rhythmic analysis of nucleosome signal and ChIP-Seq signal was performed using a Fourier analysis (Fig 3A–3D) (see Supplementary Materials and Methods for details). Differences in the amplitude of nucleosome signal rhythm (Fig 3E) were analyzed using a 2-way ANOVA. Differences in CLOCK:BMAL1 ChIP-Seq peaks genomic location were analyzed using a chi-square test (Fig 1G), and differences in the number of CLOCK:BMAL1 peaks per target genes (S5A Fig) were analyzed by a Fisher’s exact test. Differences were considered significant when p < 0.05.

We did not generate new sequencing datasets. All sequencing datasets used in our manuscript are publicly available and were downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) website. The accession numbers for these datasets are available to view in Table 1.