Genome-wide correlation analysis to identify amplitude regulators of circadian transcriptome output

Cell-autonomous circadian system, consisting of core clock genes, generates near 24-h rhythms and regulates the downstream rhythmic gene expression. While it has become clear that the percentage of rhythmic genes varies among mouse tissues, it remains unclear how this variation can be generated, particularly when the clock machinery is nearly identical in all tissues. In this study, we sought to characterize circadian transcriptome datasets that are publicly available and identify the critical component(s) involved in creating this variation. We found that the relative amplitude of 13 genes and the average level of 197 genes correlated with the percentage of cycling genes. Of those, the correlation of Rorc in both relative amplitude and the average level was one of the strongest. In addition, the level of Per2AS, a novel non-coding transcript that is expressed at the Period 2 locus, was also linearly correlated, although with a much lesser degree compared to Rorc. Overall, our study provides insight into how the variation in the percentage of clock-controlled genes can be generated in mouse tissues and suggests that Rorc and potentially Per2AS are involved in regulating the amplitude of circadian transcriptome output.

genes. Of our particular interest was Rorc, whose correlation in both the relative amplitude and the average level was one of the strongest. We also found that the level of Per2AS, a novel non-coding transcript that is expressed at the Period 2 locus, also showed a correlation. Based on these data, we propose that Rorc is involved in regulating the amplitude of circadian transcriptome output.

Methods
Microarray data processing. Microarray data were downloaded through NCBI GEO from series GSE54650 10 . Data was originated from 12 different tissues, with 24 time points from each tissue in 2-h intervals over the course of 48 h. Extracted data was normalized by Robust Multichip Average (RMA) normalization 13 and annotated by the Affymetrix Transcriptome Analysis Software package (http://www.affym etrix .com/suppo rt/techn ical/bypro duct.affx?produ ct=tac). Unannotated probesets, as well as those that had values lower than the average of all negative probesets across all timepoints in the respective tissue, were eliminated from the downstream analysis. For multiple probesets annotated to the same gene, the probeset with the highest average value was selected.
RNA-seq data processing. Mouse RNA-seq data were downloaded as fastq files through the NCBI database from SRA ID SRP036186 10 . Data contained information from 12 different tissues, with 8 time points from each tissue in 6-h intervals over the course of 48 h. Reads were mapped to the Ensembl mouse genome release 95 using STAR 2.7.0a 14 with outFilterScoreMinOverLRead = 0.3 and outFilterMatchNMinOverLRead = 0.3 options. We added the option 'Condensegenes' to select the most abundant isoform as the representative of a gene, as well as the option 'count exons' to measure only mRNA. Baboon circadian transcriptomic datasets were downloaded from series GSE54650 12 . Reads were mapped to the Ensembl baboon genome (Papio Anubis 2.0) release 90 using STAR 2.7.2b 14 with outFilterScoreMinOverLRead = 0.3 and outFilterMatchNMinOverLRead = 0.3 options. The quantification of expression level was performed by HOMER 15 using the transcripts per million (TPM) option. Any transcript with an average TPM < 0.5 across all timepoints were eliminated from the downstream rhythmicity analysis. We also used TPM to normalize the expression levels of each transcript. We eliminated white adipose data from the downstream analysis because no transcripts were rhythmic with our statistical threshold (BH. Q-value < 0.05), even though more than 13,500 transcripts were detected after applying the filter of TPM > 0.5. The expression of Per2AS was measured with the "strand-" option in HOMER. We did not apply the filter (TPM > 0.5 to call 'expressed') in quantifying the level of Per2AS, because non-coding transcripts generally have low expression levels 16-18 . Rhythmicity analysis. We used MetaCycle 19 to determine the rhythmicity of each gene. MetaCycle integrates three different algorithms ARSER, JTK_CYCLE, and Lomb-Scargle and calculates the p-value, Benjamini-Hochberg q-value (BH.Q value), period, phase, baseline value, amplitude (AMP), and relative amplitude (rAMP), which is the ratio between amplitude and baseline expression level. We defined the expression rhythmic when meta2d BH.Q < 0.05. Correlation analysis. Pearson and Spearman correlation tests were performed in R to determine the linear and non-linear correlation between the percentage of cycling transcripts in each tissue and the rhythmicity (using BH.Q value), phase, and relative amplitude of the 15 clock genes, as well as Per2AS, calculated by the MetaCycle package in R. A significant correlation was defined as a p-value < 0.05. For the transcriptome-wide correlation analysis of mouse and baboon data, we used the rcorr function from the Hmisc and tidyverse packages in R to perform Pearson or Spearman correlation tests and used the average gene expression of transcripts expressed in all 12 (microarray) or 11 tissues, excluding white adipose tissue (RNA-seq) 20,21 . Rat data expression values were Log2 normalized, and we then used JMP to perform Pearson or Spearman correlation tests. Fisher Z-scores were calculated from the Rho or R 2 with Fisher transformation. GO enrichment analysis of the significantly correlated genes was performed using the Gene Ontology Resource 22 used in this study are as follows: 36B4_Fwd 5′-CAC TGG TCT AGG ACC CGA GAAG-3′, 36B4_Rev 5′-GGT GCC TCT GAA GAT TTT CG-3′,  Rora_Fwd 5′-ACC GTG TCC ATG GCA GAA C-3′, Rora_Rev 5′-TTT CCA GGT GGG ATT TGG AT-3′, Rorc_Fwd  5′-TCT ACA CGG CCC TGG TTC T-3′, Rorc_Rev 5′-ATG TTC CAC TCT CCT CTT CTC TTG -3′. 1 μM dexamethasone was added to media for 2 h in samples used for bioluminescence recordings. After 2-h dexamethasone treatment, media was changed to phenol red-free DMEM (Cellgro 90-013-PB) supplemented with 100 μM luciferrin, 10 mM HEPES pH 7.2, 1 mM sodium pyruvate, 0.035% sodium bicarbonate, 2% FBS, 1 × Penicillin/Streptomycin, and 2 mM l-glutamine. Samples for real-time bioluminescence recordings were treated with various concentrations of nobiletin (1.5, 5, 15, or 45 μM) or SR1001 (1, 5, or 10 μM), and measurements were performed for 7 days using a LumiCycle (Actimetrics, Inc. Wilmette, IL). The first 24 h were removed from measurements to quantify amplitude using JMP software 26 .

Characteristics of circadian transcriptomic output in various mouse tissues. To gain insight
into what determines the number of clock-controlled genes in each tissue, we first retrieved existing circadian transcriptome datasets from various mouse tissues 10 . We found this particular dataset best-suited to our study, because it covered the highest number of tissues (12 total) and provided the highest time resolution (2 h intervals), compared to other studies 5,7-9 . Our in-house analysis was able to replicate the previous findings, in which the percentage of cycling genes was highest in liver, followed by kidney, lung, brown adipose, and heart, and lowest in brainstem (Fig. 1A, Supplemental Data Files 1a and 1b). The ranks are slightly different from the original study 10 , which is most likely due to the differences in the analytical methods and statistical criteria used in our study (see "Methods"). Distribution of Benjamini-Hochberg q-values from the rhythmicity analysis was also widest in liver, followed by kidney, lung, brown adipose, and heart, and was particularly narrow in white adipose and brainstem (Fig. 1B). These data indicate that the expression levels of transcripts in liver are most variable across different times of the day compared to other tissues such as white adipose or brainstem, regardless of their rhythmicity.
The difference in the percentage of rhythmic transcripts in each tissue was not due to the difference in the number of genes expressed, because the total number of transcripts detected in each tissue was comparable, irrespective of the percentage of rhythmic transcripts (Fig. 1C). It was not due to the median microarray signals per transcript in each tissue either, as it did not correlate with the percentage of rhythmic transcripts (Pearson r 2 = 0.115, p = 0.368; Spearman rho = − 0.042, p = 0.904) (Fig. 1D). When we focused only on the transcripts that were rhythmically expressed (Benjamini-Hochberg q < 0.05), the median of the relative amplitude (i.e., the ratio between amplitude and baseline expression level) ( We also performed the same set of analyses using the RNA-seq data 10 , which surveyed the same set of tissues but with a lower time resolution (Microarray: 2 h, RNA-seq: 6 h) 10 . Even though the order of the tissues was slightly different than the microarray datasets ( Fig. S1A), which was likely due to the differences in time resolution and the threshold to eliminate low-expressed transcripts (see "Methods"), the results were essentially the same: distribution of Benjamini-Hochberg q-values from the rhythmicity analysis was wider in tissues with a high percentage of rhythmic transcripts (Fig. S1B), the number of genes expressed was comparable among tissues (Fig. S1C, D), and the median of the relative amplitude or the amplitude of cycling gene expression was comparable in each tissue (Fig. S1E,F, Supplemental Data File 2).
Characterization of cycling gene expression in various mouse tissues. Because we did not observe any characteristics that had a correlation with the percentage of cycling genes at a genome-wide scale, we shifted our focus on single gene level analyses. We first analyzed a total of 15 core clock genes (Arntl, Clock, Npas2, Per1-3, Cry1-2 Rora-c, Nr1d1-2, Dbp, Nfil3), and found that most of these genes were expressed ubiquitously across all tissues, except for Rorb, whose expression was restricted to brain and brown adipose tissue (Supplemental Data File 3). A majority of these genes were also rhythmically expressed, except for the Rors: Rora was rhythmic in four tissues (liver, lung, heart, and muscle) but not in the other eight tissues (kidney, BAT, adrenal, aorta, cerebellum, hypothalamus, WAT, and brainstem), Rorb was arrhythmic in all four tissues that it is expressed in (BAT, hypothalamus, cerebellum, and brainstem), and Rorc was rhythmic in most of the peripheral tissues but not in the hypothalamus, muscle, or brainstem (Supplemental Data File 3), which was consistent with previous reports 7,9,27,28 . Clock and Cry1 were also arrhythmic in hypothalamus, making the hypothalamus the tissue with the lowest number of rhythmic core clock genes, even though hypothalamus ranked 9th out of 15 tissues in the percentage of cycling transcripts (Fig. 1A). Similar results were obtained from the RNA-seq data, although the number of rhythmic core clock genes were lower, most likely due to the lower time resolution of the RNA-seq data compared to the microarray data (Supplemental Data File 4).
As was previously reported, the phases of core clock gene expression were confined to a relatively narrow window 12,28-31 , except for a few genes such as Cry2, Rorc, and Nfil3 ( Fig. 2A). On the other hand, the relative Scientific Reports | (2020) 10:21839 | https://doi.org/10.1038/s41598-020-78851-9 www.nature.com/scientificreports/ amplitude of core clock gene expression was more variable between tissues (Fig. 2B), and nine genes, Dbp, Npas2, Nr1d1, Arntl, Per3, Per2, Rorc, Cry1, and Cry2 had their relative amplitude positively correlated with the www.nature.com/scientificreports/ percentages of cycling transcripts in either Pearson and Spearman correlation analyses (Table 1). Additional ten clock-controlled genes were expressed and rhythmic in all tissues, for which we calculated the correlation between their relative amplitude and the percentage of cycling transcripts. Among those, the relative amplitude of three genes (P4ha1, Tsc22d3, and Lonrf3), or another set of three genes (Tsc22d3, Lonrf3, and Usp2) correlated significantly with the percentage of cycling transcripts in Spearman or Pearson analyses, respectively (Table 1). Notably, the strongest correlation was observed for Rorc (Pearson r 2 = 0.852, p = 0.004; Spearman rho = 0.917, p = 0.001). In addition, the number of rhythmic core clock genes also correlates with the percentages of cycling transcripts in Spearman analysis (rho = 0.728, p = 0.007), but not in Pearson analysis (r 2 = 0.510, p = 0.0906). The sum of relative amplitudes from rhythmic core clock genes also correlates with the percentages of cycling transcripts (Pearson r 2 = 0.653, p = 0.0212; Spearman rho = 0.734, p = 0.0091). It is unclear, however, whether the higher amplitude of core clock gene expression leads to a higher percentage of rhythmic transcripts in each tissue or vice versa. We also analyzed the data from RNA-seq and performed the same analyses. However, the lower www.nature.com/scientificreports/ number of rhythmic core clock genes detected in the RNA-seq dataset significantly compromised our ability to calculate the correlations between the percentage of rhythmic transcripts in each tissue and the phase and amplitude of core clock gene expression. We also investigated the correlation between the mean levels of core clock gene expression across all time points and the percentage of cycling genes in each tissue, as non-cycling genes could contribute to the differences in the percentage of cycling genes. Interestingly, we again found that there was a positive correlation between the percentage of cycling transcripts and the mean level of Rorc, but not with any other core clock genes (Fig. 2C, Table 2). This correlation was also observed in the RNA-seq dataset (Fig. S2, Table 2). We also tested Per2AS, a newly identified non-coding RNA [32][33][34] , because the expression of Per2AS is rhythmic and antiphasic to Per2 in liver, adrenal gland, lung, and kidney 10 , and it was suspected that Per2AS is involved in regulating the circadian system 35 . Indeed, we found a linear correlation between the mean levels of Per2AS and the percentages of rhythmic transcripts in both microarray and RNA-seq datasets (Fig. 2D). Single value decomposition (SVD) analysis did not detect any clear clustering patterns for tissues (Fig. 2E). In contrast, all the clock genes were clustered together except for Rorc (Fig. 2E, Fig. S2B), indicating that the correlation between Rorc and the percentage of rhythmic transcripts is not due to innate differences between the tissues. Rather, it is most likely due to the difference in expression patterns of Rorc between tissues that are distinct from other clock genes.
To test how robust the correlation of Rorc is, we extended the analysis to the genome-wide scale. We found that among 12,024 genes expressed in all 12 tissues from microarray datasets, the mean level of 1131 and 400 genes was correlated significantly with the percentage of rhythmic genes in each tissue from the Pearson (linear) and the Spearman (non-linear) correlation test, respectively (Supplemental Data File 5). Similarly, among 8269 genes expressed in all 11 tissues from the RNA-seq datasets, we found the mean level of 925 and 1664 genes correlated significantly with the percentage of rhythmic genes from the Pearson and Spearman correlation tests, respectively (Supplemental Data File 5). Of those, 135 (Pearson) or 77 (Spearman) genes were found correlated in both microarray and RNA-seq datasets (Supplemental Data File 5), and we therefore considered those as more robustly correlated. Gene ontology (GO) analyses were then performed to assess whether a specific process contributes to the high percentage of rhythmic transcripts. No pathways were detected as statistically significant (FDR < 0.05) among those that correlated robustly in the Spearman analysis. Whereas numerous metabolic processes were enriched among those that correlated robustly in the Pearson analysis (Supplemental Data File 6). We also calculated Fisher Z-scores from each test to evaluate the relative strength of Rorc correlation, compared to other genes. Rorc was ranked 5th (Pearson) or 14th (Spearman) when we used average Z-scores from both microarray and RNA-seq datasets. These data suggest that the correlation between the level of Rorc and the amplitude of the mouse circadian transcriptome is one of the strongest.
To further gain more mechanistic insights into how Rorc contributes to the increase in the number of cycling mRNAs without driving mRNA expression, we next tested the correlation between the mean levels of Rorc and other core clock genes. Not surprisingly, we found a linear correlation between the mean levels of Rorc and Per2AS in both the microarray and RNA-seq datasets ( Table 3). The mean level of Rorc also linearly correlated with Nfil3 (Microarray) or Cry2 (RNA-seq) ( Table 3); however, the biological significance of these correlations is Table 1. Correlations between the percentage of rhythmic transcripts and the relative amplitude of core clock genes in each tissue. *Asterisks denote p < 0.05. a Rorb was excluded from our correlation analyses due to its low expression in all tissues except brain. www.nature.com/scientificreports/ unclear, as the correlations were not consistent between microarray and RNA-seq datasets. We did not detect any correlation between Rorc and Rorc-target genes such as Arntl, Cry1, Nfil3, and Nr1d1 (Table 3), whose promoter regions contain RORE motifs and the amplitude of their rhythmic mRNA expression was dampened in most of the Rorc −/− tissues [36][37][38][39] . We did not find any correlations between the level of Per2AS and other core clock genes either, except for Per2 and Npas2 (Table 4). The significance of these correlations is also unclear, because they were found only in microarray but not in RNA-seq datasets.

The effect of RORC as a transcriptional activator in regulating the circadian transcriptome.
Because Rorc directly activates the transcription of RORE-containing genes 27,29,40 , we hypothesized that, if Rorc was directly driving rhythmic gene expression leading to a high number of cycling transcripts, then the number of rhythmic genes with RORE motifs in their promoter would be higher in tissues with a higher number of rhythmic genes. To test this, we first retrieved the promoter sequence of rhythmic genes in each tissue Table 2. Correlations between the percentage of rhythmic genes and the mean expression level of core clock genes in microarray and RNA-seq datasets. *Asterisks denote p < 0.05. a Rorb was excluded from our correlation analyses due to its low expression in all tissues except brain.  Table 3. Correlations between the mean level of Rorc and the mean level of other clock genes in each tissue. *Asterisks denote p < 0.05. a Rorb was excluded from our correlation analyses due to its low expression in all tissues except brain. www.nature.com/scientificreports/ using the RNA-seq dataset (− 1000 to + 100 bp with respect to TSS), and then determined the number of DNA motifs that can be recognized by RORC in silico. We surveyed the recognition sequences of not only RORC, but also RORA, RORB, NR1D1, NR1D2, ARNTL, CLOCK, NPAS2, NFIL3, and DBP, because these proteins are all considered important to drive rhythmic gene expression 39,41 . We found that the RORC binding motif was found in approximately 8% of the rhythmically expressed genes, and this was consistent in all 11 tissues we examined (Fig. 3). Of the genes with RORC motifs in their promoter region, we found 960 genes that were expressed in all 11 tissues and rhythmic in at least one tissue. Of these 960 genes, the mean level of 291 and 252 genes correlated with the mean level of Rorc by Spearman or Pearson analysis, respectively (Supplemental Data File 7). It is plausible that RORC regulates the expression and/or rhythmicity of these downstream RORC target genes. The binding sites for NPAS2, ARNTL, NR1D1, and NR1D2 were the most highly represented (~ 10-14%), and NFIL3 and DBP were the least represented (~ 3-6%) (Fig. 3). Table 4. Correlations between the mean Per2AS TPM and the mean level of other clock genes in each tissue. *Asterisks denote p < 0.05. a Rorb was excluded from our correlation analyses due to its low expression in all tissues except brain.  www.nature.com/scientificreports/ To experimentally test the biological significance of Rorc, we utilized three independent luciferase reporter cell lines; Per2::LucSV cells (mouse embryonic fibroblasts derived from Per2::LucSV knockin mice) 42 as well as Bmal1-luc and Dbp-luc cells (derived from NIH3T3 cells that stably express luciferase genes under the control of either a Dbp or Bmal1 promoter) 43 . We treated these cells with the Rora/c agonist nobiletin or inverse agonist SR1001 to determine whether it is the mRNA level or the transcriptional activity of RORC that is important for the amplitude of core clock genes expression patterns (Fig. 4A). We first analyzed the mRNA level of Rora/c in Per2::LucSV cells, as the level of Rorc cells was under the detection limit in Dbp-luc or Bmal1-luc cells, both of which derived from NIH3T3 cells 44 . We found that the mRNA levels of Rorc were unchanged in the presence of either nobiletin or SR1001, whereas the mRNA levels of Rora decreased when treated with nobiletin (Fig. 4B). We also found that nobiletin increased the amplitude of PER2::LUCSV and Dbp-luc reporter output compared to control cells (Fig. 4C). In contrast, SR1001 decreased the amplitude of Dbp-luc reporter output. Interestingly, neither nobiletin nor SR1001 altered the amplitude of Bmal1-luc despite Bmal1 being under direct control of Rora/c (Fig. 4C) 27 . Because both nobiletin and SR1001 have a higher inhibition constant for RORC compared to RORA 45,46 , the changes in the amplitude of reporter output are likely due to their effect on RORC.

Discussion
Among the three key parameters for cycling systems (period, phase, and amplitude), the regulatory mechanisms of period and phase have been relatively well-characterized, whereas that of amplitude have remained much more enigmatic. Forward genetics or screening approaches using pharmacological or genetic perturbation have not been successful, as the variance of amplitude is much higher than that of period, compromising the statistical www.nature.com/scientificreports/ ability to distinguish true positives from false positives [47][48][49][50][51] . Furthermore, it is currently unclear whether it is a single gene, a combination of genes, one of the feedback loops, and/or topology of the network that is important for amplitude. To make things even more complicated, amplitude can be measured by various outputs, such as gene expression, firing patterns in neurons, body temperature, and locomotor activity, each of which can be under the regulation of both cell-autonomous (intracellular or local) and systemic (or extracellular) rhythms. It also remains unclear whether all the amplitude of various rhythms is regulated by the same mechanism.
In this study, we focused on the percentage of cycling genes in various mouse tissues and explored the possible mechanisms of amplitude regulation of circadian transcriptomic output. The circadian transcriptome can be influenced both by cell-autonomous and systemic cues in each tissue. However, the circadian gene expression is a common feature of the circadian clock system in various tissues and allows us to directly compare the difference in amplitude between tissues without relying on their respective physiology.
We found that 18 genes (eight core clock genes and ten clock-controlled genes) are rhythmically expressed in all tissues that we examined (Table 1, Supplemental Data File 3). This is consistent with previous observations that the rhythmicity of each gene is often tissue-specific, and while only a handful of genes are cycling in all or most tissues, others are rhythmic only in certain tissues [7][8][9][10][11] . The mechanism that drives tissue-specificity of rhythmic gene expression still remains largely unknown. Nonetheless, recent studies demonstrated that the BMAL1 DNA binding is largely tissue-specific, and the tissue-specific rhythmic gene expression can be driven by rhythmic BMAL1 binding in coordination with the activity of enhancers nearby that form chromatin looping 52,53 . It is possible that RORC drives tissue-specific rhythmic target gene expression using a similar mechanism. We also found that the relative amplitude of 13 genes (nine core clock genes and four clock-controlled genes) were correlated with the percentage of cycling genes, while the mean level of 197 genes correlated with the percentage of cycling genes. These genes are not necessarily expressed rhythmically, albeit about a half (100/197) are, and vast majority of these genes were involved in the metabolic processes (Supplemental Data Files 5,6). Given that the energy cost for cycling genes are higher than non-cycling genes 54 , it is reasonable that metabolism related genes are highly expressed in tissues that have higher number of cycling transcripts. Interestingly, metabolically active tissues, such as liver, brown fat, and skeletal muscle have rhythmic RNAs with higher amplitude (Fig. 1E), which is consistent with the previous report 53 , suggesting that the metabolic activity in each tissue affects the amplitude of rhythmic RNA expression. It is also plausible that systemic cues, including metabolites, contributed to the differences in the percentage of rhythmic gene expression in each tissue, and these genes act as mediators to connect systemic cues and rhythmic gene expression both independent and dependent on the core clock machinery. However, it is more tempting to postulate that the role of Rorc and/or Per2AS in the core-clock circuit gives it a more promising function in potentially regulating the amplitude of the circadian transcriptome, at least in the tissues where Rorc is expressed.
We also attempted to verify our findings using independent datasets. To achieve this, datasets must include gene expression data taken at multiple circadian time points (to detect rhythmicity) in at least three different tissues (to perform correlation analyses). To date, no other datasets are available in mouse. One study using rats analyzed gene expression patterns in four tissues (liver, lung, muscle, and adipose) at multiple circadian time points 11 (Supplemental Data File 8). Unfortunately, however, some of the microarray platforms used in this study did not include probes for Rorc and we were unable to calculate the correlation between the percentage of rhythmic genes and the level of Rorc. A recent study in Baboon also surveyed diurnal gene expression patterns in 64 tissues 12 . We analyzed 14 tissues from this dataset (Aorta, Adrenal Cortex, Adrenal Medulla, Bone Marrow, Heart, Hippocampus, Kidney Cortex, Kidney Medulla, Liver, Lung, Pancreas, Prostate, Smooth Muscle, and SCN) that closely mirror 12 mouse tissues (or 11 for RNA-seq) that we analyzed in our study.
When we used the MetaCycle p-value < 0.05 as a rhythmicity threshold, we did observe a positive correlation between the level of RORC with the percentage of rhythmic transcripts (Pearson: r 2 = 0.603, p = 0.029, Fig. S3A, Supplemental Data Files 9a, 9b). However, this correlation was not statistically significant when we used the JTK_CYCLE algorithm with p-value < 0.05 (Hughes et al., 2010) (Pearson: r 2 = 0.167, p = 0.587, Fig. S3B) or when we used q-value < 0.05 as a rhythmicity threshold. In contrast, the correlation between the Rorc level and the percentage of rhythmic transcripts in mouse was statistically significant with both MetaCycle (Table 2, Fig. 2C) and JTK_CYCLE (Pearson: r 2 = 0.738, p = 0.006, Spearman: rho = 0.678, p = 0.019) (Fig. S3C) with B.H. Q value < 0.05. Of note, we were also unable to detect PER2AS in any of these tissues we examined (Fig. S3D), indicating that either PER2AS does not exist in baboon, or the baboon reference genome we used (Papio Anubis 2.0) does not have high enough resolution to annotate and detect PER2AS. Overall, we concluded that it is unclear whether the RORC level and the percentage of rhythmic transcriptome correlate in baboon.
We have also found that the mRNA levels of 584 and 473 genes also correlated with the percentage of the rhythmic mRNAs in each tissue in rat and baboon, respectively (Supplemental Data File 7) 11,12 . Among these, 40 and 4 genes also showed a correlation in mouse dataset 10 whereas none was commonly detected as correlated in all three datasets (Supplemental Data File 7). The significance of these genes in regulating the amplitude of circadian transcriptome output is unclear, however. Experimental parameters (sampling resolution, number of tissues examined, transcriptomic platform) are significantly different between studies, and these have a significant impact on detecting rhythmicity of each gene as well as calculating the percentage of rhythmic transcripts (Supplemental Data Files 1a, 1b, 2, 8, 9a, 9b). It is also possible that the regulatory mechanism of amplitude is species-specific and not conserved.
The positive loop (Clock-Arntl-Rev-Ror) was originally considered to confer additional robustness to the system and, therefore, stabilizes the system. However, it is not required for circadian rhythm generation 37,55 . Recent studies have also highlighted the role of the positive loop as the central axis of amplitude regulation 45,56 . Our study is consistent with these findings and suggest that the positive loop, particularly the level of Rorc, is important in setting the amplitude of the circadian transcriptome. In addition, our study also suggested that Per2AS is involved in the positive loop, because the level of Per2AS positively correlated with the level of Rorc as www.nature.com/scientificreports/ well as the percentage of cycling genes in each tissue, even though it was originally assumed to only interact with Per2 35 . Interestingly, the mathematical model predicted the functional interaction between Ror and Per2AS, as Per2AS would restore stable circadian rhythms when they are disrupted by the overproduction of Ror or Rev-erb mRNAs 35 . It is possible that Rorc and Per2AS function synergistically in the circadian clock system. Because long non-coding RNAs, such as Per2AS, can function in trans and interacts with DNAs, other RNAs, or proteins to regulate target gene expression 57,58 , it is possible that Per2AS RNA interacts with RORC proteins to modify its transcriptional activity. Per2AS could also alter the level of Rorc by interacting with the promoter or enhancer sequences of Rorc (i.e., DNA), transcriptional regulators and/or epigenetic modifiers of Rorc or vice versa. Alternatively, the correlation between Rorc and Per2AS may simply indicate that their expression is regulated by the same or similar mechanism. It is still unclear from our study whether the relationship between Rorc and the percentage of rhythmic transcripts is simply a correlation or causation. Because the molecular clock system is quite complex, and the expression of each clock gene is dependent on the expression/activity of other clock genes directly or indirectly, we think it is critically important to keep the network intact to fully understand the function of each component. In contrast to experimental approaches in which perturbation of each genetic component (i.e., gene knock-out, knock-down, or overexpression) would often disturb the network itself, a computational approach excels in this area to shed light into the function of a component within a network. To circumvent these issues, we used a pharmacological approach and altered the transcriptional activity of RORC without changing its levels. We found that the transcriptional activity of RORC alters the amplitude of reporter bioluminescence output (Fig. 4), suggesting that the transcriptional activity of RORC is important for regulating the amplitude of the circadian clock machinery and potentially the circadian transcriptome output. Regardless, we think the level of Rorc is still biologically relevant, because the changes in the Rorc mRNA level can lead to the changes in the level of RORC and/or its activity.
It still remains unclear why Rorc, but not Rora and Rorb, correlates with the amplitude of the circadian transcriptome, as all the ROR proteins share significant sequence similarities 40,59 . Unfortunately, the physiological relevance of each ROR paralogue has never been clarified in the mammalian circadian system. One notable difference among Ror paralogues, however, is their unique expression patterns (Supplementary Data Files 1a, 1b, 2). It is possible that the systemic cues, which are, in theory, the same to all the tissues have a tissue-specific impact in regulating the level of Rorc. Understanding the difference in the regulatory mechanisms of Ror gene expression would provide insight into how their tissue-specific expression pattern is achieved and how Rorc specifically impacts the amplitude of the circadian transcriptomic output.
Overall, our study highlighted the potential role of Rorc in regulating the amplitude of the circadian transcriptome, although it is unclear whether the correlation between the Rorc and the percentage of rhythmic transcriptome is specific to mouse. Follow-up experimental studies would further complement our observations from the rich transcriptomic datasets that are publicly available and delineate the mechanisms of circadian amplitude regulation.

Data availability
The datasets generated during and/or analyzed during the current study are available in the NCBI GEO repository, from series GSE54650 (mouse), GSE98965 (baboon), or GSE8988, GSE8989, GSE20635, and GSE25612 (rat).