Oscillating and stable genome topologies underlie hepatic physiological rhythms during the circadian cycle

To elucidate the dynamics of chromatin architecture along the day-night cycle in mouse tissues, we performed 4C-seq experiments every 4h for 24h as in [23] (Material and Methods). 4C-seq probes the interaction frequencies between one “bait” DNA fragment and the entire genome [25]. We placed 4C-seq baits at gene promoters of central components of the mammalian molecular clock such as Bmal1 [3]. Bmal1 gene is rhythmically transcribed in WT mouse liver with pre-mRNA abundance peaking around ZT22 (Fig 1A) (ZT: Zeitgeber time; ZT0 corresponds to onset of lights-on; ZT12 corresponds to onset of lights-off). Note that for circadian clock-driven gene expression, the pre-mRNA accumulation is an appropriate proxy for transcription, which typically peaks several hours before the mRNA [9]. As reported for other genes [23], the 4C-seq contacts for Bmal1 were highly enriched on the cis chromosome, especially within a 2Mb region surrounding the bait position (S1A Fig). This region contained ~50% of cis counts for all time points (S1 Table) and comprised most gene regulatory interactions [16]. 4C-seq data were then normalized and analyzed applying a locally weighted multilinear regression (LWMR) using a Gaussian window (sigma=2.5 kb) centered on each fragment for local smoothing [23] (Material and Methods). Temporal analysis revealed that the Bmal1 promoter rhythmically contacted a genomic region spanning from ~40 kilobases (kb) to ~75 kb downstream of the transcription start site (TSS), with the contact frequency peaking around ZT18-20 at multiple 4C-seq peaks (Figs 1B, 1C and S1A). To characterize interacting regions, we integrated time-resolved chromatin immuno-precipitation followed by high-throughput sequencing (ChIP-seq) experiments targeting PolII and chromatin marks typical of enhancer regulatory elements such as H3K27ac and DNase1 hypersensitive sites (DHS) [8]. As expected, PolII loading was rhythmic across the entire Bmal1 gene body and peaked around ZT18-22 (Fig 1D). Multiple 4C-seq interaction sites peaking at ZT20 coincided with chromatin regions marked by rhythmic regulatory activity. For example, a preferential ZT20 contact ~45 kb downstream of the bait corresponded to a conserved intronic region marked by rhythmic H3K27ac histone acetylation and DNase1 hypersensitivity peaking around ZT20 (Fig 1D, blue region in log₂ fold-change). Another ZT20 4C-seq peak located near exon5 ~73 kb downstream of the bait coincided with a ZT20 DNase1 signal (although weak) (Fig 1C and 1D). These data suggested oscillating DNA loops between rhythmically active enhancers and the Bmal1 promoter.

Next, we explored chromatin architecture dynamics surrounding Period1 (Per1) and Period2 (Per2) genes that belong to the negative limb of the circadian molecular oscillator [3]. Per2 and Per1 4C-seq signals were largely enriched within a 2Mb window on cis chromosome (S1 Table and S1B and S1C Fig). Per2 pre-mRNA is highest around ZT14 in WT liver (Fig 2A). At the Per2 locus, a large region extending from ~35 kb to ~70 kb downstream of the TSS contacted more frequently the promoter at ZT16, with two prominent oscillating contacts at 40 kb and 65 kb downstream of the TSS (Fig 2B and 2C). The region 40 kb downstream of Per2 TSS corresponded to multiple intragenic sites near the 3’ end of Per2 containing rhythmic transcription and enhancer chromatin marks (PolII, H3K27ac, DHS) peaking around ZT16, as well as the Hes6 gene in which PolII, H3K27ac and DHS peaked at ZT16 (Fig 2D). The region 65 kb downstream of the TSS also contained H3K27ac and DHSs enhancer marks (Fig 2D). Thus, these data showed 24h rhythms in enhancer-promoter contacts accompanying Per2 transcription, and also dynamic gene-gene interactions with synchronised transcription. Furthermore, the Per1 pre-mRNA level is maximal near ZT10 in WT mouse liver (S2A Fig). Overall, the chromatin rhythms showed lower amplitudes for Per1 compared to Per2. While several sites within the 2Mb genomic region surrounding Per1 showed weakly rhythmic contacts with the Per1 promoter (S1C Fig), two proximal regions contacted the Per1 TSS rhythmically (S2B and S2C Fig). The first region was immediately upstream of the Per1 TSS and corresponded to multiple sites with PolII, H3K27ac and DHSs (S2D Fig). The temporal profile of this interaction showed highest contact at ZT8 (S2C Fig), but the phase of the harmonic fit was later near ZT12. The second rhythmic site located ~15-20 kb downstream of the TSS showed a peak phase near ZT14 and was in fact located in a valley of 4C signal. This site corresponded to the 3’ terminal region of Per1 and coincided with H3K27ac and DHS sites, as well as the promoter region of Hes7 (S2D Fig). These data suggest that sites resembling enhancer elements diurnally contacted the promoter of Per1. We also noted a weak ZT02 preferential interaction 45 kb upstream the Per1 TSS (S2B Fig). Together, the measured chromatin interaction patterns at Per1 are more complex to interpret compared to Per2, possibly because of weaker signals and limitations in the phase estimations.

Together, these time-resolved 4C-seq experiments revealed oscillating contacts between core-clock gene promoter and surrounding enhancers. The temporal dynamics of many of those DNA interactions as well as chromatin features at enhancer sites were synchronized with the transcription of the target gene, with high contact frequency and enhancer activity coinciding with the peak time of pre-mRNA synthesis.

Next, we explored chromatin interactions surrounding clock output genes. In mouse liver, two main transcriptional waves centered around ZT08 and ZT20 [26] underlie daily rhythms in physiology, including detoxification [27], glucose [28] and lipid metabolism [29], and metabolite synthesis [30–32]. Notably, clock-related TFs bind to the promoter of genes involved in carbohydrate and lipid metabolism [6,33]. Thus, we selected the promoter of genes that were rhythmically transcribed (assessed by pre-mRNA peak times) with peak times that coincided with either of the two main transcriptional waves, and we profiled chromatin interactions at the two time points (ZT08 and ZT20) [26]. In particular, we placed a bait at the promoter of Mreg that is rhythmically transcribed specifically in the liver [5]. Mreg pre-mRNA peaks near ZT22 in WT mouse liver (Fig 3A). Mreg 4C-seq signals were enriched within the 2Mb region surrounding the bait position (S1 Table). In this signal-rich region, the 4C interaction profiles showed no difference between ZT08 and ZT20 (Fig 3B, Z-scores were centered around zero, Material and Methods) with the exception of a bait-proximal region that showed a clear preferential contact at ZT20. This region spanned from ~10 kb to ~50 kb downstream of the bait and corresponded to Mreg intragenic region. The ZT20 preferential contact was the highest at ~30 kb to ~40 kb downstream of the TSS (Fig 3D). Remarkably, while PolII was rhythmic and peaked around ZT20 across the entire gene body, H3K27ac and DNase1 hypersensitivity signals were high specifically within the region spanning from ~30 kb to 40 kb downstream the TSS, and both marks peaked at ZT20 (Fig 3C and 3D). This data indicated that the region recruited to the Mreg promoter preferentially at ZT20 corresponded to DNA elements having rhythms in enhancer chromatin signature peaking at ZT20 in WT mouse liver. Thus, oscillating DNA loops can connect daily active enhancers with promoters of clock output genes.

Next we explored the dynamics of chromatin conformation surrounding genes involved in key physiological function in liver such as Nampt, that encodes the nicotinamide phosphoribosyltransferase rate-limiting enzyme in the NAD biosynthesis pathway [30,31]. BMAL1 binds to the promoter of Nampt [6] and, the NAD+-dependent histone deacetylase SIRT1 inhibits CLOCK-BMAL1 TF activity [31], illustrating the interlocking between the clockwork machinery and the metabolic state of the cell. Nampt pre-mRNA accumulates rhythmically in the liver of WT mice and peak at ZT10 (Fig 4A). Again, Nampt 4C-seq signals were mostly confined within the 2Mb region surrounding the bait at both time points (S1 Table). Unlike the dynamics observed for core clock genes and Mreg, no differential interactions between ZT08 and ZT20 were found across the entire 2Mb signal-rich region (Fig 4B). In particular, major 4C-seq interaction sites stood out at 50 kb and 125 kb upstream of the bait position (Fig 4B), suggesting that these regions contacted the promoter of Nampt at a similar frequency during day and night. Time-resolved ChIP-seq experiments showed rhythmic loading of PolII along the Nampt gene body peaking at ZT10, consistently with the rhythmic accumulation of Nampt pre-mRNAs (Fig 4A, 4C and 4D). The two upstream interacting regions were marked by rhythmic PolII loading (although the ChIP-seq signal is weak) and H3K27ac and DHS signals peaking at ZT10 (Fig 4C and 4D). We also noted secondary 4C-seq peaks 200 kb upstream and 60 kb downstream of the Nampt bait at sites marked by DHS and weak H3K27ac. These data showed DNA loops connecting the promoter of Nampt with rhythmically active enhancers, at contact frequencies that were similar at ZT08 and ZT20. The times corresponded, respectively, to the peak and trough in transcription of Nampt in WT mouse liver.

Consistently with the observations at the Nampt locus, we identified additional cases of rhythmically active enhancer-promoter pairs forming DNA loops that were stable during active and inactive transcription at the two clock output genes Pfkfb3 (S3 Fig) and Mfsd2a (S4 Fig). As shown in other model systems [34,35], this suggested that the clock-controlled transcriptional machinery can act over a frozen promoter-enhancer contact network that is insensitive to transcriptional activity.

Next, we asked if the stable chromatin topology surrounding the rhythmically transcribed Nampt gene promoter was maintained in clock-impaired animals. Therefore, we profiled chromatin conformation at the Nampt locus in livers of Bmal1 KO animals at ZT08 and ZT20. At the transcriptional level, Nampt shows lower and constant levels in livers of those arrhythmic animals compared to WT (S5A Fig). Overall, the distributions of 4C-seq signals were comparable between time points and genotypes (S1 Table). Remarkably, the Nampt promoter-enhancer loops were maintained at similar levels in clock-impaired mice compared to WT and constantly from ZT08 to ZT20 (S5B, S5D and S6A Figs). Furthermore, PolII loading and H3K27ac chromatin marks were overall arrhythmic and lower at connected promoter-enhancer regions in Bmal1 KO compare to WT, consistent with pre-mRNA profiles (S5 Fig). These data suggested that despite altered transcriptional output, the static promoter-enhancer loops detected in wild-type were unaltered in Bmal1 KO at the Nampt locus.

The above examples suggested that changes in transcription activity states could occur over a largely static conformation of chromatin. As shown in other systems [36,37], in such a model, gene-gene interactions might allow transcriptional co-regulation. A remarkable example arguing in favor of this model was observed at the locus of the liver-specific and rhythmically expressed gene Por [5]. Por encodes the cytochrome P450 oxidoreductase enzyme involved in the NADPH-dependent electron transport pathway; its rhythmic expression along the diurnal cycle contributes to detoxification in the mouse liver [27,38]. In WT livers, ZT08 and ZT20 4C-seq experiments using the Por promoter as a bait revealed chromatin contacts with a region spanning from ~50 kb upstream to ~ 100 kb downstream of the bait position, showing multiple local peaks (Fig 5B). The interaction frequency was similar between time points, as shown by the zero centered Z-scores in the cis region (Fig 5B). The 50 kb region upstream corresponded to Rhbdd2 gene, while the downstream interacting region coincided with the Tmem120a, Mir7034l, Styxl1 and Mdh2 genes (Fig 5D). PolII loading peaked at ZT10 across the entire interacting locus, coinciding with a phase-coherent accumulation of the pre-mRNAs around ZT10 (Fig 5A and 5D), indicating that the ZT08 and ZT20 time points measured chromatin contacts during the active and inactive transcription phases, respectively. H3K27ac signals marked all prominent interacting sites, showing rhythms at the Por and Rhbdd2 genes, peaking around ZT10. These data suggest that genes with synchronous transcription cycles contacted each other in WT liver, and at a constant frequency during the peak and trough of their transcription.

To further investigate the model of interaction between phase-coherent genes, we performed 4C-seq at the Por locus in livers of Bmal1 KO mice. In these animals, the Por interacting sites were also connected at similar frequencies at ZT08 and ZT20, and at a comparable level to the WT conditions (Figs 5C and S6B), suggesting that the chromatin architecture at this locus remained stable and did not depend on a functional clock. Furthermore, Por, Rhbdd2, Tmem120a, Styxl1 and Mdh2 transcripts coherently accumulated at overall dampened, time delayed, and lower levels in Bmal1 KO compared to WT mice (Fig 5A). Consistently, PolII loading and H3K27ac chromatin marks no longer oscillated across the entire interacting locus in livers of Bmal1 KO compared to WT, and the levels of histone acetylation were also reduced in clock-deficient animals (Fig 5D). These data demonstrated that the Por interacting chromatin hub connected genes sharing similar temporal dynamics of transcription across the circadian cycle, and that the chromatin hub structure was insensitive to transcriptional changes and clock-independent.

To investigate the transcriptional regulatory function of DNA loops, we analyzed ChIP-seq data of core-clock transcription factors, in particular BMAL1, REVERB-alpha and ROR-gamma. Indeed, it was suggested that, in mouse liver, BMAL1 could connect distal regulatory elements [23,24], while REVERB-alpha, through the recruitment of co-factors, disrupts chromatin loops [22]. Furthermore, we also considered CTCF ChIP-seq data in mouse liver [39], since its role in chromatin looping is well characterized [40]. The region rhythmically recruited to the promoter of Bmal1 (spanning from ~40 kb to ~75 kb downstream of the TSS, Fig 1B and 1D) coincided with multiple localized peaks in DNase1 hypersensitivity that were synchronous with H3K27ac rhythms, and bound by REVERB-alpha, ROR-gamma and CTCF (S7A Fig). At the Per2 locus, we similarly observed binding of REVERB-alpha, ROR-gamma and BMAL1 and CTCF signal in the region spanning from ~35 kb to ~70 kb downstream of the TSS and containing synchronous rhythms in DNase1 hypersensitivity and H3K27ac (S7B Fig). At the Per1 locus, BMAL1, REVERB-alpha, ROR-gamma and CTCF bound the regions rhythmically contacting the promoter (S8A Fig). These data suggested that the binding of core-clock transcription factors and CTCF at regions rhythmically recruited to the promoter of core-clock genes participate in the temporal dynamics of DNA contacts and transcription regulation. Furthermore, the rhythmic promoter-enhancer loops at the clock output gene Mreg were also bound by core-clock TF and CTCF (S8B Fig). At stable loops, we also observed the binding core-clock TFs and CTCF at connected sites (S9 and S10 Figs). For example, regions contacting the Nampt promoter were bound by REVERB-alpha and ROR-gamma (S9A Fig). At the Por and Pfkfb3 loci, for which large surrounding regions connected the respective promoters, the multiple localized 4C-seq peaks coincided with ChIP-seq sites for core-clock TFs and CTCF (S9B and S10A Figs). Together, these data support a model in which binding of specific core-clock TFs as well as chromatin architecture factors at distal region connecting clock output target gene promoters participate in 24h rhythmic transcription regulation.

All experiments were approved by the Ethical Committee of the State of Vaud Veterinary Office (authorization VD3109).

C57/BL6J (WT) and Bmal1 KO animals were maintained at the EPFL animal house facility in 12hour/12hour light-dark cycle with 4 animals per cage.

4C-seq was performed in livers of 8 to 12 weeks old male with 4 biological replicates when comparing ZT08 versus ZT20 in WT animals and 3 biological replicates in Bmal1 KO. 3 biological replicates were used in the 4h time-resolved 4C-seq experiments. Sample preparation and analyses were performed as in [23]. In brief, livers were isolated and perfused with PBS before homogenization in 4 mL of 1×PBS including 1.5% formaldehyde for 10 minutes at room temperature. 25 mL of the following ice-cold buffer (2.2 M sucrose, 150 mM glycine, 10 mM HEPES at pH 7.6, 15 mM KCl, 2 mM EDTA, 0.15 mM spermine, 0.5 mM spermidine, 0.5 mM DTT, 0.5 mM PMSF) was added to the homogenates and kept for 5 min on ice. Homogenates were loaded on top of 10 mL of cushion buffer (2.05 M sucrose, 10% glycerol, 125 mM glycine, 10 mM HEPES at pH 7.6,15 mM KCl, 2 mM EDTA, 0.15 mM spermine, 0.5 mM spermidine, 0.5 mM DTT, 0.5 mM PMSF) and centrifuged at 10⁵ g for 45 minutes at 4°C. Nuclei were washed twice in 1× PBS and resuspended in 1 mL of 10 mM Tris-HCL (pH 8.0), 10 mM NaCl, 0.2% NP-40, and protease inhibitor cocktail (Complete Mini EDTA-free protease inhibitor cocktail; Sigma-Aldrich); kept for 15 minutes on ice; and washed twice with DpnII buffer (New England Biolabs). Thirty million nuclei were incubated for 10 minutes at 65°C in DpnII buffer and triton X-100 was added to 1% final concentration. Chromatin was digested overnight with 400 U of DpnII (New England Biolabs) at 37°C with shaking. Digested chromatin was then diluted in 8-mL of ligation buffer containing 3000 U of T4 DNA ligase for 4 h at 16°C plus 1 h at room temperature. 50 μL of 10 mg/mL proteinase K was added and samples were incubated overnight at 65°C. DNA was precipitated and resuspended in TE buffer (pH 8.0) containing RNase A, and incubated for 30 min at 37°C. Libraries were digested with NlaIII using 1U/μg of DNA template (New England Biolabs) overnight at 37°C. Digested products were ligated with 2000 U of T4 DNA ligase (New England Biolabs) for 4 h at 16°C in a 14-mL final volume. Circularized products were precipitated and resuspended in TE buffer (pH 8.0). Inverse PCRs were performed on 600 ng of circular DNA template per sample as described in (23). Inverse PCR primers are mentioned in the S2 Table.

4C-seq data were analyzed as in [23]. Briefly, demultiplexed read counts were mapped to the mouse genome (mm9) using HTSstation [46]. Samples were excluded from the analysis if more than 75% of restriction fragments did not have any count on a 2Mb region surrounding bait fragment (S1 Table). The first five NlaIII fragments upstream and downstream of the bait were excluded in the analysis since they were suspected to be partially digested or self-ligated products. The 4C-seq signal was calculated using a locally weighted multilinear regression model [23]. Fragment counts for each sample were normalized by the total fragments on the cis-chromosome (excluding the five fragments upstream and downstream of the bait). To stabilize variance, the fragment counts c in each sample were log-transformed:

with p=500. For each position, the 4C-seq signals (Y) were modeled with fragment effects a_i and condition effects b_j (which can be time, tissue, or genotype). We estimated these effects by minimizing the weighted sum S of squared residuals across replicates r:

S=argmina,b∑i,j,rWi,j(Yi,j,r−ai−bj)2, with weights W_i,j are defined as W_i,j = w_g,i×w_s,j, where w_g,i is the Gaussian smoothing kernel (sigma=2500bp) at position i, and w_s,j a condition weight based on the number of samples with non-zero counts on fragment i. To compare between two conditions, we calculated p-values and condition effects at each genomic fragment using t-statistics. To detect rhythmic signal, we calculated the 24-hour Fourier coefficients of the condition effect (real and imaginary parts) from the six equally-spaced time points and used the chi-square test to test deviations from the null model that the real and imaginary parts have both a mean of zero.

Processed RNA-seq data were downloaded from [47] (GSE73554↗) and rhythms were analyzed as in [5].

Data were downloaded from [8] (GSE60430↗) and analyzed as in [23]. We binned the ChIP-seq and DNase1-seq signal (log2 counts per million) into 500 bp windows. We smoothed the signal by taking a running average across 7 bins (3 bins upstream and 3 bins downstream of the current bin). For each bin, we calculated the amplitude and phase by fitting a harmonic regression model with a 24-hour period across the 7 bins. The rhythmic signal (amplitude and phase) was mapped to a color using hue (time of maximum signal), saturation (set to 1), and value (increased with increasing statistical significance) color scheme.

To obtain smooth color transitions, the value v was calculated using a Hill function with Hill coefficient n = 5 and v=mini∈(a,p)(−log(xi)5ki5−log(xi)5), with k_a = 0.5, k_p = 4.5 and x_a, x_p being amplitude and –log₁₀(p) of the harmonic regression fit.

BMAL1 ChIP-seq data were downloaded from [6] (GS26602↗). REVERB-alpha and ROR-gamma ChIP-seq data were downloaded from [48] (GSE67973↗). CTCF ChIP-seq data were downloaded from [39].

All data have deposited on the Gene Expression Omnibus with the accession number GSE139195, and are fully accessible on the following link. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE139195↗ Analysis code and processed data can be explored through the Circadian4Cseq R package, which is available on github (https://github.com/jakeyeung/Circadian4Cseq↗).