Nascent-Seq reveals novel features of mouse circadian transcriptional regulation

To address the regulation of genome-wide transcription, we analyzed mouse liver nascent RNA expression, that is, RNA being transcribed by RNA Polymerase II (Pol II) prior to 3′ end formation (from 12 independent samples in LD, 6 time points per day done twice; see analysis of rhythmic transcription in mouse liver section). To this end, nascent RNA was extracted from purified nuclei using the high salt/urea/detergent buffer originally described by Wuarin and Schibler (1994). A very similar sequencing strategy was recently applied by Smale, Black and colleagues to macrophage nascent RNA (Bhatt et al., 2012). We then prepared illumina libraries with standard protocols for high-throughput sequencing (Nascent-Seq; Carrillo Oesterreich et al., 2010; Khodor et al., 2011). Removal of rRNA was unnecessary as approximately 65–70% of the sequences uniquely mapped to the genome (Table 1).

Seventy six percent of these uniquely mapped sequences map to introns (Figure 1A). This contrasts dramatically with more conventional RNA-Seq; it has minimal intron reads as it assays polyadenylated (pA) RNA and therefore predominantly mature (spliced) mRNA (Figure 1A). Intronic Nascent-Seq reads are also more abundant in mouse compared to Drosophila (76% vs 45%), reflecting longer intron size and less efficient mouse co-transcriptional splicing (Khodor et al., 2012; Khodor et al., 2011). Many genes exhibit a 5′ to 3′ gradient in the Nascent-Seq dataset, presumably reflecting nascent RNAs of different lengths attached to elongating Pol II (Figure 1B). In addition, Nascent-Seq signals frequently extend past the polyadenylation site, reflecting RNA not yet cleaved by the cleavage/polyadenylation specificity factor (CPSF) and/or RNA molecules still associated with Pol II after cleavage but prior to degradation by the 5′ to 3′ exoribonuclease Xrn2 (Figure 1C). These features are absent from standard RNA-Seq data, and indicate that Nascent-Seq predominantly detects nascent RNA molecules attached to elongating Pol II (Figure 1B,C).

Another feature was apparent in the comparison of Nascent-Seq and RNA-Seq: gene expression was often different between the two datasets. Some genes have a high ratio of Nascent-Seq to RNA-Seq signal (e.g., B4galt1, Figure 1B), whereas others have a low ratio (e.g., Bag1, Figure 1B). Genes with a high ratio (top 10% of all genes) are dramatically enriched for specific functions, namely, non-coding RNA (ncRNA), G-coupled protein receptors, regulation of transcription and chromatin organization (Figure 1D,E). Genes with a low ratio are enriched for genes involved in ribosome function and mitochondrial respiration (Figure 1D,E). Because these genes are associated with short or long mRNA half-lives, respectively, we compared the Nascent-Seq to RNA-Seq ratios of genes with their published mRNA half-lives (Sharova et al., 2009). Not surprisingly, genes with a high ratio have relatively short half-lives, and genes with low ratios have longer mRNA half-lives (Figure 1F). These data indicates that mRNA stability contributes to the wide range of Nascent-Seq to RNA-Seq ratios.

To identify genes that are rhythmically transcribed, we performed two independent six time-points rhythms of mouse liver Nascent-Seq. We found that some genes exhibit very high amplitude rhythms, with no detectable signal at low time points (e.g., Npas2; Figure 2A). About 15% of expressed genes manifest transcriptional rhythms (p<0.05: 6.3%, strong rhythms; 8.9% medium-strength rhythms; see ‘Materials and methods’ for analysis details; Figure 2B). Phases of maximal transcription are heterogeneous yet not uniformly distributed; very few genes peak at ZT16-20 in the mid-late night (Figure 2C,D).

Protein-coding genes dominate the cycling Nascent-Seq dataset (>85%, data not shown), but it also contains genes encoding ncRNAs, for example, pri-miRNAs and long non-coding RNAs (lncRNA) (Figure 2E,F; Figure 2—figure supplement 1–6). As previously described (Gatfield et al., 2009), pri-miRNA 122a is robustly rhythmic (Figure 2E). However, rigorous quantitation of rhythmic ncRNA transcription is precluded by the poor annotation of these transcription units.

To address the relationship of cycling transcription to cycling mRNAs, we assayed two independent six time-points profiles of pA RNA by RNA-Seq. The comparison with Nascent-Seq was restricted to genes that were sufficiently expressed in both datasets (n = 5454 genes; see ‘Materials and methods’ for details). Using the identical statistical analysis and cut-offs to determine rhythmicity, the fraction of rhythmic mRNA was higher than the fraction of rhythmic nascent RNA (22.1% and 15.1%, respectively; Figure 3A).

There was also a notably poor overlap between the two rhythmic gene sets: only 41.6% of rhythmically transcribed genes also manifest rhythmic mRNA expression (R-R gene set; 342/822; Figure 3A–C). However, the mRNA phase of these R-R genes was highly correlated to the nascent RNA phase (r = 0.92; Figure 3C), and more than half (57%) of these genes exhibit a phase difference of less than 2 hr (195/342; Figure 3D). The amplitudes of the nascent RNA rhythms were also correlated with those of mRNA (r = 0.76; Figure 3E), indicating that transcriptional regulation dominates these R-R rhythms. Not surprisingly, almost all clock genes (Figure 4) and well-characterized clock-controlled genes (e.g., Nocturnin, Por, Alas1, Upp2, Usp2, Inmt, Nfil3, etc; Figure 3—source data 1) are in this R-R gene set.

The other 58.4% of rhythmically transcribed genes (480/822) do not show robust mRNA expression (R-AR gene set; Figure 3A,F). A simple explanation is that the mRNA half-lives of these R-AR genes are relatively long and therefore mask the transcriptional oscillations. However, these genes do not have altered nascent RNA to mRNA ratios compared to the whole genome (Figure 5A) or reported longer mRNA half-lives (assessed using the dataset from Sharova et al., 2009; Figure 5B). These considerations suggest that other mechanisms account for the poor mRNA oscillations of this gene set (e.g., the rhythmic Nascent-Seq signal of 25 R-AR genes results from rhythmic transcription of an adjacent gene that reads into the R-AR gene; Figure 5C).

The opposite comparison is based on genes with rhythmic mRNA expression, of which only 28.4% (342/1204) have rhythmic transcription that meets the cycling criteria (Figure 3A,G). This surprising conclusion was similar when the analysis was restricted to genes with the strongest mRNA rhythms (121/435) and indicates that most cycling mRNAs (862/1204) likely undergo post-transcriptional regulation. This might include the circadian regulation of nuclear RNA processing, export, translational regulation and/or mRNA turnover, as described for the few circadian genes shown to be regulated post-transcriptionally (Kojima et al., 2011; Staiger and Koster, 2011; Zhang et al., 2011). Gene ontology analysis of this arrhythmic transcription-rhythmic mRNA (AR-R) gene set did not reveal any striking enrichment of particular gene functions (Figure 5D).

There is a peculiar feature of the large number of genes within this AR-R category: visualization of RNA expression with heatmap indicates that many of these transcriptionally ‘arrhythmic’ genes manifest elevated transcription at times that match their cycling mRNA peaks (Figure 3G; note that the heatmaps show that the transcription peak of many genes matches the peak phase of mRNA cycling). Further inspection of all individual expression profiles confirmed this correlation: many genes have matching peak phases despite large variations in nascent RNA expression (Figure 6A).

To substantiate this variability in nascent RNA expression, we calculated the standard deviation normalized to the mean (SD) for every gene using the two 6 time-points profiles as 12 independent samples. The reasoning was that low variability expression between time points should result in a low SD. The SD was noticeably higher in the Nascent-Seq dataset than in the RNA-Seq dataset, indicating that the transcription of most AR-R genes is indeed variable compared to the mRNA-Seq dataset (Figure 6C; note that most points are to the right of the diagonal line). Higher variability of nascent RNA expression is also observed for other categories, including the AR-AR set (data not shown), suggesting that it is a common feature of transcription vs mRNA comparisons independent of circadian considerations.

Yet a higher variability of transcription (higher SD) correlates with rhythmic mRNA expression within all arrhythmically transcribed genes (n = 4632), suggesting that transcriptional variability generally contributes to the generation of rhythmic mRNA expression (Figure 6D). This relationship is independent of nascent RNA expression levels, indicating that the correlation is not due to sequencing depth (Figure 6E). More variable transcription is also associated with higher amplitudes of rhythmic mRNA expression (Figure 6F). This correlation is valid for most AR-R genes (∼80%, Figure 6F), suggesting that this transcriptional variability or noise has a significant role in the emergence of rhythmic mRNA expression from arrhythmic transcription.

A comparison between the two replicates of rhythmic mRNAs indicates a better overlap than between one replicate and one replicate of rhythmic nascent RNAs. In contrast, the overlap between the two replicates of rhythmic nascent RNAs was no better than a single replicate rhythmic nascent RNA-rhythmic mRNA comparison (data not shown). Although the two comparisons cannot be definitive because of the limited six time point temporal resolution and resultant noise, they also support more pronounced variation at the transcriptional level than at the mRNA level.

Manual inspection of AR-R genes with high variability revealed a set of genes with high transcription at only one time point (10% of the AR-R gene set; 86/862). This was observed in both replicates and also correlates with rhythmic mRNA expression (Figure 6G). mRNA stability (half-life) regulation may contribute to the generation of rhythmic mRNA expression from what is likely a short burst of transcription. This ability of post-transcriptional regulation to generate rhythmic mRNA oscillations is selective, as not all arrhythmically transcribed genes with variable transcription exhibit rhythmic mRNA expression (Figure 6H). Moreover, not all genes exhibit variable transcriptional profiles (Figure 6B, see below).

About 20% of the AR-R genes are exceptions and exhibit higher variability at the mRNA than at the nascent level (Figure 6B,F). Because there is evidence that miRNAs can regulate mRNA levels independently of transcription, we examined whether those genes could be preferentially linked to miRNA regulation. Indeed, a higher number of predicted miRNA target sites was found for these genes compared to genes with higher transcriptional variability (using MirTarget2; p<0.01), suggesting that miRNAs contribute to mRNA cycling of genes with low transcriptional variability (Figure 6I).

The separation of these AR-R genes between genes with high or low transcriptional variability is therefore likely to be linked to different modes of gene expression regulation. Interestingly, this separation also reflects distinct biological functions, as AR-R genes with highly variable transcription were enriched for responsive genes (GO: response to stimulus) and metabolic genes (Figure 6J), indicating that their intrinsic transcriptional responsiveness is linked to their variable transcriptional profiles.

Genes within the R-R gene set include clock genes and many well-characterized clock controlled genes (see above). Because a large fraction of them are directly targeted by the core clock, we asked how CLK:BMAL1 regulates the transcription of its target genes at the genome-wide level. We also took advantage of our Nascent-Seq dataset to establish whether the phase differences between rhythmic CLK:BMAL1 DNA binding and rhythmic target gene mRNA expression (Rey et al., 2011) reflect transcriptional or post-transcriptional regulation.

To this end, we first performed a ChIP-Seq analysis of CLK and BMAL1 at a time of high DNA binding (ZT8). As expected, CLK and BMAL1 target many DNA binding sites in mouse liver (759 and 1579, respectively) and significantly overlap on 211 of these peaks (Figure 7A). Although highly significant (chi-square test, p<0.0001), the rather low fraction may indicate competition between CLK:BMAL1 and NPAS2:BMAL1 for binding sites. Importantly, about 90% of these 211 peaks have been previously characterized as rhythmic BMAL1 DNA binding sites in mouse liver (Rey et al., 2011).

CLK and BMAL1 sites overlap at their peak center, consistent with binding to DNA as a heterodimer, and they are enriched for the canonical consensus sequence CACGTG (Figure 7B,C) The data therefore indicate that almost all of the 184 direct target genes identified by the 211 CLK:BMAL1 DNA binding sites (Figure 6—source data 1) are bona fide direct target genes. They include the expected core clock genes (Figure 7E–G) as well as other interesting targets. There are for example 12 ncRNA genes, which include a cluster of four rhythmically transcribed ncRNAs (Figure 7H). These cycling ncRNAs suggest novel mechanisms by which CLK:BMAL1 impact circadian rhythms.

Although CLK:BMAL1 target genes are significantly enriched for rhythmically transcribed genes (chi-square test, p<0.0001; Figure 8A,C), there is a large discrepancy between the phases of rhythmic BMAL1 DNA binding and those of rhythmic transcription. This is because BMAL1 binding is essentially uniform at ZT3-5, whereas the transcription peaks are much more broadly distributed (Figure 8A,B). Remarkably, this is also true for the core clock genes. Whereas a small number have the expected phase similar to that of CLK:BMAL1 DNA binding (Rev-Erbα, Dbp; Figure 7E), the transcription of most target genes is significantly delayed. They include Per1 (Figure 7F), Cry1 (Figure 7G) and Per2 (8E). As recently proposed for Cry1 (Ukai-Tadenuma et al., 2011), this general delay in the peak of rhythmic transcription may be due to a collaboration of CLK:BMAL1 with other transcription factors (Rev-erbα and Dbp in the case of Cry 1).

To validate this interpretation, we assayed the pre-mRNA levels in Bmal1−/− mice of four clock genes that exhibit different phases of transcription in wild-type mice (Figure 8D). Only the levels of Rev-Erbα pre-mRNA in the mutant mice conform to expectation and are at trough levels of wild-type mice. The levels of Per1, Per2 and Cry1 pre-mRNAs in Bmal1−/− mice are higher than the trough of expression in wild-type mice, indicating that other transcriptional regulators are indeed important for the transcription of these three clock genes (Figure 8D). This notion is also supported by the enrichment of other transcription factor motifs such as HNF3/FOXA1, SP1 and E4BP4 adjacent to CLK:BMAL1 binding sites (Figure 7D).

Notably, Per2 has a prominent anti-sense transcript at times of low Per2 sense transcription (Figure 8E). Moreover, the 5′ end of this transcript coincides with peaks of Pol II (Figure 8E). This suggests that antisense transcription could be an additional mechanism responsible for the disconnect between the phase of CLK:BMAL1 DNA binding and the phase of rhythmic transcription.

3- to 6-month-old male mice housed in a 12 hr-light:12 hr-darkness (LD12:12) schedule were used. Wild-type mice (C57BL/6 strain) and Bmal1−/− mice (originally from Christopher A Bradfield; Bunger et al., 2000) were used. All experiments were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Brandeis Institutional Animal Care and Use Committee (IACUC protocol #0809-03).

Mice housed in LD12:12 were sacrificed every 4 hr (ZT0, 4, 8, 12, 16 and 20) by isoflurane anesthesia followed by decapitation. Mouse liver was then quickly removed and homogenized in 3.5 ml of 1× PBS and 3.5 ml of homogenization buffer (2.2 M sucrose, 10 mM Hepes pH 7.6, 15 mM KCl, 2 mM EDTA, 1× protease inhibitor cocktail [Roche, Basel, Switzerland], 0.15 mM spermine, 0.5 M spermidine, 0.5 mM DTT) with a dounce homogenizer (six strokes loose pestle, four strokes tight pestle). The liver homogenate was then mixed with 21.5 ml of homogenization solution and layered on the top of a 10 ml ice-cold cushion solution (2.05 M sucrose, 10 mM Hepes pH 7.6, 10% glycerol, 15 mM KCl, 2 mM EDTA, 1× protease inhibitor cocktail [Roche, Basel, Switzerland], 0.15 mM spermine, 0.5 M spermidine, 0.5 mM DTT) and centrifuged for 45 min at 2°C at 24,000 rpm (100,000×g) using a Bechmann SW27 rotor. Nuclei were resuspended in 1 ml of 20 mM Hepes pH 7.6, 150 mM NaCl, 2 mM EDTA, 1× protease inhibitor cocktail (Roche, Basel, Switzerland), 0.5 mM PMSF, 1 mM DTT, 0.5 U/ml of RNAseOUT/SUPERase-In (Invitrogen/Ambion), homogenized using a 1 ml dounce homogenizer (three times with loose pestle, two times with tight pestle), and divided into three samples of equal volume. One volume of 2× NUN buffer (50 mM Hepes pH 7.6, 2 M Urea, 2% NP-40, 600 mM NaCl, 2 mM DTT, 1× protease inhibitor cocktail [Roche, Basel, Switzerland], 0.5 mM PMSF, 0.5 U/ml of SUPERase-In [Ambion, Carlsbad, California]) was then added drop-by-drop while gently vortexing (level 2). Samples were left on ice for 20 min, then centrifuged at 24,000 rpm for 10 min at 4°C. The supernatant was removed and 1 ml of Trizol (Invitrogen, Carlsbad, California) was added to the DNA pellet. Samples were then incubated at 65°C for 15 min, and DNA pellet was then resuspended by gentle pipetting. Nascent RNA was then extracted using standard Trizol RNA extraction (Invitrogen, Carlsbad, California).

Mice housed in LD12:12 were sacrificed by isoflurane anesthesia followed by decapitation every 4 hr (ZT2, 6, 10, 14, 18 and 22). Mouse liver was then quickly removed and cut into small pieces that were frozen on dry ice. Total RNA was extracted using standard Trizol extraction (Invitrogen, Carlsbad, California).

Nascent RNA was first DNase-treated with TURBO DNase (Ambion, Carlsbad, California) using manufacturer's recommendation. Polyadenylated RNA was then removed from the nascent RNA using Dynabeads mRNA direct kit (Invitrogen, Carlsbad, California) following manufacturer's recommendations, and nascent RNA was precipitated (ethanol precipitation).

Sequencing libraries have been made using standard protocols. Briefly, 100 ng of purified nascent RNA were used to generate Illumina libraries. Nascent RNA was first fragmented using Fragmentation reagents (AM8740; Ambion, Carlsbad, California) by heating at 70°C for 5 min. Fragmented nascent RNA were then purified and used for standard Illumina library preparation. Following adaptor ligation, libraries of 200–300 bp length were size-selected on a 2% TAE agarose gel, and amplified by PCR for 15 cycles.

Strand-specific libraries were processed as above except for the following modifications that have also been described elsewhere (Levin et al., 2010). Briefly, after the first strand cDNA synthesis, dNTPs were removed by size-exclusion chromatography columns (G-50 columns; Amersham, Amersham, UK) and by ethanol precipitation using ammonium acetate. Second strand synthesis was then performed using a dNTP mixture containing dUTP instead of dTTP. After adaptor ligation and size selection (i.e., prior to the amplification), libraries were digested using Uracil-Specific Excision Reagent (USER Enzyme, NEB, Ipswich, Massachusetts) by incubating 2 units (2 μl) of USER Enzyme with 18 μl of libraries at 37°C for 30 min. Reaction was then heat-inactivated, the libraries were purified and PCR-amplified.

RNA-Seq libraries (e.g., mRNA) were made using Truseq RNA sample kit (Illumina, San Diego, California) following manufacturer's recommendations.

Adult mice housed in LD12:12 were sacrificed at ZT8 by isoflurane anesthesia followed by decapitation. Mouse liver was then quickly removed and homogenized in 3.5 ml of 1× PBS supplemented with 1% formaldehyde. After 10 min incubation at room temperature, cross-linking was stopped by mixing liver homogenate with 25 ml of ice-cold quenching solution (2.2 M sucrose, 150 mM glycine, 10 mM Hepes pH 7.6, 15 mM KCl, 2 mM EDTA, 1× protease inhibitor cocktail [Roche, Basel, Switzerland], 0.15 mM spermine, 0.5 M spermidine, 0.5 mM DTT). Homogenate was then layered on top of a 10 ml ice-cold sucrose cushion (2.05 M sucrose, 125 mM glycine, 10 mM Hepes pH 7.6, 10% glycerol, 15 mM KCl, 2 mM EDTA, 1× protease inhibitor cocktail [Roche], 0.15 mM spermine, 0.5 M spermidine, 0.5 mM DTT) and centrifuged for 30 min at 2°C and 24,000 rpm (100,000×g) using a Bechmann SW27 rotor. Nuclei were resuspended in 1 ml of 20 mM Hepes pH 7.6, 150 mM NaCl, 2 mM EDTA and sedimented at 1500×g for 1 min. Washed nuclei were resuspended in 1.2 ml sonication buffer (20 mM Hepes pH 7.6, 1% SDS, 150 mM NaCl, 2 mM EDTA) and sonicated on ice using a Fisherbrand Sonic Dismembranator at setting 2 (57 W) for 5 × 15 s to obtain chromatin fragments of about 500–1000 bp in length. The resulting chromatin was centrifuged at 15,000×g for 10 min and the resulting supernatant was aliquoted in 200 μl samples for immunoprecipitation and 25 μl samples for input.

Immunoprecipitation of chromatin was performed by mixing 200 μl of sonicated chromatin with 1.8 ml of IP buffer (10 mM Hepes pH 7.6, 150 mM NaCl, 2 mM EDTA, 0.1% NaDeoxycholate, 1% Triton X-100). Antibodies were added and samples were incubated overnight (Rabbit anti-BMAL1 antibody: 10 μl, ab3350; Abcam; Rabbit anti-CLK antibody: 10 μl, NB100-126; Novus Biologicals). Dynabeads protein G (100 μl per sample; Invitrogen, Carlsbad, California) were blocked in parallel overnight in 0.1 mg/ml yeast tRNA and 1 mg/ml BSA in IP buffer.

Following the overnight incubation, beads were washed once with IP buffer. The chromatin/antibodies mixture was then added to the beads and incubated at 4°C for an additional 2 hr. Beads were then washed once for 10 min with HSE I (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Hepes KOH, pH 7.6, 150 mM NaCl) and twice for 10 min with HSE II (0.1% SDS, 1% Triton X-100, 2 mM EDTA, 20 mM Hepes KOH, pH 7.6, 500 mM NaCl). Beads were then briefly washed with ice-cold TE and eluted with 200 µl of ChIP Elution buffer (50 mM Tris–HCl pH 8.0, 10 mM EDTA, 1% SDS, 1 mM DTT). 175 µl of ChIP Elution buffer was also added to the 25 µl input samples. Elution was performed at 65°C for 6–18 hr.

The resulting supernatant was removed, supplemented with 200 μl of TE and 8 μl of 1 mg/ml RNAse A (Ambion Cat #2271) and incubated at 37°C for 30 min. Then, 4 μl of 10 mg/ml proteinase K was added and samples were incubated at 55°C for another 2 hr. DNA was then isolated with PCR purification kit (Qiagen) and eluted with 40 μl of elution buffer.

BMAL1, CLK and input libraries were made from ChIPs performed from the same mouse liver extract. ChIP-seq libraries were made as described by Schmidt et al. (2009). Size-selected libraries of 200–300 bp length were used for Illumina deep-sequencing, whereas libraries with a 300–650 bp length were used for qPCR validation of the quality of the ChIP-seq libraries.

High-throughput sequencing has been performed as follow:

High-throughput sequencing has been performed following manufacturer recommendations and 8–12 pmol of libraries were hybridized to each lane of the flow-cells. Data were extracted and processed following Illumina recommendations. Sequences were aligned to the mouse genome (UCSC version mm9 database). Number of the sequences obtained for each library can be found in Table 1. Datasets are deposited on the Gene Expression Omnibus database under the accession number GSE36916↗ (GEO, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36916↗).

Sequences from the different libraries (fastq format) were first mapped to the mouse genome (version mm9) using bowtie (Langmead et al., 2009) with the command line: bowtie –q –a –-best –m 1. Only those that mapped uniquely to the mouse genome were used for further analysis, and their number has been used for normalization to compare signal difference between libraries.

ChIP-seq libraries were analyzed with the MACS algorithm (Zhang et al., 2008) by comparing the treatment sample (BMAL1 or CLK ChIP) to the control sample (Input) using the following criteria: effective genome size = 1.89 × 10⁹, tag size = 36, band width = 80, model fold = 5, p-value cutoff = 1 × 10⁻⁵. Significant peaks were computationally assigned to a gene. Briefly, a peak located between the transcription start site and the transcription start end of a gene was assigned to that gene, regardless of the ChIP-Seq peak position. The other peaks, referred as intergenic, were assigned to the gene with the closest transcription start site. Confirmation of this computational gene assignment was then confirmed by manual inspection for the 211 CLK:BMAL1 peaks. Visualization of the ChIP-seq signal was performed using the wig output file (from the MACS analysis) and the IGB browser.

Overlap between CLK and BMAL1 DNA binding peaks has been determined computationally using all significant peaks coordinates. Any overlap between the two peaks (even of one nucleotide) was considered significant. Quantification of the signal has been extracted from the raw data (number of reads per bp) normalized to sequencing depth of each library. For most experiments (e.g., Figure 7A,B), signal was binned using a 25 bp window.

Quantification of the number of e-boxes within BMAL1 and CLK DNA binding peaks has been performed computationally using peak fasta sequences. Enrichment has been calculated using the number of e-boxes found at a fixed position from the peak center divided by the expected number of e-boxes. The maximal window size (difference from the fixed position and the peak center) was 500 nt, as the number of expected e-boxes dropped to the background at this window size.

Motif analysis has been performed using MEME suite (http://meme.sdsc.edu/meme/intro.html↗), using a 100 bp sequence for each peak (peak center ± 50 bp). Parameters were as follow: -dna -mod anr -nmotifs 20 -minw 6 -maxw 30. The background model contained the same nucleotide distribution as the input file. Significant motifs were then analyzed using TOMTOM from MEME suite.

Gene signal (reads per base pair) was averaged for the 12 time points, and the ratio Nascent RNA/mRNA was calculated. Genes with a ratio over 2 SD from the average of all ratios were selected for gene ontology (GO) analysis; 302 genes had a Nascent-Seq/RNA-Seq ratio below 2 SD and 463 genes had a ratio over 2 SD. 13595 genes were considered for this analysis. GO analysis has been performed using GOToolBox (Martin et al., 2004) (http://genome.crg.es/GOToolBox/↗), using an hypergeometric test with Benjamini-Hochberg correction. The background model consisted of the entire list of genes.

Only genes with more than three reads per base pair for at least one time point of the Nascent-Seq dataset and two read per base pair for the RNA-Seq dataset were further considered for subsequent analysis (see above). Rhythmically expressed genes were determined based on three parameters: amplitude, F24 (24-hr spectral power, see below) and p-value of the F24.

The amplitude was calculated by dividing the highest value of the 12 time points by the lowest value. The F24 score was calculated by Fourier transformation using a R code originally described by Wijnen et al. (2005). Briefly, normalized reads per base pair from the two independent six-time points rhythms were concatenated and the 24-hr spectral power (F24) was determined for each gene. The F24 score (expressed in range 0–1) indicates the relative strength of the extracted rhythmic component.

The F24 p-value (pF24) represents the probability of observing an F24 score from randomly permuted data that is of equal or greater strength than the extracted Fourier component. It was calculated for each gene after performing 10,000 randomized permutations of the rpbp values. A pF24 was considered significant if (pF24 < 0.05) if the experimental pF24 was within the top 5% of the 10,000 pF24 calculated from randomized permutation.

Transcripts were considered to be rhythmically expressed when meeting the three following criteria: (1) pF24-values < 0.05 (i.e., experimental pF24-value is within the top 5% of the 10,000 pF24-values calculated from randomized permutation), (2) F24 > 0.45 and (3) amplitude (maximal/minimal experimental values) > 1.5. A more stringent cut-off was also used to identify strong rhythmic expression: pF24 < 0.05, F24 > 0.6 and amplitude > 1.75.

The phase information from the Fourier transformation (which indicates the peak of the cosine curve) was further used to assess phase difference between rhythmic nascent RNA and mRNA expression (Figures 3C,4B).

Total RNA from wild-type and Bmal1−/− mice was prepared from mouse liver using Trizol reagent (Invitrogen) and DNAse-treated using Turbo DNAse (Ambion) according to the manufacturer's protocols. cDNA derived from RNA (using Invitrogen Superscript II and random primers) was utilized as a template for quantitative real-time PCR performed with the Rotor-Gene 3000 real-time cycler (Qiagen). The PCR mixture contained Platinum Taq polymerase (Invitrogen), optimized concentrations of Sybr-green (Invitrogen) and specific primers for either Rev-Erbα, Per1, Per2 or Cry1 pre-mRNAs. Quantitative PCR using Actg1-specific primers was used as a loading control. Cycling parameters were 95°C for 3 min, followed by 40 cycles of 95°C for 30 s, 55°C for 45 s, and 72°C for 45 s. Fluorescence intensities were plotted vs the number of cycles by using an algorithm provided by the manufacturer. mRNA levels were quantified using a calibration curve based upon dilution of concentrated cDNA.