Introduction Circadian clocks provide higher organisms with cell-autonomous and organ-based metronomes that control temporally gated and tissue-specific gene expression or metabolic programs [1]–[4]. In the liver, such programs have been implicated in detoxification [5], glucose homeostasis [6],[7], cholesterol biosynthesis [8],[9], and gating of the cell cycle [10],[11]. The mammalian clock depends on a cell-autonomous [11] core oscillator that is built around interlocked transcriptional feedback loops. These use a variety of transcriptional regulators: the basic helix-loop-helix (bHLH) PAS domain proteins CLOCK, NPAS2, and BMAL1 [12],[13], orphan nuclear receptors of the REV-ERB [14] and ROR families [15], and the DEC bHLH repressors [16]. In addition, important co-regulators such as PER and CRY proteins mediate negative feedback by repressing their own transcriptional activators, BMAL1/CLOCK [17]–[20]. Among all these regulators, the Bmal1 gene is the only single gene in the circadian network whose knockout results in arrhythmicity [21],[22]. BMAL1 functions as a heterodimeric complex, BMAL1/CLOCK, that activates transcription of its targets via E-boxes [12],[23],[24]. The DNA-binding activity of BMAL1/CLOCK is thought to cycle because of circadian changes in post-translational modifications [25],[26]. The core oscillator exerts its function by controlling temporally gated outputs, notably metabolic functions [5],[7],[27]. Transcriptional regulation of circadian output is known to occur both directly via the core clock transcription factors and indirectly, as, for example, via the PAR-bZip regulators DBP/TEF/HLF, which are themselves controlled by BMAL1/CLOCK [28]. Thus, circadian output function is controlled via a hierarchical network of transcription regulators that drives vast programs of tissue-specific gene expression both in the suprachiasmatic nucleus [29] and in peripheral tissues [29]–[34] in the mouse. Notably, these transcript rhythms cover the full range of expression phases, which thus begs the question about the mechanism behind phase-specific circadian gene expression. It has been proposed that virtually any peak expression phase can be achieved by suitably tuned regulatory sequences that integrate a small number of phase-specific core regulators [35]. Here we investigate the degree to which BMAL1 recruitment to the genomic DNA is itself rhythmic and to what extent peak binding carries phase information for downstream circadian mRNA expression. To address these questions and further dissect the hierarchical structure of circadian clock networks, we perform a time series chromatin immunoprecipitation (ChIP) analysis for the master clock regulator BMAL1 in mouse liver. This allows us to identify a comprehensive set of direct BMAL1 targets in a circadianly controlled tissue, to model the DNA-binding specificity of BMAL1 in vivo, and to determine how tightly the phase of mRNA output follows rhythmic protein-DNA interactions. Our results reveal the pervasiveness of circadian protein-DNA interactions in a mammalian tissue by showing widespread rhythmic and phase-specific binding of BMAL1 to coding and non-coding genes. This enables us to characterize the cooperative interactions of BMAL1/CLOCK complexes at tandem E-box elements (E1-E2), and to emphasize the complexity of circadian phase control that involves transcriptional and post-transcriptional mechanisms. Results BMAL1 Binds Rhythmically to Thousands of Genomic Regions in Mouse Liver To obtain a time-resolved and genome-wide map of BMAL1 target sites, we performed ChIP in mouse liver at 4-h time intervals during one light-dark cycle. Following initial testing of ChIP efficiency by quantitative PCR (qPCR) (Figure S1), two independent BMAL1 ChIP time courses were subjected to ultra-high-throughput sequencing to yield about 20 million tags per time point (Table S1) and were analyzed via a bioinformatics pipeline that combines existing and novel methods. Briefly, we used the MACS software [36] to detect regions with enriched BMAL1 binding compared to an input chromatin sample (see Materials and Methods). To efficiently reject spurious signals and accurately estimate the location of binding sites, we developed a model-based deconvolution method for ChIP combined with deep sequencing (ChIP-Seq) data (see Text S1). We identified 2,049 bona fide BMAL1 binding sites in mouse liver. Among the top 200 sites, more than 90% are significantly rhythmic (Fisher test, p 0.2) are colored. (B) BMAL1 targets in the insulin signaling pathway. (C) BMAL1 targets in the Pparα signaling pathway. These graphs were generated using KEGG Mapper (http://www.genome.jp/kegg/tool/color_pathway.html). (1.31 MB PDF) Click here for additional data file. Figure S6 Weight matrix and structure of the HMM used for sequence analysis. (A) Logo of the E-box PSWM used for autocorrelation analysis. At each position of the PSWM, the most probable letter has p = 0.96875, while the others have p = 0.03125. (B) Structure of the HMM. E1 and E2 model, respectively, the collection of hidden states of the first and second E-box. M states allow for filtering of spurious signal, namely GTGT repeats. B1 and SP represent, respectively, background and spacer states. For simplicity, the reverse complement of the HMM is not shown here. (0.35 MB PDF) Click here for additional data file. Figure S7 Pre-mRNA and mRNA measurements of longer lived transcripts. (A) mRNA transcript stability may explain lag and relative amplitude between pre-mRNA and mRNA accumulation in the Gys2, March8, and Qdpr transcripts. Experiments were performed as described in Figure 7A–7C. Approximate half-lives for March8 and Qdpr are 5.4 h and >10 h, while that for Gys2 is not available (see [B]). (B) mRNA half-lives from mouse embryonic stem cells [61] and mouse fibroblasts [62] for the genes in Figure 7A–7C and (A). When several measurements from the same cell line were available, we took the mean. (0.51 MB PDF) Click here for additional data file. Figure S8 The anti-BMAL1 antibody recognizes specifically BMAL1. Ponceau staining (A) and Western blot (B and C) of nuclear extracts (15 ug) from wild-type and Bmal1 −/− mouse liver at ZT6. The nuclear extracts were electrophoresed on a 12% SDS-PAGE gel, transferred onto a nitrocellulose membrane, and detected using anti-POLII Cter (ab817-100, Abcam) (B) and anti-BMAL1 antibodies (C). The sequence of the peptide used for the immunization is located at the C-terminal of the mouse BMAL1 protein: LEADAGLGGPVDFSDLPWPL. (3.48 MB PDF) Click here for additional data file. Table S1 Sequencing data: number of sequenced and non-redundant tags at each time point. (0.05 MB PDF) Click here for additional data file. Table S2 Functional annotation clustering of putative BMAL1 targets using DAVID tools. These annotations link the sites to the closest gene irrespective of the distance. In total, 1,551 out of 2,049 sites have a functional annotation. For details regarding the positions and binding strength of these sites, see Text S2. For the small clusters, we list the gene symbols in the most significant subcategory. (0.12 MB PDF) Click here for additional data file. Table S3 Putative BMAL1 targets with transcription factor activity (DAVID, GO:003700). Additional columns include the rank of BMAL1 binding strength, the p-value for cyclic mRNA expression (data as in Figure 6; significant values, p<0.05, are in bold), and phase of mRNA expression. For the nuclear receptors, we also indicate results for mRNA expression patterns by real-time PCR in mouse liver [27]. According to those analyses, all 18 bound receptors are expressed and 9/18 show circadian accumulation. (0.12 MB PDF) Click here for additional data file. Table S4 Enriched KEGG pathways identified with DAVID (p<0.05). The BMAL1 putative targets in the three most significant pathways are shown in Figure S5. (0.23 MB PDF) Click here for additional data file. Table S5 PSWM for the E1-E2 motif. E1 goes from position 1 to 13, position 14 corresponds to the spacer, and E2 goes from position 15 to 27. (0.05 MB PDF) Click here for additional data file. Table S6 TaqMan probes for ChIP-PCR measurements. (0.04 MB PDF) Click here for additional data file. Table S7 TaqMan probes for mRNA measurements. (0.05 MB PDF) Click here for additional data file. Table S8 Annealing primers for EMSA. (0.04 MB PDF) Click here for additional data file. Table S9 Annealing primers for transactivation assays. (0.05 MB PDF) Click here for additional data file. Text S1 Supplementary methods. (0.09 MB PDF) Click here for additional data file. Text S2 List of all BMAL1 sites with annotations and binding strength at each time point. (0.20 MB TXT) Click here for additional data file. Text S3 List of BMAL1 sites near RCGs. (0.01 MB TXT) Click here for additional data file.
