To investigate the extent to which ∼12-h rhythms of gene expression are present and functionally preserved in the absence of the circadian clock in vivo, we performed a post hoc analysis of a recently published temporal RNA-seq data collected from the liver of male whole body BMAL1 knockout mice (3 months of age) (Bunger et al., 2000) that were kept under constant darkness after being entrained under 12-h/12-h light/dark cycle for 14 days (Aviram et al., 2022). This specific dataset was selected due to the following reasons. First, Bmal1 is the only gene that when knocked out in mice leads to a complete abolishment of circadian locomotor and feeding cycles, and as a result BMAL1 knockout mice is a widely-accepted model for circadian clock ablation (Bunger et al., 2000). Secondly, in this study, BMAL1 knockout mice were kept in constant darkness, thus ruling out the possibility that any ultradian rhythms hereby uncovered are driven by light cues. Thirdly, of all the available temporal transcriptome dataset from BMAL1 knockout mice, it has the best temporal resolution (at 4-h interval), duration (for a total of 48 h) and sample size (n of four), thus permitting rigorous testing for ∼12-h ultradian rhythms.

To unbiasedly reveal all oscillations from BMAL1 knockout mice liver, we initially applied the eigenvalue/pencil method (Zhu et al., 2017; Antoulas et al., 2018), which is a nonparametric spectrum analytical algorithm that deconvolutes raw temporal data into a linear combination of exponentials (sinusoidal waveforms with a decay factor) plus noise. Consistent with the role of BMAL1 as a master regulator of circadian rhythms, the prevailing population of circadian oscillations observed in wild-type mice are lost in the absence of BMAL1 (Figure 1A). By contrast, 3,589 genes cycling with periods between 10–13 h were identified (Figures 1A, B; Supplementary Table S1). Since the eigenvalue/pencil algorithm is a non-statistical signal processing method, we used a permutation-based method that randomly shuffles the time label of gene expression data to determine the false discovery rate (FDR) for the ∼12-h genes (Rey et al., 2018; Pan et al., 2020). This way, the FDR for the ∼12-h genes was estimated to be 0.31, a range in line with prior studies of circadian rhythms (Mauvoisin et al., 2014; Robles et al., 2014; Rey et al., 2018). Furthermore, the calculated FDR is likely an overestimation since amplitude and phase information for the ∼12-h rhythms is not considered in the permutation dataset.

To assess the robustness of the uncovered ∼12-h oscillations in BMAL1 knockout mice, we applied an orthogonal rhythm-identification algorithm, RAIN (Thaben and Westermark, 2014), which detects rhythms with arbitrary waveforms exhibiting a single pre-specified period and has been successfully used to uncover ∼12-h ultradian rhythms from high temporal resolution dataset in prior studies (Meng et al., 2020; Pan et al., 2020). In line with the eigenvalue/pencil method, 12-h were also the dominant oscillations identified by the RAIN method, with 390, 993 and 1,557 12-h genes identified with FDR (or q value) cut-off of 0.05, 0.1 and 0.15, respectively (Figures 1C, D; Supplementary Figures S1A–C; Supplementary Table S2). Importantly, the robustness of 12-h oscillations for genes under different FDR cut-off does not look visually different (Figure 1D; Supplementary Figures S1A, B), and ∼12-h genes identified by both methods further showed large convergence (Figure 1E).

We subsequently performed GO analysis on ∼12-h genes identified by both methods using all mouse genes or all hepatically expressed genes as background, and identified commonly enriched pathways of lipid metabolism, proteostasis (including ribosome assembly, protein transport, protein processing, protein folding/stability/degradation, autophagy, and Golgi organization), and to a lesser degree, mRNA metabolism (transcription, RNA splicing and processing) (Figure 1F; Supplementary Figure S1D). Among these pathways, protein and mRNA metabolism were previously observed to be highly enriched in hepatic ∼12-h transcriptome in wild-type mice as well (Meng et al., 2020; Pan et al., 2020). Collectively, these results indicate the presence of prevalent ∼12-h ultradian rhythms of gene expression in the absence of the circadian clock. Furthermore, since the locomotor activity and feeding behaviors of BMAL1 knockout mice are completely arrhythmic under constant darkness condition (Bunger et al., 2000), it suggests that these ∼12-h rhythms of gene expression are highly robust and resistant to feeding cues.

We next aim to determine the extent to which ∼12-h transcriptome is conserved between wild-type and BMAL1 knockout mice. Those commonly found in both genotypes are likely controlled by the 12-h oscillator, while “de novo” ∼12-h genes only observed in BMAL1 knockout mice may reflect compensatory gain of new ultradian rhythms secondary to circadian clock ablation. We selected the hepatic RNA-seq dataset published by Pan et al. (2020) for comparison due to the two following reasons. First of all, this dataset has the highest temporal resolution (at 2-h interval), duration (48 h) and sample size (duplicates at each time point) among all temporal hepatic RNA-seq dataset for wild-type mice kept under constant darkness condition (Brooks et al., 2022). Secondly, the strain (C57BL/6), sex (male) and age (8–12 weeks) of these mice are comparable to those of BMAL1 knockout mice.

Comparing ∼12-h genes uncovered by the eigenvalue/pencil method in both datasets revealed a statistically significant (p-value of 1.54e-6) overlap of 1,162 commonly shared genes, indicating preservation of ∼12-h gene programs between wild-type and BMAL1 knockout mice (Figures 2A, B; Supplementary Figure S2A). The relative amplitude (absolute amplitude normalized to average gene expression) of shared ∼12-h genes is comparable between the two genotypes, albeit with a slight (but statistically significant) increase observed in BMAL1 knockout mice (Figure 2C). GO terms (regardless of the background gene sets used) associated with those 1,162 commonly shared ∼12-h genes are highly similar to those of ∼12-h genes in wild-type mice, namely, proteostasis (protein transport, Golgi organization, translation regulation, UPR, tRNA processing, and ubiquitin-mediated protein degradation) and mRNA metabolism (transcription regulation, mRNA processing and splicing and mRNA transport). The convergence of ∼12-h genes and their functionality between wild-type and BMAL1 knockout mice are further validated with the RAIN method (FDR cut-off of 0.15) (Figures 2D, E; Supplementary Figures S2A, B).

In addition to commonly shared ∼12-h genes, many genes with ∼12-h rhythms exclusively identified in either wild-type or BMAL1 knockout mice were present (Figure 2B). A close examination of their period change in the opposite genotype reveals that 2,718 ∼12-h genes only identified in wild-type mice often exhibit periods slightly lower than 10-h in BMAL1 knockout mice, thus being labelled as not having ∼12-h oscillations (Supplementary Figure S3A). These genes are also enriched in protein and mRNA metabolic pathway as expected regardless of the background gene sets used (Supplementary Figures S3B, C). On the contrary, many of the 2,218 ∼12-h genes exclusively observed in BMAL1 knockout mice exhibit periods in the circadian range in wild-type mice and these genes are strongly enriched in lipid metabolism, mitosis, complement activation and Notch signaling, pathways that are not strongly associated with ∼12-h transcriptome in wild-type mice (Supplementary Figures S3A–C). The latter group can be categorized as de novo ∼12-h genes that are gained in BMAL1 knockout mice secondary to the ablation of circadian clock. Taken together, these findings provide additional evidence to support the existence of a 12-h oscillator that controls ∼12-h rhythms of gene expression of protein and mRNA metabolism in mice liver.

∼12-h genes commonly identified in both wild-type and BMAL1 knockout mice are likely under 12-h oscillator control. In agreement with this hypothesis, their ∼12-h rhythms of gene expression are impaired in XBP1 (LKO) mice (Pan et al., 2020) (Figures 2A, D). Conversely, as previously reported, core circadian clock and the majority of circadian output genes are maintained in XBP1 LKO mice (Meng et al., 2020; Pan et al., 2020), and abolished in BMAL1 knockout mice, as expected (Figure 3A). While the repertoire and functionality of circadian and ∼12-h genes in wild-type mice are largely separate (Figures 3B, C; Supplementary Figure S4A), we did observe hundreds of genes having both circadian and ∼12-h rhythms superimposed (Figure 3C). The phenomenon of superimposition of circadian and ultradian rhythms in the same genes are widespread and have been previously reported in many different organisms by other groups (Ananthasubramaniam et al., 2018; Yang et al., 2022; Psomas et al., 2023). Using a very stringent criterion where the superimposed circadian and ∼12-h rhythms need to be abolished in BMAL1 KO and XBP1 LKO mice, respectively, but maintained in XBP1 LKO and BMAL1 KO mice, respectively, we identified 141 genes under robust dual control of BMAL1-dependent circadian clock and XBP1s-dependent 12-h oscillator (Figure 3C; Supplementary Table S3). Examples of such genes include chaperone Cct3, mitochondria import receptor Tomm40l (Tom40 b), master transcriptional regulator of fatty acid oxidation Ppara, and RNA helicase Ddx1 (Figure 3D; Supplementary Figure S4B). In all cases, eigenvalue/pencil method revealed coexistence/superimposition of both circadian and ∼12-h rhythms (sometimes with additional ∼4.8-h or ∼8-h rhythms as well) in wild-type mice, and more importantly, the linear addition of these rhythms closely matches the raw temporal gene expression profiles, indicating the robustness of deconvolution (Figure 3D; Supplementary Figure S4B). Interestingly, while in BMAL1 KO and XBP1 LKO mice, circadian and ∼12-h rhythms are no longer identified by the eigenvalue/pencil method, respectively, ultradian rhythms with even smaller periods (4.8∼8-h) are often observed (Figure 3D; Supplementary Figure S4B), indicating the likely existence of additional oscillators for high frequency rhythms.

To infer the gene regulatory network governing ∼12-h oscillator, we performed motif analysis on the commonly shared ∼12-h genes between wild-type and BMAL1 knockout mice identified by either the eigenvalue/pencil or RAIN methods. Regardless of the method used, we identified top motifs associated with KLF, ETS, NFY and basic leucine zipper (including ATF6, XBP1, CREB3, CREB3L2, and CREBRF) families of transcription factors (Figure 4A), consistent with previous results (Pan et al., 2020; Zhu, 2020). XBP1 itself exhibited a ∼13-h oscillation in the liver of BMAL1 knockout mice (Figure 4B). We further cross-referenced enriched motifs with transcription factors exhibiting ∼12-h rhythms in both wild-type and BMAL1 knockout mice, but not in XBP1 LKO mice, and identified KLF7, ELF1 and ATF6B as strong putative transcriptional regulators of 12-h oscillator. Interestingly, KLF7 was previously shown to also exhibit a ∼12-h rhythm of expression in human dorsolateral prefrontal cortex (Scott et al., 2023). In addition, we further identified NYFB and CREBRF as two transcriptional factors exhibiting ∼12-h rhythms in all three genotypes of mice (Figure 4C). Of these transcription factors, hepatic XBP1s, NFYA, NFYC and GABPB1 expressions further exhibit ∼12-h rhythms at the protein level in wild-type mice, as previously measured (Wang et al., 2018) and quantified (Dion et al., 2022), while the expression levels of GABPA/B2, ATF6/6B, NFYB, KLF7 and ELF1 were not reported from the same study.

To ensure that our analysis is robust to datasets used, we further compared the ∼12-h hepatic transcriptome of BMAL1 knockout mice with an additional wild-type dataset published by the same group (Aviram et al., 2022). The relatively low resolution (4-h interval for 48-h under constant darkness) renders this dataset not ideal for ultradian rhythms detection by statistical methods like RAIN due to the superimposition of strong circadian rhythms for many genes. Nonetheless, using eigenvalue/pencil method, we uncovered hundreds of overlapping ∼12-h transcriptome between wild-type and BMAL1 knockout mice that has 1) similar amplitude but distinct phases, 2) strong enrichment in proteostasis pathways, and 3) over-representation of KLF, bZIP, ETS or NYF TFs motifs in their promoters (Supplementary Figures S5A–H), all consistent with the results described above. Collectively, our analysis of ∼12-h rhythms in wild-type, BMAL1 knockout and XBP1 LKO mice provided novel mechanistic insights on the regulation of mammalian 12-h oscillator.

Having established the preservation of ∼12-h rhythms of gene expression in the absence of circadian clock in mice, we wanted to further corroborate this finding in additional animals that also lack a canonical circadian clock. We focused upon the fruit fly Drosophila melanogaster, where ∼12-h rhythms, to the best of our knowledge, have never been reported before. In Drosophila melanogaster, the model of the circadian clock regulation involves transcription factors CYCLE (CYC) and CLOCK (CLK), the homologs of BMAL1 and CLOCK in mammals, which in turn control the transcription of several repressive clock genes including period (per) and timeless (tim) (Panda et al., 2002). Besides a subset of cells located in the principal pacemaker driving daily activity rhythms and physiology, most other cells including Drosophila Schneider 2 (S2) cells [originally derived from a primary culture of late-stage embryos (Schneider, 1972)] do not express these core circadian clock genes and therefore do not have the necessary apparatus to form a TTFL to drive circadian oscillations (Saez and Young, 1996; Darlington et al., 1998; Rey et al., 2018). We thus deduce that if ∼12-h rhythms of gene expression with functionality like those identified in mice can be observed in S2 cells, it will provide additional evidence for the existence of a cell-autonomous 12-h oscillator independent from the canonical circadian clock.

We analyzed a recently published temporal RNA-seq data (Rey et al., 2018) collected from temperature entrained S2 cells (2 days of temperature cycles with 12-h at 28°C, 12-h at 23°C before free-running at 28°C) at 3-h interval for a total duration of 60 h. Using both eigenvalue/pencil and RAIN methods, we identified tens of hundreds of transcripts cycling with an ultradian period between 11 and 15 h (Figures 5A–D; Supplementary Tables S4, S5). This period range is longer than what was observed in mice, likely due to the likewise longer circadian period (22–28 h) identified in S2 cells (Figures 5A–C). This way, both ∼12 and ∼24-h rhythms are still cycling at harmonics frequencies with each other in S2 cells. GO analysis of these ∼12-h ultradian genes revealed top enriched pathways of protein transport (such as Rab1) and mRNA splicing (such as Pnn), distinct from those associated with circadian rhythms but similar to those associated with ∼12-h genes in mice (Figures 5E, F; Supplementary Figure S6). Among the putative/confirmed transcription regulators governing ∼12-h rhythms in mice, we found the fly orthologs of Gabpa, Son, Atf6b and Elf1 exhibiting a period between 11 and 14 h in S2 cells (Figure 6). The fly orthologs of Nfya and Xbp1 exhibit a longer period of 15–16 h in S2 cells (Figure 6). In sum, we herein provided additional evidence supporting the existence of ∼12-h rhythms of gene expression of mRNA and protein metabolism in clock-less fly cells. More importantly, the fact that ∼12-h rhythms of gene expression implicated in protein and mRNA metabolism have been found in divergent species of sea anemone (Sorek et al., 2018), C. elegans (van der Linden et al., 2010), fly, mice and human (Scott et al., 2023) argue for the strong evolutionary conservation of the ∼12-h oscillator.

For BMAL1 knockout or WT mice RNA-seq data (GSE171975), raw counts were first normalized to total sequencing depth to get counts per millions (CPM) values, which were used for cycling transcripts identification by either the eigenvalue/pencil or RAIN methods. For the eigenvalue/pencil method (Zhu et al., 2017; Antoulas et al., 2018), the mean expression at each time point were calculated and used as input. A maximum of two oscillations were identified for each gene. Criterion for ∼12-h genes: period between 10 h and 13-h, decay rate between 0.8 and 1.2. The relative amplitude was calculated by dividing the amplitude by the mean expression value for each gene. To determine FDR, we used a permutation-based method that randomly shuffles the time label of gene expression data and subjected each permutation dataset to the eigenvalue/pencil method applied with the same criterion (Rey et al., 2018). These permutation tests were run 5,000 times, and FDR was estimated by taking the ratio between the mean number of rhythmic profiles identified in the permutated samples (false positives) and the number of rhythmic profiles identified in the original data. Analyses were performed in MatlabR 2017A. RAIN analysis was performed as previously described in Bioconductor (3.4) (http://www.bioconductor.org/packages/release/bioc/html/rain.html) (Thaben and Westermark, 2014). FDR (also known as q value) was calculated using the Benjamini-Hochberg procedure.

For RNA-seq in S2 cells (GSE102495), FPKM values reported in the original study were used for cycling transcripts identification. Eigenvalue/pencil analysis was performed as previously described except for that genes with periods between 11 and 14-h are defined as ∼12-h genes. For RAIN analysis, the FPKM data was first subject to polynomial detrend (n = 2). For both datasets, heat maps were generated with Gene Cluster 3.0 and TreeView 3.0 alpha 3.0 using either Z score or log2 (fold change over mean) values.

DAVID (Version 2021) (Huang da et al., 2009) (https://david.ncifcrf.gov) was used to perform Gene Ontology analyses. Briefly, gene names were first converted to DAVID-recognizable IDs using Gene Accession Conversion Tool. The updated gene list was then subject to GO analysis using either all Mus musculus (for mice) or Drosophila melanogaster (for S2 cells) genes or hepatically expressed genes (for mice) or all genes expressed in S2 cells as background and with Functional Annotation Chart function. GO_BP_DIRECT was used as GO categories. Only GO terms with a p-value less than 0.05 were included for further analysis.

Motif analysis was performed with the SeqPos motif tool (version 0.590) embedded in Galaxy Cistrome using all motifs within the mouse reference genome mm9 as background. Genome region spanning 0.5 kb upstream and 05 kb downstream of transcription start were used as input.

Data were analyzed and presented with GraphPad Prism software. Plots show individual data points and bars at the mean and ± the standard error of the mean (SEM). One-tailed t-tests were used to compare means between groups, with significance set at p < 0.05. For representative temporal transcritptome data, the curves connecting individual data points were fitted with Cubic spline function in Graphpad. For all period-identification analysis, only the raw data points (but not the fitted Cubic curves) were used.