Mouse genotypes drive the liver and adrenal gland clocks

Circadian rhythms regulate a plethora of physiological processes. Perturbations of the rhythm can result in pathologies which are frequently studied in inbred mouse strains. We show that the genotype of mouse lines defines the circadian gene expression patterns. Expression of majority of core clock and output metabolic genes are phase delayed in the C56BL/6J line compared to 129S2 in the adrenal glands and the liver. Circadian amplitudes are generally higher in the 129S2 line. Experiments in dark – dark (DD) and light – dark conditions (LD), exome sequencing and data mining proposed that mouse lines differ in single nucleotide variants in the binding regions of clock related transcription factors in open chromatin regions. A possible mechanisms of differential circadian expression could be the entrainment and transmission of the light signal to peripheral organs. This is supported by the genotype effect in adrenal glands that is largest under LD, and by the high number of single nucleotide variants in the Receptor, Kinase and G-protein coupled receptor Panther molecular function categories. Different phenotypes of the two mouse lines and changed amino acid sequence of the Period 2 protein possibly contribute further to the observed differences in circadian gene expression.

The majority of organisms on Earth have evolved a robust inner body circadian (circa = approximately, dian = day) clock that ticks away in most of their cells. It controls physiological, behavioural and cellular processes, ranging from body temperature, sleep/wake and feeding rhythms to metabolism and hormone secretion 1 . In mammals, the hierarchical structure of the circadian system necessitates that the master clock located in the suprachiasmatic nucleus (SCN) of the hypothalamus not only synchronises peripheral clocks but also entrains its own clock to the day/night cycle 1 . While the SCN clock seems to entrain to the outside world primarily through photic information via the retina and the retinohypothalamic tract, peripheral tissues entrain to cues through direct (humoral and neural connections) and indirect SCN signals (rest/activity cycles, feeding time, body temperature) 2,3 . Despite differences in the anatomical structure of the circadian systems, the molecular mechanisms generating oscillations are identical in all cells and consist of a transcriptional-translational feedback loop. The positive part of this loop, composed of a BMAL1 dimer with either CLOCK or NPAS2, enables E-box mediated transcriptional activation of core clock and clock controlled genes (CCG). 43% of all protein coding genes show circadian rhythms in transcription somewhere in the body, largely in an organ-specific manner 4,5 .
Inbred mouse strains are invaluable tools towards understanding the molecular biology and influence of rhythm perturbations on disease development in vivo 2,6 Especially for mutant strains, conducting research on one mouse strain alone may lead to missing diverse phenotypes 7 , as had been shown for the development of metabolic syndrome in Clock mutant mice 8,9 . Other aspects of circadian physiology are also affected by genotype including free-running activity rhythms 10,11 , feeding cycles 12 and phase-shifting effects 13 . Among rare understood molecular effects are mutations in the Hiomt gene 14 that result in different melatonin production in mouse strains. Recent whole genome sequencing projects exposed a large number of structural and single nucleotide variants between the inbred mouse strains 15 . Mouse knock-out models of core clock genes were frequently derived by gene targeting of embryonic stem cells of 129 mouse strains 16,17,18 and the initial analyses performed on mice with a mixed genetic background, where results can be confounded by background effects 19 .
We discovered crucial changes in circadian gene expression in the C57BL/6JOlaHsd and the 129S2/ SvPas × C57BL/J6 strains in peripheral organs. Whole exome sequencing and in silico analyses were applied to expose the nucleotide variants which could contribute importantly to the observed differences in the circadian or diurnal patterns of genes expression (Fig. 1).

Results
Circadian gene expression in liver and adrenal glands of C57BL/6 and 129S2 mouse lines. We evaluated gene expression patterns of 10 core clock and 15 metabolic genes collected under DD and LD conditions (Supplementary Figs S1 and S2). In liver we screened expression of the genes from cholesterol synthesis (Cyp51, Cyp7a1, Cyp8b, Por, Hmgcr) and regulators of metabolism (Ppara, Ppargc1a, Car, Rxr). In adrenal glands we focused on steroid metabolism (Cyp51, Cyp17a1, Cyp11a1, Cyp11b2, Cyp21a1). A detailed list is in Supplementary Table S1. Time series gene expression data were analysed by cosinor analysis to determine the presence/absence of circadian oscillation (Table 1, Supplementary Table S2). Under LD conditions all core clock genes displayed a clear diurnal rhythm in liver and adrenal glands of both mouse lines (Table 1). Under DD, Dec2 and Cry2 in the liver and Per1 in adrenal glands are arrhythmic in C57BL/6 (Supplementary Table S2). Metabolic genes are significantly more affected by genotype since many are circadian (in DD) and diurnal (in LD) in one but not the other mouse line under the same condition and tissue (Supplementary Table S2). Overall, a higher number of investigated genes displays circadian or diurnal expression in the 129S2 genotype (87%) compared to the C57BL/6J (76%) ( Table 1). In Table 2 and Supplementary Table S3 we compared circadian phases and amplitudes of gene expression in both mouse lines. Based on mathematical models of circadian rhythms constructed from in vivo data we regarded phases as more important since changes in phases lead to large perturbations of the whole circadian system 20 . The comparative analysis shows that phases of several core clock genes are affected by the genotype, including Bmal1, Dec1 and Cry1 in the liver, and Per2, Dec1, Cry1 and Rev-ErbA in adrenal glands (Figs 2 and 3). Interestingly, whereas in liver the circadian phases of gene expression differ between mouse lines in complete darkness (DD) (Figs 2A and 3D), the differences in phases in adrenal glands are seen primarily under LD conditions (Figs 2B and 3B). As shown in Fig. 3, the peak expression of core clock genes is phase delayed in C57BL/6J compared to the 129S2 with the exception of Rev-ErbA in adrenal glands (Fig. 3A).
Generally, metabolic genes were more variable concerning the presence/absence of circadian expression which makes comparisons of some phases and amplitudes impossible. More importantly, in both light conditions and tissues all metabolic genes for which circadian or diurnal rhythmicity was confirmed showed a phase delay in the C57BL/6J genotype (Figs 2 and 3). The phase delay is on average 2.8 h for core clock genes and 6.1 h for metabolic genes for both tissues and light conditions. Similarly, as for the phases, genotype dependent rules were observed also for amplitudes. For the liver metabolic genes, the amplitudes were less affected by the genotype compared to core clock genes (Supplementary Table 3) while in adrenal glands the amplitudes were affected only when light had been present (LD). Generally, amplitudes were in most cases higher in the 129S2 genotype (Cry1 and Cyp11b1 are exceptions) ( Table 2,  Supplementary Table S3), suggesting differences in entrainment and in transmission of the light signal to the adrenal gland between C57BL/6J and 129S2 lines.  Table S4, S5). A total of 90,994 SNVs (85,224 in the 129S2 and 5,770 in the C57BL/6J lines) were discovered in exomes when compared to the NCBI reference mouse line (NCBIM37/mm9). The 3,426 SNVs that were identical in both mouse lines were removed from further analysis. As shown in Table 3 Table 1. Assessment of circadian rhythmicity in gene expression. The presence of circadian (DD) or diurnal (LD) oscillation for each gene was determined with the cosinor analysis of expression data gathered with qRT-PCR. A p-value of 0.01 was determined as a cut of value above which genes were not considered to have a circadian or diurnal expression pattern. The number and percentage of genes displaying a circadian or diurnal pattern is shown below for each of the 8 possible conditions based on the number of strains (2), tissues (2) and light conditions (2). A list of core clock and metabolic genes can be found in Supplementary Table S1.  Table 2. Comparison of amplitudes and phases (peak expression) of core clock and metabolic genes in 129S2 line and C57BL/6J. The amplitudes and phases were compared based on the cosinor analysis (a 95% confidence interval for amplitudes and phases was calculated for each gene). Genes that did not show a circadian or diurnal expression were not included in the analysis. A p-value of 0.01 was determined as a cut of value above which genes were not considered to be differentially expressed between stains. The number and percentage (%) of genes with different circadian or diurnal patterns for each condition and tissue are reported. within exon regions (Fig. 4A). These nucleotide variations can affect the resulting phenotype, especially 6,938 (36.5%) variants causing nonsynonymous substitutions including loss or gain of stop codons. These variants were distributed in 2,553 genes of which the majority (86%) carried multiple (up to 4) variants. The rest (356 genes) showed a higher number of nonsynonymous SNVs. To predict which of the variants could have deleterious effect, Grantham Matrix Score (GMS) was determined for each nonsynonymous SNV. The majority of nonsynonymous SNVs belong to conservative (38.31%) and moderately conservative groups (45.07%) while the most interesting are the 5% and 11.6% SNVs from the radical and moderately radical groups (Fig. 4B). 287 genes from the radical group were annotated: significant association was found with Panther molecular function categories including Receptor (p < 5.7e-4), Kinase (p < 1.9e-2) and G-protein coupled receptor (p < 2.4e-2) (Fig. 4C). DAVID Analysis of SNVs that result in gain or loss of stop codons also showed significant association with Panther molecular function categories: Receptor (p < 4.7e-2) and G-protein coupled receptor (p < 6.5e-2). 42 out of 80 variants discovered by exome sequencing were included in DAVID analysis. The remaining SNVs were removed after manual inspection with Variant Effect Predictor (VEP) since Ensembl had shown them not to be true gain or loss of stop codons variants (Supplementary Table S6). Since the biggest differences in gene expression was detected in adrenal glands in LD conditions, we screened the literature for known genes and pathways involved in transmission of light information from the retina to the adrenal gland. A list of over 200 genes was created including genes from light reception, SCN light transduction, SCN neuropeptides, steroid metabolism and regulation, etc. (Supplementary Table S7). These genes were analyzed for the presence of SNVs and 74 genes were identified. The most interesting were melanopsin (Opn4), a G-protein coupled receptor and a photopigment that is important for proper circadian entrainment 21 ; vasoactive intestinal peptide receptor (Vipr2) important for synchronization of pacemaker neurons in the SCN 22 ; Neuronal PAS domain-containing protein 2 (Npas2), a functional substitute of the CLOCK protein in the SCN 23 and the steroidogenic factor 1 (Sf1), essential for expression control of most steroidogenic genes 24 .
40 SNVs were discovered in intronic, UTR and exonic regions of Opn4 with 4 non-synonymous SNVs located in the C-terminal cytoplasmic region with several amino acids that can undergo phosphorylation 25,26 . Three SNVs (c.C1504T:p.P502S, c.G1496C:p.C499S and c.A1402T:p.T468S) induce a conversion of the original amino acid to serine in the 129S2 mouse line, increasing potential phosphorylation sites and possibly affecting signaling. In humans mutations of the Opn4 gene result in seasonal affective disorder 27 .
Vipr2 gene contains 25 intronic and 1 nonsynonymous exonic SNV (c.G1309C:p.A397P). The later induces conversion of alanine to proline in C terminal region, potentially changing the secondary structure. Vipr2 −/− mice SNV location/type   lack a normal circadian rest/activity behavior and have no circadian expression of the core clock genes Per1, Per2, and Cry1 in the SCN. They also failed to show acute induction of Per1 and 2 by nocturnal illumination 28 . 20 SNVs were located in Npas2, one in the 5' UTR and 4 in exons, two of these nonsynonymous (c.T1938C:p. L481P, c.A2189G:p.T562A). Npas2 is expressed in neural tissues and is critical for adaptation to restricted feeding 29 .
Sf-1 (Nr5a1) contains 2 SNVs, one is nonsynonymous (c.G715T:p.A172S) and lies near the phosphorylation site (p.203) which is important for maximal SF1-mediated transcription 30 . The newly formed serine can also be phosphorylated, suggesting that this SNV could affect the protein activity in the 129S2 mouse line. Sf-1 expression requires E-box element located in the promoter region 31 . SF-1 regulates Cyp11a1 and Cyp21a1 expression that show different amplitudes in LD in adrenal glands (Supplementary Tables S2 and S3). Described SNVs and others mentioned in Discussion could affect the gene expression in adrenal gland.
Differences in the liver expression patterns could be influenced also by metabolic signals (majority of differences between strains were detected in DD). Literature screen and SNV analysis identified some candidate genes that could potentially serve as sensors of hepatocyte external stimuli, but only Ppargc1a and Npas2 harbor nonsynonymous SNVs in coding regions. Ppargc1a encodes a transcriptional coactivator (PGC1a) that coordinates gene expression of mitochondrial biogenesis, respiration, hepatic gluconeogenesis and is considered as one of the most important integrators of external stimuli. SNV analysis revealed a change in amino acid sequence (c.G2163A:p.R675H), that is important for interaction with other transcription factors.
Even if both mouse lines had the access to the same food ad libitum, the phenotypic differences in food and water consumption could participate to hepatic circadian expression differences 32 . The hepatic clock can easily be re-set with the restricted feeding (RF). According to the literature, one of the most described signaling pathways that could trigger hepatic re-setting is the insulin signalling cascade 33 . We chose 50 genes involved in the insulin signalling 34 and checked for the presence of SNVs (Supplementary Table S8). The insulin receptor gene Insr contains only synonymous SNVs (6 exonic, 1 intergenic, 4 in 3UTR, 14 intronic) so it is difficult to speculate on their function. A nonsynonymous SNV (c.T590G:p.V197G) was found in Irs2 (insulin receptor substrate protein) in the PTBi domain which is important for the phosphorylation status dependent peptide binding. Since the valine to glycine change is near the phosphoinositide binding site the activity of the IRS2 might be modified. Non-synonymous SNVs reside in other feeding related genes: in protein tyrosine phosphatase gene Ptpn14 which is involved in mediation of insulin sensitivity by de-phosphorylating tyrosine residues in the insulin receptor and reducing its activity. Among 21 SNVs one is nonsynonymous (c.A2660G.p.H887R) in the protein region with yet unknown function. SNVs are found also in other genes: in the gene c-Jun-N-terminal kinase Mapk10 (c.C1370T:p.A457V) at the site with the unknown function, in Cdc42bpa (c.G12562A:p.R521Q), and in Cdc42se2 (c.C128T:p.S43F) near the GTPase interaction site. CDC42 binds to p85a PI3K isoform and modulates insulin sensitivity of the cell. Nine SNVs are found in Cap1 gene involved in the CAP-Cbl pathway collaborating with the PI3K pathway in stimulation of GLUT4 translocation. One SNV (c.T724C:p.S242P) is located in the region with the trypsin-like protease activity.
At the molecular level, the circadian system is maintained by intertwining feedback loops, so we wondered whether the SNVs from the exome sequencing lie in our genes of interest and could therefore affect the differential circadian patterns in mouse lines. 77 SNVs were discovered in core clock (Bmal1, Per1, Per2, Dbp) and metabolic genes (Cyp21a1, Cyp17a1, Cyp51, Ppargc1a, Car and Ppara) ( Table 4 and Supplementary Table S9) with more than 75% in Bmal1, Per2, Ppargc1a and Car (in average 14 SNVs per gene). Nonsynonymous mutations reside in the core clock gene Per2 and in transcriptional coactivator Ppargc1, with GMS scores of 58 and 29, respectively. In Per2 serine is changed to threonine at position 172 (G515C:p.S172T) just before the PAS 1 domain important for the dimerization. The amino acid substitutions for the Ppargc1a gene are described above.
Genome variations in binding sites of core clock proteins. We evaluated promoters of genes of interest and checked for SNPs of clock protein binding sites. DNA binding sites of BMAL1, CRY1, CRY2, PER1, PER2, CLOCK and NPAS2 were taken from ChIP-Seq experiments performed by Koike et al. 35 . Since the majority of binding regions were outside the exome sequencing area, we applied data from the Imputed Mouse SNP resource 36 that included six 129 lines: 129P1/ReJ, 129P3/J, 129S1SvlmJ, 129S6, 129T2/SvEmsJ and 129 × 1/SvJ, for which SNPs are determined compared to the C57Bl/6J reference genome. We searched for the presence of SNPs within clock protein binding regions (Supplementary Table S10). Figure 5A  When determining the density of SNPs within each clock protein binding region, we compared it with the location of ChIP-Seq peaks (Fig. 5B). In most cases the location of ChIP-Seq peak coincides with a peak in SNP density so the differential regulation of circadian expression is expected.
To evaluate the SNP density in transcription factor binding regions compared to the rest of the DNA, we assessed it in the context of the DNase hypersensitive sites (DHS) form open chromatin regions 37 (Supplementary  Table S11). DHS regions where clock related proteins bind contain one SNP per 219.6 nucleotides in comparison to all standard DHS sites which contain one SNP per 224.8 nucleotides on average; the difference is statistically significant (p-value 0.045).

Discussion
The relatively high degree of anatomic, physiological and genomic similarities between humans and mice (approx. 95% of genes are shared) have made mice important models for the study of various biomedical issues 6,38 . To Scientific RepoRts | 6:31955 | DOI: 10.1038/srep31955 eliminate misleading results due to genetic variability of mice, lines with the lowest degree of intra-strain variability are selected. This precluded us from seeing the effects of mouse genetic backgrounds on studied processes 7 .

Symbol
Entrez ID  Table 4. SNVs in genes with genotype related expression differences. The table shows the number of SNVs that were discovered between the C57BL/6J and 129S2 strain by exome sequencing in genes whose circadian or diurnal expression profiles were affected by genotype. Here we show a clear genotype effect on circadian or diurnal gene expression, with differences in phases and amplitudes of core clock and metabolic genes in C57BL/6J and 129S2 mouse lines. In liver of the C57BL/6J genotype two core clock (Cry2 and Dec2) and several metabolic genes (Car, Cyp7a1, Pgc1a and Pparα) lacked circadian expression under DD in contrast to the 129S2 genotype (Supplementary Table S2). While the absence of circadian expression in the two clock genes is not sufficient to abolish rhythmicity due to the redundancy of the clocks molecular network 39 , as shown by mouse knock-out models of Cry2 16 and Dec2 40 , it could lead to changes in circadian expression of clock controlled genes 41 . This is supported by the absence of circadian rhythm in two important regulators of hepatic energy metabolism, Pparα and Pgc1α. Although expression of both genes is regulated by the clock, they also reciprocally regulate expression of Bmal1 and in turn circadian rhythms (Fig. 6A) 42,43 . The lack of circadian peak of Pparα and Pgc1α in the C57BL/6J genotype could lead to reduced and delayed peak of Bmal1 expression (Supplementary Figure S1, Supplementary Tables S2 and S3). This might also be the reason for the reduced amplitudes of Per1, Cry1, Cry2, Dbp and Rev-ErbA as well as the observed delayed peaks of Dec1 and Cry1 (Figs 2A, 3D and 6A, Supplementary Table S3). PPARα is a nuclear receptor that binds exogenous as well as endogenous substances. The availability of those ligands can alter the PPARα activity which is not necessary following the mRNA expression. Fatty acids (arachidonic, polyunsaturated FAs and their derivates) are among the most important ligands of PPARα 44 . In our experimental setup both mouse lines were kept under the same laboratory conditions and had the access to the same amount and type of food. However, the ingestion, absorption and metabolism of ligands can differ due to genetic background. Similar interpretation could be used also for the activity of PGC-1α which is a master integrator of external stimuli: AMPK phosphorylation and by SIRT1 regulated epigenetic status can alter PGC-1α activity according to the metabolic state of the organism (energy needs, cold, fasting, exercise etc) 45 .

Number of Mutations
The largest influence of genotype was, however, observed in adrenal glands under LD conditions. Here the majority of measured clock genes were differentially expressed between mouse lines either in amplitude or phase (Figs 2B and 3A, Table 2). This large genotype effect may be due to entrainment of the circadian network in peripheral organs. The main pacemaker, the suprachiasmatic nucleus (SCN) of the hypothalamus, entrains peripheral tissues either directly through humoral and neural factors or indirectly through behaviour (locomotor activity, fasting/eating) 46 . In LD conditions the daily entrainment takes place with the onset of light. It was shown that adrenal glands are in the first line of peripheral organs that respond to light signals 47 . Several layers of signal Figure 6. The effect of genotype on gene expression. Circadian (DD) or diurnal (LD) expression of core clock and metabolic genes was shown to be affected by mouse genetic background. A mechanism based on gene expression data only that could explain the differences observed in liver is proposed under A. (A) The noncircadian expression of Ppargc1a and Ppara in the C57BL/6J strain could lead to a reduced and phase delayed expression of Bmal1 (confirmed by expression data) which in turn could lead to a reduced and delay expression of other core clock (confirmed for Cry1 and Bhlhe40). (B) In adrenal glands the largest differences observed were under LD conditions, leading us to believe that entrainment pathways could be affected by mouse genotype. Several SNVs were discovered in genes known to be involved in circadian entrainment such as Opn4, Vipr2 and Npas2. The picture depicts several neural pathways between different brain regions known to be involved in entrainment as well as hormonal connections important in adrenal regulation that could be affected by SNVs present between strains. transduction and regulation exists between the SCN and the adrenal gland (the hypothalamus-pituitary-adrenal -HPA axis and neural innervation 46,21 ).
SNVs found by exome sequencing in genes with differential circadian or diurnal expression could not fully explain the observed expression differences. Functional characterization of genes with highest variability between the 129S2 and C57BL/6J mouse lines, followed by enrichment analysis of genes involved in adrenal and liver circadian entrainment, returned genes from Receptor, Kinase and G-protein coupled receptors functional categories. These are key pathways of circadian photoreception in the SCN and light signal transduction to the periphery 48 (Fig. 6B). In addition to the genes that have evident role in circadian events, such as Opn4, Vipr2 and Npas2 described in detail in the exome sequencing results section, other genes with nonsynonymous SNVs are: neural transmitters transporters (Slc6a4), receptors (Chrnb4, Nmur2), signal proteins (Ssfa2), proteins involved in cholesterol to steroid hormone synthesis (Tspo, Stard13) and cholesterol to bile acid conversion enzyme (Cyp7a1). 26-31% of nonsynonymous SNPs are expected to have effect on protein function 49 , thus it is possible that some of these SNVs could influence entrainment pathways leading from SCN to the adrenal gland and result in subtle changes in circadian gene expression. One of the mouse lines used in our experiments was C57BL/6JOlaHsd. It is generally known that a related subline C57BL/N carries a mutation in Crb1 gene known as "rd8" that causes retinal lesions. Zurita et al. 50 performed SNP analysis of ten different C57BL/6 strains and found no "rd8" mutation in C57BL/6JOlaHsd subline. Therefore this mutation cannot contribute to different light perception.
Adrenal gland and liver are both governed by the SCN clock and are innervated by the splanchnic nerve, but both organs differ greatly regarding the signals for entrainment. It was shown before that the dominant factor for re-entrainment differs importantly among organs and tissues. For the liver the restricted feeding overruns HPA-axis derived glucocorticoid signals 51 . Adrenal gland is from this point of view more strongly affected by the SCN mediated light perception. Under LD conditions organism is re-entrained daily with the onset of light, so it is expected that differences in adrenal gland LD expression will be more profound than differences in the liver. Ishida et al. 47 clearly showed that the re-entrainment induced by light affected adrenal gland Per1 expression followed by a corticosterone peak, but did not affect the liver Per1 expression. In our experimental set-up both mouse lines were housed and fed under the same conditions so the difference in the liver daily re-entrainment (in LD conditions) are not expected unless some great difference in processing of the metabolic re-entrainment signals is present due to the genotype.
The mouse lines 129S2 and C57BL/6 differ in several phenotype facts. For example, the 129 line has a lower body weight, higher daily food intake and lower water intake compared to C57BL/6 32 . Additionally, restricted feeding, glucocorticoid levels, or the temperature cycles, can all influence the circadian entrainment of peripheral organs. These phenotypic differences, even if mechanistically not well defined, could also affect metabolic regulation of the hepatic circadian rhythmicity. The insulin signalling is among better described mechanisms of hepatic clock synchronization. We found several SNVs in insulin -related genes but the effectiveness of these SNVs awaits determination. Berglund et al. 52 indeed showed that 129 mice have the greatest insulin secretion and C57Bl/6J respond the fastest to hyperinsulinemia after 5h fasting. Among 4 mouse strains used in their study C57BL/6J has the highest blood glucose and glucagon levels while 129 the lowest blood glucose and insulin and highest insulin clearance.
It is not negligible that 129/S2 and C57BL/6J strains differ in the average length of free-running period (23.93 ± 0.07 h and 23.77 ± 0.02 h, respectively) 10 . They differ as well in locomotor activity with C57BL/6J mice being more active than the 129 line 55 . The behavioural activity can change the circadian clock by suppression of SCN neuronal activity 56,57 , so such mechanism can represent a direct influence of phenotypic differences in the circadian clock. The body temperature rhythm is lower in C57BL/6 compared to three other mouse lines 58 . Buhr et al. 59 suggested that temperature is a universal resetting cue. During the restricted feeding regime, the body temperature cycles are slightly changed and this can activate the Heat-shock-protein pathway 60 . Heat shock factor 1 (HSF1) activity is strongly circadian and phosphorylated form can bind the HSE elements in Per2 promoter. HSF1 off-sync activation can cause Per2 expression in an immediate early fashion which can re-set the clock in peripheral organs [61][62][63] . We checked if some SNVs are present in the genes important for the activation of HSP pathway: in Hsp90aa1, Hsp90ab1 and Hsf1 no SNVs were found. We checked also the two HSE regions of the Per2 promoter 63 and also found no SNVs between the C57Bl/6 and 129/S2 lines.
A large contributor to the genotype effects could also lie in SNPs in the core clock protein binding sites. Mutations in binding sites determined by ChIP-Seq can largely reduce the circadian amplitude of gene expression 53 . However, since the density of SNPs is enriched at the ChIP-Seq peaks of clock protein binding regions (Fig. 5B) it is expected that these SNPs indeed have an important role in regulating the circadian expression. Changes at the genomic level alone will not be able to explain all differences observed. As we have shown in our previous work, the circadian amplitude of Cyp17a1 in adrenal glands is affected by methylation of its promoter region 54 .
In conclusion, our data show the effect of mouse genetic background on circadian expression of core clock and metabolic genes. The genotype influenced the tissue (liver and adrenal glands) and light dependent circadian (DD) or diurnal (LD) expression. While we uncovered potential SNVs that could explain these differences between the C57BL/6 and 129S2 mouse lines, confirmation of their causative effect is difficult due to complexity of the circadian system 46 . It is important to comprehend that the C57BL/6 and 129S2 mouse lines vary in tissue circadian expression patterns where the phases in C57BL/6 are generally delayed and the amplitudes lower. The experiments together with sequencing and data mining proposed the differences of clock-related transcription factors binding regions and the entrainment mechanisms of the peripheral organs as major causes for the observed circadian or diurnal differences. Methods Animal experiments, Circadian Collection of Samples and Isolation of materials. C57BL/6JOlaHsd and 129S2/SvPasCrlf × C57BL/6JRj male mice with free access to food and water were housed under a 12:12 h LD or under DD and sacrificed every 4h over a 24h period. Liver and adrenal glands were collected. RNA and DNA were isolated by standard protocols. Animal experiments were carried out in accordance with approved guidelines: European Convention for the protection of vertebrate animals used for experimental and other scientific purposes (ETS 123); National Institutes of Health guidelines for work with laboratory animals. All experimental protocols were approved by the Veterinary Administration of the Republic of Slovenia (license number 34401-9/2008/4, 34401-38/2009/2, 34401-44/2009/2).
Gene Expression Data Analysis. RT-q PCR was performed by gene-specific primers applying the SYBR Green method and data normalized as described earlier 64 . Normalized gene expression values were analysed by the Cosinor analysis 65 , calculated and visualized in program language R 66 .
Whole Exome Sequencing and Data Analysis. Whole exome sequencing was performed on the Illumina Hiseq2000 platform. The Burrows-Wheeler Aligner was used to align the mouse strains sequences to the NCBI reference sequence (NCBI37/mm9). The final BAM files were used as input for SOAPsnp in order to identify single nucleotide variations (Supplementary Data S1 and Supplementary Data S2). In order to predict which of the variants could have the most deleterious effect on the protein activity the Grantham Matrix Score (GMS) was determined for each nonsynonymous SNV (Fig. 4B) 67 and DAVID applied to annotate genes belonging to the radical group. Other procedures were performed in R and Bioconductor. Whole exome sequencing data is available at the European Nucleotide Archive under the Study code: PRJEB11861 (http://www.ebi.ac.uk/ena/ data/view/PRJEB11861).
In silico analysis. Genomic location of genes whose circadian or diurnal expression was measured by RT-qPCR was extracted from ChIP-seq data obtained from Koike et al. 35 . After obtaining genomic locations of binding sites the presence of SNPs was determined with the aid of the Imputed Mouse SNP Resource 36 (Supplementary Table S10). The NCBI Build 37 (mm9) was used as a reference for genomic locations for both databases.
Imputed Mouse SNP Resource 36 was used to count the number of SNPs in DHS sites from Ling et al. 37 supplemental Table S5A containing a merged list of Standard DHS sites. SNP density was obtained by dividing the number of SNPs by DHS site length. DHS sites of clock related protein were identified using Koike et al. 35 Supplemental Table S2 containing master peak list for each transcription factor. Average SNP density for DHS regions was calculated for each clock related protein individually (i.e. BMAL1, CLOCK, NPAS2, PER1, PER2, CRY1, CRY2) and compared to an average SNP density over all standard DHS sites 37 using a single-sample t test (Supplementary Table S11).