Technical Comments

Comment on “Circadian rhythms in the absence of the clock gene Bmal1

See allHide authors and affiliations

Science  16 Apr 2021:
Vol. 372, Issue 6539, eabe9230
DOI: 10.1126/science.abe9230


Ray et al. (Reports, 14 February 2020, p. 800) report apparent transcriptional circadian rhythms in mouse tissues lacking the core clock component BMAL1. To better understand these surprising results, we reanalyzed the associated data. We were unable to reproduce the original findings, nor could we identify reliably cycling genes. We conclude that there is insufficient evidence to support circadian transcriptional rhythms in the absence of Bmal1.

Recently, Ray et al. (1) reported transcriptional rhythmicity in mouse tissues lacking BMAL1. BMAL1 is a core component of the circadian molecular oscillator (2) whose deletion is associated with loss of physiological and molecular rhythms (3). Transcriptional circadian rhythms are thus unexpected in Bmal1–/– knockouts.

Several experiments were described in (1) using mouse skin fibroblasts from Bmal1–/– knockout and wild-type mice after synchronization with dexamethasone (DEX), including:

(i) Identification of cycling genes, with samples taken every 3 hours for 72 hours at 37°C, starting 48 hours after DEX;

(ii) Investigation of temperature compensation, with samples taken every 2 hours for 48 hours at 27°, 32°, or 37°C, starting 48 hours after DEX;

(iii) AM/PM phase shift experiment, with samples taken every 3 hours for 48 hours at 37°C, starting 48 hours (“AM”) and 60 hours (“PM”) after DEX.

RNA sequencing data were made publicly available on GEO (GSE111696, GSE134333) as both raw data and processed FPKM (fragments per kilobase per million mapped reads) values.

A common condition among these three experiments (37°C, 48 hours after DEX) provides an opportunity to test whether cycling transcripts could be reproducibly detected. An additional check of consistency comes from experiment iii, in which we expected that both conditions would reveal a common set of cycling genes, phase-shifted by 12 hours. Genes under true circadian control should be reliably detected as cycling across these experiments.

We examined the reproducibility of rhythmic genes in these experiments. Because lists of statistically significant cycling transcripts were not published, we reanalyzed Ray et al.’s data using the authors’ stated protocols: “Data were filtered by removing any transcript that had a zero FPKM value in any of the time points, and the FPKM values were subsequently log2 transformed . . . RAIN (default mode, independent method) [(4)] was used to detect rhythmic transcripts . . . using the following parameters: period = 24-hour, period-delta = 0 or 6, FDR < 0.1 or p-value < 0.05” (1, 5). The default ABH method was used for RAIN’s within-gene adjustment (4), followed by Benjamini-Hochberg FDR (false discovery rate) adjustment across the multiplicity of genes (6).

We analyzed the results in three ways. We compared the number of cycling transcripts to those reported in (1). We compared the cycling statistics across the three experiments to investigate the reproducibility of the cyclers. Finally, we tested whether cycling genes in experiment iii exhibited the expected phase shift.

Our analysis failed to reproduce the findings in (1). We identified far fewer cycling transcripts in experiment i than were reported in (1) [Fig. 1A; compare to figure 1F of (1)]. For example, we found 185 genes with FDR < 0.1 in Bmal1–/– versus 745 reported in (1). Although there was some ambiguity regarding their methods, no reasonable choice of parameters reproduced their counts. Because our analysis used the published FPKM values, discrepancies between our analysis and (1) cannot be attributed to differences in alignments.

Fig. 1 Reanalysis of transcriptional cycling in Bmal1+/+ and Bmal1–/– mouse skin fibroblast cells across three experiments performed by Ray et al.

(A) Genes detected as cycling in mouse skin fibroblasts at various significance levels for experiment i. P values were adjusted for multiple waveform tests using the ABH method of RAIN, followed by FDR correction across the multiplicity of genes. Relative to the values reported in figure 1E of (1), we identify far fewer cycling transcripts in both the wild type and the Bmal1–/– knockout. (B) Upset plot showing overlap in genes detected as cycling at FDR < 0.1 in all arms of experiments i to iii; blue and green bars indicate wild type and Bmal1–/–, respectively. Darker bars indicate the common 37°C, 48 hours post-DEX control condition; lighter bars indicate the experimental conditions (27° and 32°C in experiment ii; PM in experiment iii). The minimum set overlap is truncated at 25 for readability. Note that most genes detected as cycling are unique to a particular dataset and not reproduced by others. Inset: Untruncated upset plot for the common 37°C, 48 hours post-DEX control condition. (Overlaps between control arms may include genes that also overlap in non–control arms, shown in the main plot.) Note again that few genes are detected reproducibly among all three datasets despite the common conditions. Upper right: Scatterplots (upper triangle), histograms (diagonal), and rank correlation coefficients of FDR-adjusted P values for the common genes in the 37°C condition of the three experiments. Red lines in the upper triangle indicate FDR = 0.1, with points displayed on a –log10 scale; genes above the lines on this plot have FDR < 0.1. In the lower triangle, conditional probabilities of having FDR < 0.1 along the y axis given FDR < 0.1 along the x axis (and vice versa) are also shown. (C) Relationship among AM FDR, PM FDR, and AM-PM phase difference in experiment iii. Significant genes at FDR < 0.1 lie above 1 on the –log10 scale; the most significant cyclers lie close to the upper right corner. Highly significant genes tend to have large phase differences in Bmal1+/+; this is not observed in Bmal1–/–.

The three experiments have a common control condition (constant darkness, 37°C, sampled beginning 48 hours after DEX treatment), for which only the length and frequency of sampling differ. Although sampling choices can have an impact on cycling detection, the sampling schemes in the three experiments have comparable performance (7). Hence, genes under circadian control should be reproducibly detected in all three experiments. However, few genes were consistently detected as cycling in all three experiments: Only 26 transcripts in the wild type and four transcripts in Bmal1–/– passed FDR < 0.1 (Fig. 1B), of which only one is an annotated gene (Clk4). This suggests a lack of reproducibility for the vast majority of “rhythmic” transcripts. The overlap of results can be further quantified by the conditional probability that a probe is rhythmic (at FDR < 0.1) in one dataset if it is rhythmic in another (8). Perfect reproducibility would be 1, whereas reproducibility no better than chance would be 0.1. The observed conditional probabilities are low, ranging from 0.01 to 0.3 with a median of 0.12 (approximately chance). We also find low rank-correlation among the FDR-adjusted P values, implying that the highly ranked cyclers in one experiment are generally not highly ranked in the others. This nonreproducibility is observed in both the Bmal1–/– knockout and the wild type; notably, the wild-type and Bmal1–/– P values are no more concordant within genotypes than across the genotypes (Fig. 1B, inset).

In experiment iii, we expected that circadian genes would be detected in both the AM and PM conditions, although phase-shifted (reflecting the earlier DEX pulse). Although Ray et al. show “representative” rhythmic genes displaying a phase shift [figure 2, F and G, of (1)], no systematic analysis was reported. We calculated the number of genes passing FDR < 0.1 and FDR < 0.05 in both AM and PM, as well as their phase differences (Table 1 and Fig. 1C). There was relatively little overlap in rhythmicity between AM and PM for either genotype, despite expectations that these should be concordant. Of genes detected as cycling at FDR < 0.1 in both AM and PM, the majority exhibited much smaller phase differences than expected in Bmal1–/– (Table 1). This effect was even more pronounced at FDR < 0.05, in which strong cyclers tended to have the expected large phase differences in Bmal1+/+ but not in Bmal1–/– (Table 1 and Fig. 1C). This suggests that the oscillations observed in Bmal1–/– are not synchronized by the DEX pulse.

Table 1 Distribution of phase differences between AM and PM conditions for genes that are significantly cycling at FDR < 0.1 (no brackets) and FDR < 0.05 (brackets) in both the AM and PM conditions.

The conditional probabilities that a gene is cycling with FDR < 0.1 in AM given cycling with FDR < 0.1 in PM (and the reverse, PM given AM) were less than 0.2 in both genotypes, suggesting relatively little overlap in cycling genes generally between AM and PM. The majority of the strong cyclers exhibit the expected large phase differences in the wild type but small phase differences in Bmal1–/– (median phase difference among FDR < 0.05 genes is 9 hours in the wild type versus 3 hours in Bmal1–/–). The distributions differ significantly between Bmal1+/+ and Bmal1–/–, with phases in Bmal1+/+ higher than in Bmal1–/– (P < 0.0001, Fisher’s exact test).

View this table:

Note that the original version of (1) indicated that the RAIN “longitudinal” method was used and did not mention filtration of genes; this was later updated (5). Results using the original settings may be found in our code repository; although the output differs, the key findings (lack of reproducibility and internal consistency) are unchanged.

Our analysis suggests that the various genes detected as “cycling” in any given study are likely false positives (7, 8) or experimental artifacts. The fact that the majority of rhythmic genes are not reproducibly cycling is consistent with previous observations that RAIN has a high false-positive rate (7, 8); however, we observed a similar lack of reproducibility and internal consistency using JTK_CYCLE (9) and ARSER (10). Spurious cycling may also arise from interactions with the environment or between cells (11). Although we cannot exclude the possibility that there may be circadian dynamics in Bmal1–/–, we do not find any evidence of consistently observed circadian rhythms in these data.


Acknowledgments: We thank M. Rosbash and J. Hogenesch for sharing their independent analyses of the Ray datasets, and the reviewers for helpful comments. Funding: Supported by NSF grant 1764421 and Simons Foundation grant 597491-RWC as part of the NSF-Simons Center for Quantitative Biology at Northwestern University. Author contributions: E.N.C., R.A., and R.B. conceived the study; E.N.C. and R.B. carried out the analysis. All authors contributed to the manuscript. Competing interests: None. Data and materials availability: Data were obtained from GEO (GSE111696, GSE134333). R code to reproduce our results (along with additional analyses) is available from

Stay Connected to Science

Navigate This Article