Rhythmic transcription of Bmal1 stabilizes the circadian timekeeping system in mammals

Bmal1 gene is essential for circadian oscillation^8,27 (Supplementary Fig. 1), and two RREs in the Bmal1 5’-UTR region are highly conserved among mammals (Supplementary Fig. 2a). Luciferase reporters containing the RREs has been widely used for monitoring cellular bioluminescence rhythms^20,24. Deletion or mutation of these two RRE sequences in a luciferase reporter abolished the bioluminescence rhythm (Supplementary Fig. 2b–e) as was shown previously^20,24,28, indicating that the two RREs are important for generating rhythms of the reporter. To examine contribution of the two RREs to the transcriptional rhythm of endogenous Bmal1 gene, we generated mutant NIH3T3 cells and mutant mice that are deleted of the two RREs by using CRISPR-Cas9 system. Specifically, we isolated three independent lines of mutant cells and two lines of mice (named ΔRRE mutants), all of which had a double allelic deletion in the RREs (Fig. 1b and Supplementary Fig. 3). The ΔRRE mutant mice were born in the expected Mendelian ratio with normal morphology. In quantitative RT-PCR analysis, we observed clear Bmal1 expression rhythm in wild-type (WT) NIH3T3 cells, and the Bmal1 rhythm was abrogated in the ΔRRE mutant cells (Fig. 1c). Similarly, the Bmal1 expression rhythm was abolished in the mutant mouse liver, kidney, heart, and epididymal white adipose tissue (eWAT) (Fig. 1c). These results demonstrated that transcriptional regulation through the two RREs located in the Bmal1 5ʹ-UTR region is essential for the circadian expression of endogenous Bmal1 gene.

The amounts of Bmal1 mRNA in the ΔRRE mutant NIH3T3 cells were maintained at a high level that was comparable to the peak of the Bmal1 mRNA rhythm in WT, whereas the mutant mouse tissues showed constant expression levels that were in the middle or lower range (Fig. 1c). These data suggest that the balance in transcription factors between the activators (RORs) and the repressors (REV-ERBs) of the RRE-mediated transcription is different among cells and tissue types. This idea is supported by the fact that ROR isoforms are expressed at quite different levels among tissues, whereas REV-ERBs show relatively ubiquitous expression in almost all tissues^25,27.

Even in the absence of the rhythmic Bmal1 transcription, we found that BMAL1 phosphorylation rhythm was still maintained (Fig. 3c). In particular, the ΔRRE mutant and WT models simulate nearly identical rhythms of the BMAL1 phosphorylation. To examine this prediction, we performed immunoblot analysis of BMAL1 protein phosphorylation, which can be monitored by its band-shift in SDS-PAGE gels²³. Then we observed time-of-day-dependent band-shifts of the BMAL1 protein bands in the ΔRRE mutant cell lysates and the mutant liver nuclei, with temporal profiles that were quite similar to those observed in WT (Fig. 3d). Taken together, functional rhythm of BMAL1 protein represented by its phosphorylation rhythm^21–23,32 was maintained even in the absence of the Bmal1 mRNA rhythm, and this should contribute to the apparently normal rhythmicity of the other clock genes in the ΔRRE mutant (Fig. 2c–h).

We investigated the molecular mechanism underlying the phosphorylation dynamics (functional rhythm) of BMAL1 in the ΔRRE mutant by using the mathematical model. We first examined the possibility that rhythmic transcription of Clock gene through its RRE³ may drive the functional rhythm of BMAL1 in CLOCK-BMAL1 heterodimer, because the heterodimerization is known to regulate the BMAL1 phosphorylation^22,23,33. In the model, we eliminated the RRE-mediated transcriptional regulation of Clock gene in addition to that of Bmal1 gene. Even under the constant expression of Clock and Bmal1 mRNAs, the BMAL1 phosphorylation rhythm was unaffected (Supplementary Fig. 8a–c), suggesting that the Clock expression rhythm plays a minimal role for driving functional rhythm of BMAL1 in the ΔRRE mutant. Experimentally, shRNA-mediated Clock knock-down reduced the amplitude of the cellular rhythms in both the ΔRRE mutant and WT cells, and the degree of the amplitude reduction was similar to each other (Supplementary Fig. 8d, e). It was also reported that Clock gene is almost constantly expressed in the mouse SCN²⁰. Together, these observations support our mathematical model predicting that the expression rhythm of Clock gene is dispensable for driving the BMAL1 functional rhythm.

This idea led us to investigate whether the ΔRRE mutant and WT differently respond to perturbation of Cry transcription. In the mathematical model, we increased the transcriptional activity of REV-ERBs that is known to suppress transcription of Cry1 gene through its RRE elements⁶. This completely abolished the rhythmicity of CRY1 protein and phosphorylated BMAL1 levels in the ΔRRE mutant model but not in the WT model (Fig. 4d). Experimentally, NIH3T3 cells were treated with SR9009, an agonist of REV-ERBs³⁵, and indeed, we observed that the amplitude of the bioluminescence rhythm was more severely reduced in the ΔRRE mutant cells than in the WT cells (Fig. 4e, f). Furthermore, the SR9009 treatment lengthened the circadian period, and the effect was significantly stronger in the ΔRRE mutant cells (Fig. 4g and Supplementary Fig. 10c). These observations demonstrated that the ΔRRE mutation increased susceptibility of the circadian oscillation to perturbations of CRY1 protein.

Taken together, we focused on the BMAL1 phosphorylation rhythm in the absence of the Bmal1 transcriptional rhythm, and found that the CRY1 protein rhythm is a key contributor to maintaining the apparently normal circadian oscillation in the ΔRRE mutant cells. The RRE-mediated transcriptional rhythms of the two clock genes Bmal1 and Cry1 cooperate to generate the functional rhythm of BMAL1 protein and hence leads to the robust circadian rhythm of E-box-mediated transcription (Fig. 4h).

To further explore the advantage of the interlocking structure composed of the two feedback loops, we examined the effects of the ΔRRE mutation on global changes of transcriptome (Supplementary Fig. 13, Supplementary Data 1), proteome (Supplementary Fig. 14, Supplementary Data 2, 3) and phospho-proteome (Supplementary Fig. 15, Supplementary Data 4, 5). We collected livers from the ΔRRE mutant and WT mice at 6 time points (4-hr intervals) throughout the day (Fig. 5b). In the liver transcriptome analysis, we found that the number of rhythmic genes was reduced in the mutant (Fig. 5c). Such reduction was observed particularly in a population of genes peaked at CT6-CT12 (Fig. 5d, Supplementary Fig. 13). Because the circadian expressions of E-box-regulated genes are known to peak at around CT6-CT12 in the mouse liver³, we next compared the rhythmicities of the E-box-regulated genes between the ΔRRE and WT livers. Among the E-box-targeted 2186 genes identified in our previous CLOCK-ChIP-Seq study³⁶, 1717 genes were detected as expressed in the present RNA-seq analysis. The heat map of the 1,717 genes showed that the temporal expression profiles of the E-box-targeted genes were markedly affected by the ΔRRE mutation (Fig. 5e). The circadian expressions of the rhythmic genes peaked at CT6-CT12 (including many E-box-regulated genes) might be perturbed by the impaired E-box-mediated transcriptional rhythms in the ΔRRE mutants. This observation raised the possibility that the RRE-mediated feedback loop through the Bmal1 rhythm affected the circadian rhythm of E-box-mediated transcription through the functional rhythm of BMAL1. In the simulation, the ΔRRE mutation caused a lengthened phase relationship between peak times of phosphorylated BMAL1 rhythm and mRNA rhythm of Rev-erbα, an E-box-regulated gene (Supplementary Fig. 16). This altered relationship would be due to the slower decrease in abundance of active BMAL1 protein (not sequestered by the repressors) in the ΔRRE mutant model. Consequently, the E-box-mediated transactivation driven by active BMAL1 may be delayed in the ΔRRE mutant, leading to lengthening of the circadian period observed in the cellular bioluminescence rhythm (Fig. 2f, g and Supplementary Fig. 6a, b) and the mathematical model (Fig. 3b).

Collectively, the mathematical modeling and experiments on the ΔRRE mutant demonstrated that the RRE-mediated feedback loop plays pivotal roles in stabilizing the circadian timekeeping system and organizing its appropriate output. The present study suggests that the structure consisting of the two interlocking feedback loops provide evolutionary benefit to the circadian clocks.

NIH3T3 cells (RIKEN Cell Bank) were maintained at 37 °C under 5% CO₂, 95% air in Dulbecco’s modified Eagle’s medium (DMEM) (Sigma) supplemented with 10% fetal bovine serum (FBS) and Pen Strep (Gibco) containing final 100 units/mL penicillin and final 100 µg/mL streptomycin.

A pX330 vector containing SpCas9 and a chimeric guide RNA (Addgene, cat. no. 42230)^45,46 was digested by BbsI. For deleting the RRE elements located in the Bmal1 5’-UTR region, two sgRNA sequences (Table S1) were designed using an online software, CRISPR direct⁴⁷ (Supplementary Fig. 3). The oligonucleotides containing the sgRNA sequences were synthesized (Sigma), phosphorylated, annealed, and ligated with the pX330 vector. NIH3T3 cells were transfected with the plasmids by using Lipofectamine3000 Reagent (Thermo Fisher scientific). After the transfection, cells were seeded at a density of one cell per well in 96 well plates. Wells containing a single colony derived from a single cell were selected. The genome DNA sequences of the target site were checked by direct Sanger sequencing.

We performed real-time monitoring of luciferase expression as described previously⁴⁸ with minor modifications. Briefly, NIH3T3 cells plated on 35-mm dishes were transiently transfected with Bmal1-luc/pGL4.12. After 24 hrs, the cells were treated with 0.1 µM (final) dexamethasone (Sigma) for 2 hrs to synchronize circadian phase of the cellular clocks in the dish. Then the media were replaced by recording media: phenol red-free DMEM (Sigma) supplemented with 10% FBS, 3.5 g/L glucose, Pen Strep (Gibco) containing 100 units/mL penicillin and 100 µg/mL streptomycin, 0.1 mM luciferin, and 10 mM 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (Hepes)-NaOH (pH 7.0) at the final concentrations. The time point of this medium change was defined as Time After Dex (TAD) 0. The bioluminescence signals were continuously recorded for 6 to 8 days in air with Dish Type Luminescencer, Kronos Dio (ATTO) or LumiCycle (Actimetrics) at 37 °C, unless otherwise mentioned.

For knockdown of Clock, NIH3T3 cells were transiently transfected with 1 ng or 3 ng of Clock shRNA in combination with Bmal1-luc/pGL4.12. Target sequence 5’-GAACA TCAGG CTATG ATTAC T-3’ was inserted into pSilencer 3.1-H1 puro vector (Thermo Fisher scientific), and its empty vector was used as a negative control.

To perturb CRYs protein rhythms, SR9009 (Cayman Chemical Company) or KL001 (Cayman Chemical Company) was dissolved in DMSO and mixed with the recording media. For the control, the same volume of DMSO was added as vehicle to the recording media.

Raw data of the bioluminescence rhythms were smoothed by 2-hr moving averages and were detrended by subtracting 24-hr centered moving averages. The highest and lowest levels of bioluminescence in each cycle were defined as the peak and trough, respectively. The circadian periods were calculated by using three peaks and three troughs in each dish. For calculating the amplitudes of the rhythms, the areas under the 4 curves (from the 1st trough to the 2nd peak) of the detrended data were divided by the length of time in this range.

NIH3T3 cells plated on 35-mm dishes were treated with 0.1 µM (final) dexamethasone (Sigma), and after 2 hrs, the media were replaced with a fresh one. The cell lysates were collected every 4 hr from 14 hr after the medium change and lysed with TRIzol reagent (Invitrogen). Total RNAs were prepared with the RNeasy Mini Kit (Qiagen) according to the manufacturer’s protocol.

Mice (C57BL/6 background) were handled in accordance with the Guidelines for the Care and Use of Laboratory Animals at The University of Tokyo. Mice were housed in cages with free access to commercial chow (CLEA Japan) and tap water. The animals were maintained in a light-tight chamber at a constant temperature (23 ± 1 °C) and humidity (55 ± 10%). Unless otherwise mentioned, the ΔRRE mutant mice #1 were used for the experiments in this study.

Cas9 protein (Guide-it recombinant Cas9 protein, TaKaRa) and sgRNA (Fasmac) were delivered by electroporation to C57BL/6 J embryos at the pronuclear stage as described previously⁴⁹ with slight modifications. Briefly, the embryos were washed three times with Opti-MEM I (Thermo Fisher Scientific) supplemented with 0.1% polyvinylalcohol (PVA) and once with 0.1% PVA-Opti MEM I containing Cas9 protein (100 ng/µL) and sgRNAs (50 ng/µL each). Then the embryos were placed in a line in the gap of an electrode (LF501PT1-10, BEX, Tokyo, Japan) filled with 0.1% PVA-Opti-MEM I (total 5 µL in volume) containing Cas9 protein and sgRNAs, and electroporation was performed using a CUY21EDIT II electroporator (BEX). After electroporation, survived embryos were cultured in modified Whitten’s medium for a night and were transferred into oviducts of 0.5-day-post-coitum recipients.

Male mice (5–22 wk old) were individually housed in cages each equipped with a running wheel and were entrained to the 12-hr light/12-hr dark (LD) cycles. Wheel revolutions were recorded under the LD cycles for two weeks and subsequently DD condition for 24 days or longer. Their spontaneous locomotor activities were recorded as the number of the wheel revolutions in 5-min bins and were analyzed with ClockLab software (Actimetrics). The circadian period of the activity rhythm under the DD condition was determined via a chi-square periodogram of the wheel running activities from animals that showed rhythmicity with P < 0.001, based on the activity in days 11–24 after the start of DD condition.

Tissues from 8- to 24-wk-old male mice were collected every 4 hr from 26 hrs after the beginning of the DD condition (expressed as projected CT). The livers, kidneys, hearts, and eWATs with 3–4 biological replicates were placed in Lysing Matrix D tubes (MP Biomedicals) with TRIzol reagent (Invitrogen), and homogenized by using Fast-Prep 24 instrument (MP Biomedicals). Total RNAs were prepared with the RNeasy Mini Kit (Qiagen) according to the manufacturer’s protocol.

For quantification of gene expression, RNAs prepared from NIH3T3 cells and the mouse tissues were reverse transcribed by Go Script Reverse Transcriptase (Promega) with both anchored (deoxythymine) 15 primer and random oligo primers. For quantification of mRNA levels, the cDNAs were subjected to StepOnePlus Real-Time PCR Systems (Applied Biosystems) by using GoTaq Master Mix (Promega) with specific primers of target genes (Table S2). The mRNA levels of clock genes were normalized by that of Rps29. The maximum value of WT was normalized to 1.0.

Bioluminescence signals derived from luciferase expression in the SCN, lung and eWAT cultures were recorded as described previously⁵⁰ with some modifications. The Bmal1-ΔRRE mice were bred with PER2::LUC knockin mice⁵¹. The coronal SCN slices (200 µm thick) at the middle of the rostro-caudal axis were prepared from the PER2::LUC-WT and the PER2::LUC-ΔRRE mice by using a vibrating blade tissue slicer (Neo LinearSlicer MT, Dosaka EM). The lung and eWAT samples were prepared by using a razor blade. The tissues were cultured on membranes (Millicell-CM, Millipore) in 35-mm dishes with the recording media. The bioluminescence signals were continuously recorded for 6 to 8 days at 37 °C in air with LumiCycle (Actimetrics).

Nuclear proteins were isolated as previously described²³. Briefly, the mouse liver (1 g wet weight) was washed with ice-cold phosphate-buffered saline and homogenized on ice with 9 mL of ice-cold buffer A (10 mM Hepes-NaOH [pH 7.8], 10 mM KCl, 0.1 mM EDTA, 1 mM dithiothreitol [DTT], 1 mM phenylmethylsulfonyl fluoride [PMSF], 4 µg/mL aprotinin, 4 µg/mL leupeptin, 50 mM NaF, and 1 mM Na₃VO₄). The homogenate was centrifuged (700 ×g for 5 min), and the resultant precipitate was rinsed twice. The precipitate was resuspended in 2 mL of ice-cold buffer C (20 mM Hepes-NaOH [pH 7.8], 400 mM NaCl, 1 mM EDTA, 5 mM MgCl₂, 2% [vol/vol] glycerol, 1 mM DTT, 1 mM PMSF, 4 µg/mL aprotinin, 4 µg/mL leupeptin, 50 mM NaF, and 1 mM Na₃VO₄). After being gently mixed at 4 °C for 1 hr, the suspension was centrifuged (18,000 × g for 30 min), and the supernatant was used as the nuclear fraction.

Immunoprecipitation was performed as previously described²³ with some modifications. The nuclear fraction of the mouse liver was diluted with two volumes of buffer D1 (20 mM Hepes-NaOH [pH 7.8], 5.5 mM NaCl, 1 mM EDTA, 6.5% [vol/vol] glycerol, 1.5% [vol/vol] Triton X-100, 1 mM DTT, 1 mM PMSF, 4 µg/mL aprotinin, 4 µg/mL leupeptin, 50 mM NaF, and 1 mM Na₃VO₄). Both the lysates of NIH3T3 cells and diluted nuclear fractions were incubated at 4 °C for 2 hrs with anti-CLOCK antibody²³ (CLNT1; D334-3; Medical & Biochemical Laboratories). Protein G-Sepharose 4 Fast Flow (Amersham Biosciences) was then added to this mixture, and mixed gently at 4 °C for 1 hr. The beads were collected by centrifugation for 2 min at 360 × g as the immunoprecipitants.

Proteins separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS/PAGE) were transferred to polyvinylidene difluoride membrane (Millipore). The blots were incubated in a blocking solution: 1% [wt/vol] skim milk in TBS (50 mM Tris-HCl, 140 mM NaCl, 1 mM MgCl₂ [pH 7.4]) for 1 hr at 37 °C and then incubated overnight at 4 °C with a primary antibody diluted 1:1,000 in the blocking solution. The signals were visualized by an enhanced chemiluminescence detection system (PerkinElmer Life Sciences). The blot membrane was subjected to densitometric scanning, and the band intensities were quantified by using ImageQuant TL (GE Healthcare). To evaluate the phosphorylation rhythm of BMAL1 protein, we calculated the intensity of the up-shifted band (i.e., phospho-BMAL1) relative to the total band intensity at each time point. Uncropped and unprocessed scans of all the blots were provided in Supplementary Figs. 17, 18. The primary antibodies used in this study were as follows: anti-CLOCK²³ (CLSP3; D333-3; Medical & Biological Laboratories), anti-ARNTL²³ (B1BH2; D335-3; Medical & Biological Laboratories), anti-PER1 (PM091; Medical & Biological Laboratories), anti-PER2 (PM083; Medical & Biological Laboratories), anti-CRY1 (PM081; Medical & Biological Laboratories), anti-CRY2 (PM082; Medical & Biological Laboratories), anti-DBP (PM079; Medical & Biological Laboratories), anti-NR1D1 (PM092; Medical & Biological Laboratories), and anti-ATF2 (C-19; sc-187; Santa Cruz Biotechnology).

Enzymatic digestion of mouse liver proteins was performed as previously described⁵² with some modifications. The livers from 9- to 24-wk-old male mice were collected every 4 hr from 26 hrs after the beginning of the DD condition (expressed as projected CT) and the liver blocks frozen in liquid nitrogen were lysed with ice-cold PTS buffer (12 mM sodium deoxycholate [DOC], 12 mM sodium N-lauroylsarcosinate [SLS], 100 mM NH₄HCO₃) containing 1% [vol/vol] Phosphatase Inhibitor Cocktail 2 and 3 (Sigma) followed by extensive sonication in an ultrasonic disruptor (UR-21P; TOMY; output power 5 for 2 min). After centrifugation at 18,000 ×g, protein concentrations in the supernatants were quantified by Pierce 660 nm Protein Assay Reagent, and the supernatants were diluted into 2.5 mg/mL proteins with the PTS buffer. Proteins (2.5 mg/mL, 100 µL) in the solution were reduced with final 10 mM DTT at 60 °C for 30 min, and then alkylated by incubation with final 22 mM iodoacetamide [IAA] at 37 °C for 30 min in the dark. The resultant protein sample was diluted with 100 mM NH₄HCO₃ solution up to 1 mL and digested with trypsin (Sigma) at 1:100 [wt/wt] by incubation at 37 °C for 18 hrs in the dark. After the digestion, an equal volume of ethyl acetate was added to the sample, and the mixture was acidified with final 0.5% [vol/vol] TFA and then well mixed in order to transfer the detergents into the organic phase. After the sample was centrifuged at 18,000 × g for 1 min at room temperature, the aqueous phase containing peptides was collected. The sample was concentrated by a centrifugal evaporator (EYELA) and desalted using a MonoSpin C18 column (GL Sciences). After 10% of the eluate was dried by the evaporator, it was analyzed by LC-MS/MS as a sample labeled “total peptides.” For the phospho-proteomics, the remaining eluate was dried and applied to High-select Fe-NTA phosphopeptide enrichment kit (Thermo Fisher Scientific). The enriched sample, labeled “phospho-peptides,” was dried and subjected for LC-MS/MS analysis.

LC-MS/MS-based proteome analysis was performed as previously described⁵² with some modifications. The dried and desalted peptides were dissolved in distilled water containing 2% acetonitrile and 0.1% TFA. The LC-MS/MS analyses were performed using a mass spectrometer (Q Exactive Plus, Thermo Fisher Scientific) equipped with a nano ultra-HPLC system (Dionex Ultimate 3000; Thermo Fisher Scientific). The peptides were loaded to the LC–MS/MS system with a trap column (0.3 × 5 mm L-column ODS; Chemicals Evaluation and Research Institute) and a capillary column (0.1 × 150 mm L-column ODS; Chemicals Evaluation and Research Institute) at a flow late of 20 µL/min. The loaded peptides were separated by a gradient using mobile phases A (1% formic acid in distilled water) and B (1% formic acid in acetonitrile) at a flow late of 300 nL/min (0% B for 5 min, 0–30% B for 150 min, 30–50% B for 10 min, 50–95% B for 0.1 min, 95% B for 9.8 min, 95–0% B for 0.1 min, and 0% B for 5 min). The eluted peptides were electrosprayed (2.0 kV) and introduced into the MS equipment (positive ion mode, data-dependent MS/MS). Each of the most intense precursor ions (up to the top 10) was isolated and fragmented by higher collision energy dissociation (HCD) with the normalized collision energy (27%). For full MS scans, the scan range was set to 350–1500 m/z at a resolution of 70,000, and the automatic gain control (AGC) target was set to 3e6 with a maximum injection time of 60 ms. For MS/MS scans, the precursor isolation window was set to 1.6 m/z at a resolution of 17,500, and the AGC target was set to 5e5 with a maximum injection time of 100 ms. The Orbitrap mass analyzer was operated with the “lock mass” option to perform shotgun detection with high accuracy. The raw spectra were extracted using Proteome Discoverer 2.2 (Thermo Fisher Scientific) and searched against the mouse SwissProt database (TaxID 10,090 and subtaxonomies, v2017-10-25) with following settings. The parameter of the cleavage was set to trypsin, and the missed cleavage was allowed up to 2. The mass tolerances were set to 10 ppm for the precursor ion and 0.02 Da for the fragment ion. As for protein modifications, we set carbamidomethylation (+57.021 Da) at Cys as static (fixed) modifications for peptide, oxidation (+15.995 Da) at Met and phosphorylation (+79.966 Da) at Ser and Thr as dynamic (non-fixed) modifications for peptide, and acetylation (+42.011 Da) at amino-terminus as a dynamic modification for protein terminus. The amount of each peptide was semi-quantified using the peak area with Precursor Ions Quantifier in Proteome Discoverer 2.2. According to the manufacturer’s procedure, relative levels of proteins were quantified based on unique peptides without shared peptides. Comparison of the phosphorylated amino acids between the present study and a previous study⁵³ revealed similar proportions of phosphorylated residues, i.e., 85.2% phospho-Serine (pS) and 14.8% phospho-Threonine (pT).

GO analysis was performed with open source software ShinyGO v0.61⁵⁴ with the following settings: search species, mouse; P value cutoff (FDR), 0.05; gene sets, GO Biological Process. A hierarchical clustering tree was obtained from the top 20 significant terms.

KSEA analysis was performed as previously described⁴⁴ with some modifications. The resources for kinase-substrate relationship dataset used in this analysis were PhosphoSitePlus (ver. 6.5.9.3) (https://www.phosphosite.org/homeAction↗), RegPhos2.0 (http://140.138.144.141/~RegPhos/index.php↗) and the Universal Protein Resource (Uniprot) (https://www.uniprot.org↗). For quantitative analysis, we used all the phosphorylation sites that were identified in our phospho-proteome analysis on the mouse livers. The amount of each substrate peptide was normalized by an average of the amounts of the substrate peptide in WT samples, and the normalized values were represented as log2 fold changes (positive and negative values indicate an increase and decrease in abundance of the phosphorylated peptide, respectively). These values were subjected to a Python package, Kinact (ver. 0.3)⁵⁵, in order to estimate protein kinase activities by using the kinase-substrate relationship dataset. A kinase activity score was defined as the average of the log2 fold changes of substrate peptides.

To analyze the circadian rhythmicity, datasets were run through BIO_CYCLE (ver0.9.3)⁵⁶ obtained from CircadiOmics web portal. The range of the circadian period was set to 20-28 hr according to the default settings. Unless otherwise mentioned, rhythmic genes were clarified with three criteria: (1) BIO_CYCLE P-value is less than 0.3, (2) When the expression levels are averaged at each time, the highest expression level (Peak) is 1.5 times higher than the lowest expression level (Trough), and (3) The largest SD value among each time point (SD_MAX) is within 1.5 times the deference between the Peak and Trough (Amplitude).

RNA-Seq analysis was performed as previously described²⁶ with minor modifications. The livers from the mutant and WT male mice harvested at six time points throughout the day (CT2, 6, 10, 14, 18, and 22; n = 2) were placed in Lysing Matrix D tubes (MP Biomedicals) with TRIzol reagent (Invitrogen), and homogenized by using Fast-Prep 24 instrument (MP Biomedicals). Total RNAs were prepared with the RNeasy Mini Kit (Qiagen), followed by ethanol precipitation. Poly(A)-tailed RNA was isolated from the total RNA by using NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, #E7490) and NEBNext Ultra II Directional RNA Library Prep Kit (NEB, #E7760). The library was sequenced on a NovaSeq 6000 (150 bp, pair end).

We used a mathematical model to investigate how the circadian rhythms of clock molecules can be maintained even in the absence of the RREs of Bmal1 genes. Their rhythms were simulated by using the Kim–Forger model¹⁴. The Kim-Forger model is a detailed mathematical model of the intracellular mammalian circadian clock, which describes the reactions among clock molecules (e.g., phosphorylation and binding) in the SCN by using ordinary differential equations based on mass action kinetics (181 variables and 75 parameters)¹⁴. The model accurately captured various phenotypes of circadian mutations (e.g., Rev-erbα^−/−). To simulate the rhythms in the ΔRRE mutant cells, the RRE-mediated transcriptional regulation of Bma1 is eliminated in the model. That is, we adjusted the model to make Bmal1 mRNA be expressed constitutively whose level is the same to the averaged Bmal1 mRNA level in the WT model (Fig. 3a) based on the experimental data (Fig. 1c). Similarly, to simulate the rhythms in the Bmal1- and Clock-ΔRRE mutant cells in Supplementary Fig. 8a-c, the RRE-mediated transcriptional regulation of both Bmal1 and Clock is eliminated. That is, we adjusted the model to make Bmal1 mRNA and Clock mRNA be expressed constitutively whose levels are the same with the averaged levels of Bmal1 mRNA and Clock mRNA (Supplementary Fig. 8a).

To simulate the inhibition of FBXL3-dependent degradation of CRYs in Fig. 4a and Supplementary Fig. 9, we decreased the degradation rate constant for CRY1 and CRY2, uro and urt, in the nucleus, respectively, as done in the previous work³⁴. To simulate the increase in transcriptional activity of REV-ERBs in Fig. 4d, we decreased the unbinding rate constant for REV-ERBs to Bmal1 RREs and Cry1/Clock RREs, unbinrevb and unbinrev, respectively, and thus decreased the dissociation constant of REV-ERBs to Bmal1 RREs and Cry1/Clock RREs.

In Fig. 5a, we randomly perturbed all the values of the model parameters, which represent the transcription rates, translation rates, degradation rates, unbinding rates, phosphorylation rates, and nuclear localization rates, by 50%. Then, we checked whether the WT model and the ΔRRE mutant model simulate rhythmic expression of clock molecules. Specifically, to check whether their level oscillates, we measured the relative amplitude (i.e., maximum abundance minus minimum abundance divided by the maximum abundance) of a core clock gene, Per2. The rhythmicity criteria were set from weak rhythmicity (the relative amplitude >0.1) to strong rhythmicity (the relative amplitude >0.4). Note that as the expressions of clock molecules are tightly interlocked in the model, the result of the analysis changes little even if other clock gene expression profiles are used to calculate the relative amplitude.

Further information on research design is available in thelinked to this article. Nature Research Reporting Summary

Supplementary Information Reporting Summary Description of Additional Supplementary Files Dataset 1 Dataset 2 Dataset 3 Dataset 4 Dataset 5