Analysis of Gene Regulatory Networks in the Mammalian Circadian Rhythm

Identification of a Common Set of Circadian Genes in Mouse

We searched for circadian oscillating genes in 21 circadian time series microarray data covering 14 tissues in mouse (Table S1) by fitting them to cosine functions with different phases, and extracted circadian phase information for circadian oscillating genes. We identified 9,995 known genes showing circadian oscillations in at least one tissue (Table S2). The number of genes showing circadian oscillation in multiple tissues decreases rapidly as the number of tissues increases, whereas the consistency of their circadian phases across tissues as measured in p-values of circular range tests improves rapidly (Figure 1). We identified 41 common circadian genes, defined as the genes showing circadian oscillation in at least 8 out of 14 tissues in mouse (Table 1). 13 out of 19 previously known key circadian genes were among the common circadian genes that we identified in this study. Other known key circadian genes: Rorb, Cry2, Rora, Npas2, and Hlf were found to be circadian oscillating in one, three, three, four, and five tissues, respectively. Bhlhb3 was not found to be circadian oscillating in any tissue. 39 of these common circadian genes showed significant consistency (p<1/3 in circular range test) of their circadian phases across all tissues.

Comparison between Tissues

We surveyed tissue-specific gene expression profiles in a mouse tissue gene expression atlas [7] for the circadian oscillating genes in different tissues. To cross-validate the circadian phase data with the tissue gene expression data, we created a binary matrix of 1 or 0 to denote the presence or absence of circadian oscillations in 14 tissues in circadian phase data and compared it to the gene expression matrix in 61 tissues from the tissue gene expression atlas. For each pair of tissues from the two matrices, we calculated a correlation coefficient. The circadian data in liver, kidney, skeletal muscle, adrenal gland, and white adipose tissue correctly correlated best with their corresponding tissues in the tissue gene expression atlas, whereas SCN correlated equally well with preoptic and hypothalamus, and brown adipose tissue correlated equally well with adipose tissue and brown fat. These results reflected the fact that sufficiently high gene expression levels are the prerequisite to be detected as circadian oscillating in our collection of microarray datasets.

To investigate if the differences in the circadian phases of circadian oscillating genes across tissues are caused by the differences in their gene expression levels, we calculated the variances of circadian phases and the variances of gene expression for circadian oscillating genes across the seven tissues common to our circadian datasets and the tissue gene expression atlas. There is no significant correlation (r = 0.01, p = 0.71) between these two variances. For example, the gene expression level of Per2 is 27 times higher in adrenal gland than in skeletal muscle, but this has no effect on the consistency of circadian phases of Per2 between the two tissues. In fact, the common circadian genes have significantly higher variances of gene expression across the 61 tissues than those from the same number of randomly selected genes. We observed that the correlation coefficients r_ij between the tissue gene expression data of the common circadian gene pairs (i,j) negatively correlated with their circadian phase differences (r = −0.22, p<10⁻⁸). The gene pairs positively correlated in their tissue gene expression patterns had a significantly lower circadian phase difference than expected by random, whereas the gene pairs negatively correlated in their tissue gene expression patterns had a significantly larger circadian phase difference than expected by random (Figure S1). Therefore, the common circadian genes with similar gene expression patterns across tissues also tend to have similar circadian phases. The circadian gene regulation may share a similar mechanism that gives rise to tissue-specific gene expression.

We clustered the 21 circadian phase datasets using hierarchical clustering. The datasets from the same tissue or biologically closely related tissues were clustered together, suggesting that the differences in circadian phases between tissues resulted from their biological differences (Figure 2). To ensure that these differences between tissues were also reproducible between experiments, we used circular ANOVA to identify the circadian oscillating genes shared between two tissues but associated with significantly different circadian phases between these tissues. There were 12 circadian oscillating genes shared between two SCN datasets and at least two liver datasets. Among them, Per1, Per2, Nr1d2, and Avpr1a showed a significant (p<0.01) advance of about 6 hours in their circadian phases in SCN datasets compared to liver datasets, whereas Dnajb1, Hmgb3, Hsp110, and Pdcd4 showed no significant differences in their circadian phases between SCN and liver (Figure 3). To test if such differences also exist between SCN and whole brain tissues, we also compared SCN with 3 whole brain datasets. There were 12 circadian oscillating genes shared between two SCN datasets and at least two whole brain datasets. Per2, Nr1d2, and Tuba8 again showed a significant advance of about 6 hours in their circadian phases in SCN datasets compared to whole brain datasets, whereas Hmgb3, Hsp110, Sgk, and Fabp7 showed no significant differences in their circadian phases between SCN and whole brain. Further examination validated that the known key circadian genes including Per1, Per2, Cry1, Arntl, Nr1d1, and Nr1d2 all showed around 6 hour advances in circadian phases between SCN and non-SCN tissues in general, whereas heat shock proteins showed consistent circadian phases across all tissues. There were 15 circadian oscillating genes shared between 3 heart datasets including whole heart, atria, and ventricle and at least 3 liver datasets. Comparing the heart datasets with the liver datasets, Bhlhb2 (p<0.001) and Tspan4 (p = 0.006) had circadian phase 5–6 hours earlier in heart than liver whereas Dscr1 (p = 0.002) had circadian phase 8 hours later in heart than liver. Other known key circadian genes such as Per1/Per2, Arntl, and Nr1d1/Nr1d2 showed consistent circadian phases between heart and liver. Comparing the whole brain datasets with the liver datasets, Tfrc, St3gal5, and Tspan4 had circadian phases more than 4 hours earlier in whole brain than liver, whereas Hist1h1c, Tsc22d1, Myo1b, Litaf, and BC004004 had circadian phases more than 4 hours later in whole brain than liver.

Comparison between Mammalian Species

Among the 1,269 rat genes identified as circadian oscillating genes in rat liver, 1,137 of them had homologues in mouse. 232 of them overlapped with 944 mouse liver circadian oscillating genes in at least 2 mouse liver datasets. We used the circular ANOVA test to identify the circadian oscillating genes shared in both mouse and rat livers but with significantly different circadian phases. 10 genes had significantly (p<0.01) different circadian phases between mouse and rat livers. The circadian phases of BC006779, Cdkn1a, Svil, Uox, Ak2, Nr1d1, Mtss1, Nudt16l1, and Gss were 4–6 hours later in rat liver than mouse liver, whereas Hsd17b2 was in anti-phase between mouse and rat livers (Figure S2).

Among 803 rat skeletal muscle (SKM) circadian oscillating genes, 703 of them had homologues in mouse and 64 of them overlapped with 440 mouse SKM circadian oscillating genes. Among the overlapping genes, 34 of them did not show circadian phase differences larger than 4 hours between mouse and rat SKM. 22 of them had circadian phases more than 4 hours later in rat SKM than mouse SKM. Cpt1a, Pdk4, and Ucp3, involved in lipid metabolism, showed a 5–8 hour delay in their circadian phases in rat SKM compared to mouse SKM. 8 genes had circadian phases more than 4 hours earlier in rat SKM than in mouse SKM. Among them, Fkbp5 and Sgk, which are controlled by the glucocorticoid receptor element (GRE), had about 6 hour advance in their circadian phases in rat SKM compared to mouse SKM. There were 11 circadian oscillating genes common to mouse liver and SKM, and rat liver and SKM. The 4–5 hour delay in circadian phases in rat compared to mouse was observed in both liver and SKM for all 11 circadian genes except Dynll1.

Among 603 rhesus macaque adrenal gland circadian oscillating genes, 560 had homologues in mouse and 170 overlapped with 4,162 mouse adrenal gland circadian oscillating genes. We found significant differences in circadian phases also between these two species. Among the overlapping genes, 47 did not show circadian phase differences larger than 4 hours between mouse and macaque, whereas 66 had circadian phases more than 4 hours later in the macaque adrenal than in the mouse adrenal. Known key circadian genes, Arntl, Dbp, Nr1d1, and Bhlhb2, showed about 8 hour delay in their circadian phases in the macaque adrenal compared to the mouse adrenal. Although Per2 did not satisfy our criteria (p<0.01) to be a circadian oscillating gene in macaque adrenal, this gene has a circadian phase at CT21 (p = 0.03), which is also about 8 hours later than that in mouse. Similarly, heat shock proteins, Hsp110, Hspa8, Dnaja1, and Dnajb6, had circadian phases around CT16 in the mouse adrenal but around CT0 in the macaque adrenal. Cold inducible protein (Cirbp) had a circadian phase around CT7 in the mouse adrenal but around CT16 in the macaque adrenal, in anti-phase with heat shock proteins in both mouse and macaque. On the other hand, there were also 57 genes showing circadian phases more than 4 hours early in the macaque adrenal than in the mouse adrenal.

In the human circadian SKM microarray study, there were only two circadian time point measurements: CT1 and CT13. Hence we can only roughly estimate the circadian phases to be either CT1 or CT13 in human SKM. Among the common circadian genes, Per1, Per2, Nr1d2, and Dbp had circadian phases around CT1, whereas Arntl and Cry1 had circadian phases around CT13 in human SKM. Our estimates of circadian phases for Per1 and Per2 in human SKM were in good agreement with the study in human peripheral blood mononuclear cells where a 2 hour sampling time was used throughout 72 hours [8]. The heat shock proteins, Dnaja1, Dnajb4, and Hspa4, had circadian phases around CT13, consistent with the peak of common body temperature at CT10 in human [8].

Next, we made a three-species comparison of circadian phases in the SKMs of mouse, rat, and human. We found 12 circadian oscillating genes common to SKM in all three species (Table 2). After we rounded the circadian phases in mouse and rat to their closest time points, CT1 or CT13, we observed that Per2, Arntl, Dbp, Ppp1r3c, and Ablim1 had conserved circadian phases between mouse and rat, but were 12 hours away from those of human. Epm2aip1, G0S2, and Maf had conserved circadian phases between mouse and human but 12 hours away from those of rat. Finally, D19Wsu162e, Myod1, Pfn2, and Ucp3 had conserved circadian phases among all three species.

Biological Functions of the Circadian Rhythm

We searched for the Gene Ontology (GO) categories significantly over-represented in circadian oscillating genes in each mouse tissue using GOminer program [9]. We further tested the associations of GO categories with any specific circadian phase intervals using Fisher's test with a rotating window method. The list of significant biological processes associated with circadian phases in different tissues is shown in Table S3. The most common of these biological processes were steroid biosynthesis, heat shock response, and protein folding. Steroid biosynthesis was associated with CT22 in liver, kidney, adrenal, brown adipose tissue (BAT), and white adipose tissue (WAT). Heat shock response or protein folding were associated with CT16 in SCN, liver, kidney, adrenal, aorta, BAT, WAT, calvarial bone, and whole brain, due to a large number of heat shock proteins consistently showing circadian phases near CT16 in most tissues. In liver, carbohydrate and amino acid metabolism were associated with CT17 and CT15 respectively, consistent with the rise of activities after light off in mouse. In BAT, WAT, and adrenal, lipid metabolism was associated with CT22. Negative regulation of protein kinase activities was associated with CT17 in prefrontal cortex and CT21 in whole brain. There were also notable differences in the circadian phases of some biological processes between tissues. For example, protein translation was associated with CT20 in SCN but CT9 in WAT. Organ development was associated with CT22 in heart and BAT but CT10 in adrenal.

Promoter Analysis

To test the association of transcription factor (TF) regulation with the circadian oscillation of gene expression, we predicted the TF binding sites on the mouse promoters of circadian oscillating genes in each tissue using positional weight matrix (PWM) based methods. We first tested whether there was a significant over-representation of TF PWM binding sites on the promoters of circadian oscillating genes using the Fisher's exact test. Among the significant TF PWMs, we again tested their associations with any specific phase intervals using the Fisher's test with a rotating window method. To remove the redundancy in TF PWMs, we grouped the TF PWMs into TF families and averaged the associated circadian phases of significant TF PWMs within the same TF families. The results are shown in Table S4. EBOX, AP-2, CRE, SP1, and EGR were the top 5 TF families associated the circadian phase in most tissues. However, unlike the consistent circadian phases of the common circadian genes across tissues, the associated circadian phases of the significant TF families varied considerably among different tissues. EBOX was associated with CT12 in the majority of tissues including SCN, liver, aorta, adrenal, WAT, brain, atria, ventricle, and prefrontal cortex, but it was associated with CT0 in skeletal muscle, BAT, and calvarial bone. CRE was consistently associated with CT11 in SCN, liver, aorta, heart, adrenal, calvarial bone, prefrontal cortex, and ventricle, but with CT20 in atria. Two other known TF families related to circadian rhythm, RRE and DBOX, were detected to be associated with circadian phase only in two tissues. RRE was associated with CT0 in liver and WAT. DBOX was associated with CT16 in aorta and adrenal.

Identification of Gene Regulatory Interactions

We obtained microarray data from TF knockout or mutants for Clock, Arntl, Npas2, Nr1d1, Rora/Rorc, Egr1/Egr3, Dbp/Hlf/Tef, and Ppara in various mouse tissues, together with Cebpa/Cebpb/Cebpd/Cebpe transfection microarray data in NIH3T3 cells. To study the systematic effects of glucocorticoids, cAMP, and temperature on the circadian rhythm, we included microarray data from Nr3c1 (glucocorticoid receptor), Pka, and Hsf1 knockouts or mutants in response to DEX (glucocorticoid agonist), cAMP, and heat stimulation, respectively, compared with wild type mouse. We also included microarray data from a light response mouse model in order to identify light sensitive genes in mouse SCN [10]. The complete list of knockout or mutant microarray experiments used in this study is shown in Table S5. We assumed that the target genes of TFs will be significantly down-regulated in the knockout or mutant compared with the wild type mouse in the case of activators, and up-regulated in the case of repressors, such as Nr1d1. To identify the direct targets of TFs in knockout or mutant experiments, we required that the significantly affected genes in the knockout or mutant must have at least one putative binding site of their corresponding TFs in the promoter regions. Under these criteria, we identified 320 EBOX, 295 RRE, 43 DBOX, 492 EGRE, 455 CRE, 326 GRE, 122 HSE, 607 CEBP, and 516 PPRE controlled genes respectively (Table S6). For these genes, we extracted their mean circadian phases if they have consistent circadian phases across multiple tissues (p<1/3, circular range test). We observed that EBOX was significantly associated with CT12 (p<10⁻⁶, Fisher's exact test), RRE with CT1 (p<10⁻⁶), DBOX with CT15 (p<10⁻⁵), HSE with CT17 (p<10⁻⁶) (Figure S3).

Circadian Gene Regulatory Network

Based on these regulatory interactions, we constructed the gene regulatory network for the circadian oscillating genes in mouse. In Figure 4, we show a network consisting of the circadian oscillating genes identified in at least 7 mouse tissues. Among the 81 circadian oscillating genes identified in at least 7 tissues, 53 of them can be included through 88 regulatory interactions with 9 cis-regulatory elements in our network. Their circadian phases were represented by different colors in the color wheel. We were able to identify almost all known transcription regulatory interactions for common circadian genes in the literature, except EBOX → Per1, EBOX → Nr1d1, EBOX → Ppara, RRE → Nr1d1, and RRE → Cry1. To further complete our network, we supplemented these missing gene regulatory interactions with known protein interaction information (Per/Cry Arntl/Clock and Fkbp:Hsp90 Nr3c1) and protein phosphorylation information (Csnk1d → Per/Cry and Gsk3b → Nr1d1) from the literature. These relationships are shown in red color in Figure 4.

Two well-known negative feedback loops can be reconstructed from this analysis: Arntl/Clock → EBOX → Per1/Per2 Arntl/Clock and Nr1d1/Nr1d2 RRE → Arntl/Clock → EBOX → Nr1d1/Nr1d2. Two feedforward loops are attached to the negative feedback loops through Arntl/Clock → EBOX → Dbp → DBOX → Per1/Per2 acting as an alternative route of Arntl/Clock → EBOX → Per1/Per2 and Nr1d1/Nr1d2 RRE → Nfil3 DBOX → Per1/Per2 Arntl/Clock acting as an alternative route of Nr1d1/Nr1d2 RRE → Arntl/Clock. Bhlhb2 inhibiting EBOX is also regulated by EBOX and Nr1d1 inhibiting RRE is also regulated by RRE, therefore forming two auto-regulatory loops.

The effects of food and light act on common circadian genes directly through GRE and CRE respectively. GRE controls Per1 and Per2, while CRE controls Per1, Rora, Nr1d2, and Nfil3. As shown in Figure 4B, the effect of temperature acts on common circadian genes rather indirectly through the route HSE → Hsp90aa1 → Fkbp/Hsp90 Nr3c1 → GRE → Per1/Per2. Nr3c1 and the Fkbp/Hsp90 complex are also components of another negative feedback loop, Nr3c1 → GRE → Fkbp5 → Fkbp/Hsp90 Nr3c1, which may play an important role in glucocorticoid stimulation. Nr3c1 is also under the control of CRE and therefore may be responsive to light stimulation. Nr3c1 and the Fkbp/Hsp90 complex feed into EBOX by regulating Per1/Per2 through GRE. In turn, EBOX controls both components of the Fkbp/Hsp90 complex, i.e., Fkbp5 directly and Hsp90aa1 indirectly through EBOX → Ppara → PPRE → Hsp90aa1. Therefore, Nr3c1 and Fkbp/Hsp90 play central role of integrating the regulatory inputs from diverse environmental signals into circadian genes in our network (Figure 4B).

Circadian Microarray Data

We collected all available circadian microarray data from different laboratories for mouse, rat, rhesus macaque (Macaca mulatta), and human. The total mouse data consisted of 21 datasets covering 14 tissues including two datasets in SCN, five datasets in liver, three datasets in whole brain, one dataset in kidney, aorta, heart, skeletal muscle (SKM), adrenal gland, brown adipose tissue (BAT), white adipose tissue (WAT), calvarial bone, prefrontal cortex, atria, and ventricle. The three datasets in whole brain were from three different mouse strains: C57BL/6J, AKR/J, and DBA/2J. The rat data consisted of one dataset in liver and one dataset in skeletal muscle. The macaque data consisted of one dataset in adrenal gland. The human data consisted of one dataset in skeletal muscle. The complete list of all circadian microarray datasets used in this study is shown in Table S1. Most circadian microarray experiments were conducted in a time series of every 4 hours. The human microarray experiment was only conducted at CT1 and CT13. For simplicity, we did not distinguish the light conditions, i.e., 12 h light:12 h dark (LD) or 12 h dark:12 h dark (DD), under which the animals were kept during the experiments. In order to have a more complete and consistent analysis of the data from different experiments, we decided to re-analyze all the datasets by our own method rather than simply taking the gene lists from the original publications. For the datasets where the CEL files were available, we normalized the data by RMA method in “affy” package. For the datasets where only normalized data were available, the normalization step was skipped.

We used the method similar to that described in [2] to analyze all microarray data. Namely, cosine functions A_ij(t) = cos(2πt/T_i−ϕ_j) where T_i = 20+i, ϕ_j = 2πj/60, 0 ≤ i ≤8, and 0 ≤ j ≤59 were used as the reference time series of circadian oscillation. The gene expression time series of each probe set on the microarray were fitted to each cosine function time series A_ij(t) and the cosine function with highest correlation coefficient was chosen. A p-value <0.01 in the regression for the best cosine function was used as the criterion for circadian oscillation, and we estimated a false positive rate of about 10% for this cutoff using a random permutation test. When the experimental replicas at each time point were available, we further carried out a one-way ANOVA test on the time series using time points as factor and p-value <0.05 as an additional criterion. For the probe sets satisfying the criteria for circadian oscillation, the gene expression time series were again fitted to the cosine functions with fixed 24 hrs period but changing phases, B_j(t) = cos(2πt/T−ϕ_j), where T = 24, ϕ_j = 2πj/144, and 0≤ j ≤ 143. The circadian phase was calculated from the best fitted B_j(t) as ϕ_{j *}24/2π. We were unable to obtain the microarray data in [2] so we only extracted circadian gene lists with their circadian phase information. In the human SKM study, vastus lateralis muscles were taken from exercised and non-exercised legs of 4 patients at CT1 (8AM) and CT13 (8PM). We used circadian time and exercise state as two factors in two-way ANOVA. A p-value <0.05 in the circadian time comparison was used as the criterion for circadian oscillation. We estimated the circadian phase to be either CT1 or CT13, depending on when the average expression value was the highest in human SKM.

The R package “Circular” was used to analyze the circadian phases obtained from circadian microarray datasets. For each circadian microarray dataset, the probe sets were annotated by R package “annaffy” and only the probe sets corresponding to known genes were used in the analysis. The probe sets that passed circadian oscillation criteria and that corresponded to the same genes were merged by the following procedure. First, a circular range test was used to assess the consistency of phases estimated from the different probe sets for the same genes, where p<1/3 was used as the criterion to take into account the 4 hour intrinsic errors in phase estimation as the animals were sampled every 4 hours in most experiments. Then, a circular mean function was used to calculate the mean circadian phases from the consistent probe sets. The same procedure was used to combine the different datasets for the same tissue, i.e., five datasets for liver, two datasets for SCN, three datasets for whole brain. In liver and whole brain, we only selected the genes identified as circadian oscillating in at least two out of five liver datasets or two out of three whole brain datasets, respectively. In SCN, we selected the genes identified as circadian oscillating in one out of two SCN datasets considering the small number of circadian oscillating genes in Ueda et al.'s SCN dataset [2]. We identified 9,955 circadian oscillating genes in at least one out of 14 tissues (Table S2). The number of circadian oscillating genes in different number of tissues was plotted using the “barplot” function in R and is shown in Figure 1A. The circular range test was also used to describe the consistency of phases of circadian oscillating genes across tissues. The distribution of p-values of circular range tests in different number of tissues was plotted using the boxplot function in R and is shown in Figure 1B. We defined the 41 circadian oscillating genes identified in at least 8 out of 14 tissues as common circadian genes, and these are shown in Table 1.

Tissue-Specific Gene Expression of Circadian Oscillating Genes

The microarray data of 61 mouse tissues after gcrma normalization were downloaded from the mouse tissue gene expression atlas website: http://symatlas.gnf.org [7]. We selected the probe set with the highest average expression value across tissues to represent the genes with multiple probe sets. To remove the non-detected probe sets, we filtered out the probe sets with gene expression values lower than 100 in all 61 tissues. We obtained the expression profiles for 19,168 genes across 61 tissues. For 9,955 circadian oscillating genes identified in at least one tissue, we created a matrix of 1 or 0 to denote the presence or absence of circadian oscillation in 14 tissues. For 8,029 genes having both circadian data and tissue expression data, we calculated the correlations between the circadian 1 or 0 matrix with the matrix of log₂(gene expression) in 61 tissues in tissue gene expression atlas. We searched for the tissues in tissue data having the highest correlation coefficient with the tissues in circadian data. Liver (r = 0.29, p<10⁻¹⁵), kidney (r = 0.23, p<10⁻¹⁵), skeletal muscle (r = 0.10, p<10⁻¹⁵), adrenal gland (r = 0.06, p = 10⁻⁷), and white adipose tissue (r = 0.18, p<10⁻¹⁵) in circadian data have the highest correlations with their corresponding tissues in the tissue data, whereas SCN in circadian data correlates equally well with preoptic and hypothalamus (r = 0.22, p<10⁻¹⁵) in tissue data and BAT correlates equally well with adipose tissue and brown fat (r = 0.19, p<10⁻¹⁵). For the seven tissues having both circadian data and tissue data: liver, heart, BAT, WAT, kidney, adrenal gland, and SKM, we calculated the variances of circadian phases in circadian data using the “circular var” function for the circadian oscillating genes identified in at least two tissues, and the variances of log₂(gene expression) in tissue data across the tissues where the circadian oscillations have been identified in circadian data. The correlation coefficient of these two variances is 0.01 (p = 0.71). For the 37 common circadian genes identified in at least 8 tissues having tissue data, the median of variances of log₂(gene expression) across 61 tissues was 2.28. In comparison, the expected median of variances of log₂(gene expression) for the same number of randomly selected genes was 0.54 based on 10⁶ random simulations. The correlation coefficients r_ij between the tissue gene expression profiles of the common circadian gene pairs (i,j) were negatively correlated with their circadian phase differences d_ij (r = −0.22, p<10⁻⁸). To further demonstrate the relationship between r_ij and d_ij, we defined two functions y ₊(x) = median(d_ij(r_ij>x)) and y ₋(x) = median(d_ij(r_ij<x)) for −1 ≤ x ≤ 1. We plotted y ₊(x) and y ₋(x) in Figure S1. y ₊(x) for x>0 is significantly lower than the median of d_ij for all gene pairs (5.84) and reaches the minimum 3.068 at x = 0.64, whereas y ₋(x) for x<0 is significantly higher and reaches the maximum 9.939 at x = −0.26. These results indicated that the common circadian genes with positive correlations in their tissue gene expression profiles tended to have closer circadian phases, whereas those with negative correlations in tissue gene expression profiles tended to have larger differences in their circadian phases.

Comparison between Tissues and Species

We used the median of phase differences of circadian oscillating genes shared by two tissues as the distance measure of global phase dissimilarity between two tissues. We use these distances to cluster the phases of circadian oscillating genes in all 21 datasets using hierarchical clustering with complete linkage (Figure 2). For the mouse tissues where multiple datasets were available, i.e., liver, SCN, whole brain, and heart (whole heart, atria, and ventricle), we conducted pair-wise comparisons of the phases across tissues, using the “circular ANOVA” function for the genes identified as circadian oscillating in at least two datasets in each tissue under comparison. The same method was used to compare mouse liver data with rat liver data. To compare the circadian oscillating genes across species, rat, macaque, and human gene symbols were converted to mouse orthologs using the HomoloGene database of NCBI (build 56, http://www.ncbi.nlm.nih.gov/HomoloGene).

Gene Ontology Analysis

Gene symbols of circadian oscillating genes identified in each tissue in mouse, rat, macaque, and human were uploaded to Gominer [9] for Gene Ontology (GO) annotation and enrichment analysis. We selected the biological processes significantly over-represented in circadian oscillating genes in each tissue using False Discovery Rate (FDR) less than 0.05 as the criterion. For the circadian oscillating genes in each enriched biological process, we further tested their associations with any specific phase intervals using the Fisher's exact test with a rotating window method. In each 1,000 equally spaced phase intervals of size 4 hours between CT0 and CT24, the Fisher's test was applied to test the association between the biological process and the phase interval. The smallest p-value among Fisher's tests in all intervals was obtained to represent the significance of the association. The significant biological processes (p<0.005) in each tissue were colored using a color circle to represent their associated circadian phases. We visualized the significant biological processes as GO maps created by Cytoscape program (version 2.5). The significant biological processes were represented by the nodes and their hierarchical GO relationship was represented by the directed edges between them so that close-related biological processes were clustered together. All GO maps in different tissues can be found in our website (http://www.picb.ac.cn/circadian/). We manually selected the most representative biological processes for each GO cluster and summarized the result in Table S3.

Promoter Analysis

Transcriptional start sites (TSSs) information of mouse and human were integrated from three databases: DataBase of Transcriptional Start Site (DBTSS) [13],[14], the CAGE (Cap-Analysis Gene Expression) database of Fantom3 (Functional annotation of mouse) project [15], and the NCBI RefSeq database [16]. The criteria to select the TSSs were as follow: for DBTSS TSSs, the proportion of confident cDNA clones (non-exonic start clones, i.e., the clones mapped to the non-exonic regions of the genome) was not less than 0.75; for CAGE TSSs, the total number of corresponding CAGE tags was not less than 2 and can be mapped around the 5′ end of a known mRNA. If no TSS can be found for the gene from either DBTSS or CAGE under the above criteria, the 5′ end of the mRNA in RefSeq (human build 36.1, mouse build 36) was used as the TSS of the gene. Human (hg18 or NCBI build36.1) and mouse (mm8 or NCBI build 36) genome sequences were downloaded from UCSC. The 3000bp flanking sequences of each TSS were extracted from the genome as the promoter regions. As the CAGE database was based upon the older versions of human and mouse genome, i.e., hg17 and mm5, we mapped the CAGE TSSs to the new version of genomes using liftOver program in UCSC. In addition, orthologous promoter regions of mouse (mm8) vs. human (hg18) genome alignment results were also downloaded from UCSC.

Three positional weight matrix (PWM) based motif searching programs, match [17], motifscan [18], and profilestas [19], were used to identify the putative transcriptional factor binding sites (TFBS) on the extracted promoter regions. All vertebrate PWMs in TRANSFAC 11.2 were used as inputs in these programs. For the match program, we used the cut-off profile that minimizes the false positive rate, i.e., minFP profile in TRANSFC 11.2. For the motifscan program, we used a third-order background model by the CreateBackgroundModel program [20] to distinguish between the motifs that occurred frequently throughout the genome and the ones that were specific to the promoter regions. For the profilestas program, we first used the profilestas package to generate the scoring matrix and scoring threshold that minimized the false positive rate for each PWM. Then we used the patser program [21] to scan the promoter sequences and select the TFBSs above the scoring thresholds. The putative TFBSs predicted from all three programs have been compared and yielded very similar results. For simplicity, all the promoter analysis results presented in this paper were based on the match program.

We first tested the significant over-representation of putative TFBSs among a total of 568 PWMs on the promoters of circadian oscillating genes using the Fisher's exact test (p<10⁻⁴), using the promoters of all known genes as the background. Among the significant PWMs, we again tested their associations with any specific circadian phase interval in each tissue using the Fisher's exact test with a rotating window method as described above in the GO analysis (p<0.005). To remove the redundancy in PWMs, we grouped the PWMs into TF families according to their classifications in the TRANSFAC database and we averaged the associated circadian phases of significant TF PWMs in the same TF families using the “mean” function in the R “circular” package. The results are summarized in Table S4. The detailed information about TF enrichment and their associations with any specific circadian phase intervals can be found in our website (http://www.picb.ac.cn/circadian/).

Knockout or Mutant Mouse Microarray Data

We collected microarray data in different tissues or cell types from knockout or mutant mice, including liver and skeletal muscle in a Clock mutant, atrium and ventricle in a cardiomyocyte-specific Clock mutant, liver in a liver-specific conditional Nr1d1 mutant, aorta in Arntl and Npas2 knockout or mutant, liver in a Rora/Rorc knockout, liver and kidney in a Dbp/Hlf/Tef knockout, liver in a Ppara-null mice on Sv129 background treated by the Ppara agonist Wy14643, NIH 3T3 cells under Cebpa/b/d/e transfection, S49 cells in a Pka knockout under cAMP stimulation, cortex and thymus in a Egr1/Egr3 knockout, liver and primary chrodrocytes in a Nr3c1 (glucocorticoid receptor) knockout treated by the glucocorticoid agonist deamethasone (DEX), and embryonic fibroblast in a Hsf1 knockout under heat shock (Table S5). We also included the microarray experiment in the SCN of mouse exposed to 30 minute light pulse at 1 hour after the light off period compared to a dark pulse [10]. For the knockout or mutant mice microarray data where time series were available, we applied a two-way ANOVA using genotypes and time series as factors. The p-values and fold changes in the genotype comparison were used. For the knockout or mutant mice microarray data where external treatments such as Wy14643, cAMP, DEX, and heat were available, we applied a two-way ANOVA using genotypes and treatments as factors. Here the p-values and fold changes in cross-interactions between two factors were used. For Dbp/Hlf/Tef and Egr1/Egr3 knockout or mutant and Cebpa/b/d/e transfection experiments, we applied one-way ANOVA using genotypes as factor. For Rora/Rorc, Arntl, and Npas2 knockout or mutant experiments, we applied the LIMMA program using genotypes as the factor. In the Rora/Rorc knockout or mutant experiment, Rora knockout, Rorc knockout, Rora/Rorc double knockout were treated as the same genotype. In Pka knockout or mutant experiment, only the data at 0 hr and 2 hr of cAMP stimulation were used to include the directly affected genes in the cAMP signaling cascade. In Dbp/Hlf/Tef knockout or mutant experiments, the averages of log₂(p-value) and log₂(fold change) in three experiments: triple knockout vs. wild type in liver, triple knockout vs. triple heterozygotes in liver, and triple knockout vs. wild type in kidney were used as the overall log₂(p-value) and log₂(fold change). Ppara knockout data were obtained from the third and fourth study in [22]. To combine the results in third and fourth studies, we extracted the probe sets with consistent log fold changes of Ppara knockout effect of both studies. The maximum of p-values of both studies and mean of log fold changes were used. For Egr1/Egr3 knockout in cortex and thymus and Nr3c1 knockout in liver and primary chrodrocytes, the significantly affected gene lists were simply merged in two tissues or cell types. In all knockout or mutant data, a p-value less than 0.01 and a |log₂(fold change)|>0.5 were used to identify the significantly up- or down-regulated genes in the knockout or mutant.

To reliably identify Arntl/Clock and Nr1d1/Rora/Rorc controlled genes, we combined the evidences from multiple datasets. Arntl/Clock controlled genes were identified as those satisfying two out of the five conditions: down-regulated in the Clock knockout in liver, down-regulated in the Clock knockout in skeletal muscle, down-regulated in the cardiomyocyte-specific Clock knockout in atria, down-regulated in the cardiomyocyte-specific Clock knockout in ventricle, and down-regulated in the Arntl knockout in aorta. As Nr1d1, a repressor, was significantly down-regulated in the Arntl or Clock knockout or mutant, the significant up-regulation in the Arntl or Clock knockout or mutant was also considered to be the evidence for Nr1d1 controlled genes. Thus, Nr1d1/Rora/Rorc controlled genes were identified as those satisfying one out of the seven conditions: up-regulated in the Clock knockout in liver, up-regulated in the Clock knockout in skeletal muscle, up-regulated in the cardiomyocyte-specific Clock knockout in atria, up-regulated in the cardiomyocyte-specific Clock knockout ventricle, up-regulated in the Arntl knockout in aorta, up-regulated in the Nr1d1 conditional knockout, and down-regulated in the Rora/Rorc knockout. Dbp/Hlf/Tef, Ppara, Egr1/Egr3, Pka, Nr3c1, and Hsf1 controlled genes were identified as those that were significantly down-regulated in a knockout or mutant mouse compared to the wild type mouse. CEBP controlled genes were identified as those that were significantly up-regulated in Cebpa/b/d/e transfected cells compared to the control cells.

We identified 380 Arntl/Clock, 1,166 Nr1d1/Rora/Rorc, 53 Npas2, 53 Dbp/Hlf/Tef, 627 Cebp, 536 Ppara, 710 Egr1/Egr3, 464 Pka, 341 Nr3c1, and 425 Hsf1 controlled genes from the knockout or mutant experiments. To identify the direct target genes of transcription factors in knockout or mutant experiments, we required that the significantly affected genes in a knockout or mutant must have at least one putative binding site of their respective transcription factors in the promoter regions. We identified 320 EBOX, 295 RRE (Rev-erb/Ror element), 47 Npas2-regulated element, 43 DBOX, 607 CEBP, 516 PPRE (peroxisome proliferator responsive element), 492 EGRE (Egr element), 455 CRE (cAMP response element), 326 GRE (Glucocorticoid response element), and 122 HSE (Heat shock element) directly controlled genes after combining with the promoter analysis (Table S6).