key: cord-0011776-lfm34fsw
authors: Li, Yan; Shan, Yongli; Kilaru, Gokhul Krishna; Berto, Stefano; Wang, Guang-Zhong; Cox, Kimberly H; Yoo, Seung-Hee; Yang, Shuzhang; Konopka, Genevieve; Takahashi, Joseph S
title: Epigenetic inheritance of circadian period in clonal cells
date: 2020-05-27
journal: nan
DOI: 10.7554/elife.54186
sha: 9e88d65d7ea23b4036ea7ee0f9ff9d1cca44cf38
doc_id: 11776
cord_uid: lfm34fsw

Circadian oscillations are generated via transcriptional-translational negative feedback loops. However, individual cells from fibroblast cell lines have heterogeneous rhythms, oscillating independently and with different period lengths. Here we showed that heterogeneity in circadian period is heritable and used a multi-omics approach to investigate underlying mechanisms. By examining large-scale phenotype-associated gene expression profiles in hundreds of mouse clonal cell lines, we identified and validated multiple novel candidate genes involved in circadian period determination in the absence of significant genomic variants. We also discovered differentially co-expressed gene networks that were functionally associated with period length. We further demonstrated that global differential DNA methylation bidirectionally regulated these same gene networks. Interestingly, we found that depletion of DNMT1 and DNMT3A had opposite effects on circadian period, suggesting non-redundant roles in circadian gene regulation. Together, our findings identify novel gene candidates involved in periodicity, and reveal DNA methylation as an important regulator of circadian periodicity.

Circadian oscillations maintain daily rhythms to control multiple physiological and behavioral processes, including metabolism, cell growth, immune response, and the sleep-wake cycle. Disruptions of the circadian clock have been linked with various disease processes and aging (Takahashi et al., 2008; Kondratova and Kondratov, 2012) . Circadian oscillations display remarkable fidelity in their periodicity even in the absence of environmental cues. This precision of the internal biological clock arises from a complex gene network. In mammals, the core of this network is composed of an autoregulatory transcriptional negative feedback loop involving Clock, Bmal1, Per1/Per2, and Cry1/Cry2, and there are additional feedback loops interlocked with the core (Takahashi et al., 2008; Mohawk et al., 2012; Takahashi, 2017) . Interestingly, although the cell-autonomous clock is ubiquitous, individual cells often do not maintain a perfect 24 hr circadian period, and within cell populations there are heterogeneous autonomous oscillations with a broad distribution of period length (Nagoshi et al., 2004; Welsh et al., 2004; Leise et al., 2012) . The heterogeneity in intrinsic period of hypothalamic suprachiasmatic nucleus (SCN) neurons confers important functions of phase liability and phase plasticity (Welsh et al., 1995; Liu et al., 1997; Ko et al., 2010; Mohawk et al., 2012) . However, it is still unclear how heterogeneous circadian periodicity is established and maintained under physiological conditions, or how much of this heterogeneity is heritable.

The origin of heterogeneity is complex, but may be driven by genetic variation, epigenetic modifications, and/or transcriptional noise (Jaenisch and Bird, 2003; Raser and O'Shea, 2005; Raj and van Oudenaarden, 2008; Burrell et al., 2013; Kelsey et al., 2017; Cavalli and Heard, 2019; Liu et al., 2019) . We have recently shown that nonheritable noise is the predominant source of intercellular variation in circadian period within clonal cell lines (Li et al., 2020) . However, it is still unclear what heritable factors contribute to period variation among different clonal cells. DNA methylation has been recognized as a chief contributor to gene expression states, and it is essential for mammalian embryonic development, with genome-wide methylation patterns changing during differentiation (Greenberg and Bourc'his, 2019) . There are three canonical cytosine-5 DNA methyltransferases that catalyze the addition of methylation marks. DNMT3A and 3B, the de novo methyltransferases, set up DNA methylation patterns during early development. Once established, DNMT1 will copy those patterns onto the daughter strand during DNA replication ensuring methylation maintenance (Jaenisch and Bird, 2003) . DNMT dysfunction has been associated with various diseases, and DNMT-deficient mice exhibit embryonic lethality (Greenberg and Bourc'his, 2019) . Numerous studies have supported the role of DNA methylation in gene silencing; however, more recent work suggests that DNA methylation can also be involved in transcriptional activation (Rinaldi et al., 2016; Yin et al., 2017b; Harris et al., 2018; Lyko, 2018) . Interestingly, despite high fidelity in mitotic inheritance, DNA methylation is variable across individuals, tissues, and cell types (Jaenisch and Bird, 2003; Jones, 2012; Varley et al., 2013) . Thus, we hypothesized that differential DNA methylation could contribute as a heritable factor underlying heterogeneous circadian oscillations in clonal cell lines.

Here, by examining phenotype-associated high-throughput multi-omics profiles in clonal cell populations, we identified and validated a pool of novel candidate genes regulating circadian period length and uncovered complex gene co-expression networks highly enriched in stress response and metabolic pathways. We next explored the origins of heterogeneous gene expression and found differences in global DNA methylation patterns that were associated with both silencing and activation of differentially expressed genes. Using gene knockdown studies, we also found that DNMT1 and DNMT3A have opposite effects on period length. Together, our findings demonstrate the important role of DNA methylation in the regulation of circadian period.

To assess cellular phenotypic heterogeneity, we utilized an immortalized mouse ear fibroblast cell line carrying a PER2::LUCsv bioluminescence reporter generated from Per2::lucSV knockin mice (Chen et al., 2012; Yoo et al., 2017) . We recently showed that these cells express persistent, robust, and cell-autonomous circadian oscillations over a 2 week period. Moreover, clonal cell lines generated from the parent culture had period distributions similar to those seen with single cells, indicating that circadian period is a heritable phenotype ( Figure 1A Li et al., 2020) . Here, we used the clonal cell lines to address the underlying molecular mechanism for heterogeneous circadian periodicity. To examine the stability of this heritability, twenty clonal cell lines were randomly selected and cultured continuously for 20 passages and tested for circadian period every five passages. Although two-way ANOVA revealed significant effects (p<0.01) of both cell line and passage, there was no interaction (p=0.09). Moreover, cell line was the dominant source of variation (74.70%), while passage only contributed 2.64% of the total variation. Multiple comparisons within each clonal cell line across passages identified a significant difference (adjusted p<0.05) for only~5% of comparisons (11 out of 200), which is consistent with 5% false positive rate. These results indicate that circadian period of clonal cell lines is stable and transmissible for at least 20 cell passages ( Figure 1C ).

To explore potential underlying mechanisms, we selected two groups of clonal cell lines from the two tails of the period distribution ( Table 1 , 5 short period (SP) and five long period (LP) clones) (Li et al., 2020) and performed RNA-seq analysis (Figure 2-source data 1). We compared their transcriptomic profiles and identified 5,137 period-correlated differentially expressed (DE) genes, with 2,782 genes upregulated and 2,355 genes downregulated in the LP group (Figure 2A To narrow down the target pool further and identify candidate genes more directly responsible for periodicity differences, we selected four additional groups of subclones established from two representative clonal cell lines with different periods: a shorter period subgroup and a longer period subgroup from short period clone#33 (SSP and LSP), or long period clone#114 (SLP and LLP), respectively ( Figure 2B ; Li et al., 2020) . These subclones and the original 10 clonal cell lines constituted a continuous period spectrum beneficial for identifying period-correlated genes ( Figure 2C , Table 1) .

We identified 535 additional period-correlated DE genes from subclones originating from SP clone#33 and 1,352 additional DE genes from subclones originating from LP clone#114 ( Figure 2D -E, Figure 2 -source data 1). By comparing the three RNA-seq datasets, 67 overlapping DE genes were identified ( Figure 2F ). From these, we selected 14 genes based on the strength of the correlation between their expression and circadian period length from all 88 samples and performed knockdown experiments to validate their function in circadian periodicity. Out of 7 positively correlated DE genes, knockdown of Ak3 and Trim3 significantly shortened period, whereas knockdown of Cpeb1, Lrrfip1, Rbfa, and Dars lengthened period ( Figure 3A -C, Figure 3 -source data 1). Out of 7 negatively correlated DE genes, knockdown of Ipo13 and Tmem165 significantly lengthened period, whereas Slc8a3, Jun, Med23, and Cpa4 knockdown shortened period ( Figure 3D -F, Figure 3-source data 1). Knockdown of two other genes, Eif4e2 and Rfx5, did not alter period length. We also examined the effect of knockdown of five representative genes in 10 clonal cell lines and found that they all showed the same period alterations as that seen in the parent culture . These results suggest that multiple genes function together to determine circadian period length and that there were no unique (clone-specific) effects on the direction (long or short) of the period changes. Since the majority of the DE genes identified here have never been reported as having effects on circadian period, these data provide a new pool of candidate genes functioning in circadian periodicity.

Large-scale gene networks are associated with period heterogeneity Because functionally related genes are usually co-expressed (Heyer et al., 1999) , we further characterized the period-correlated DE genes by examining their co-expression patterns. Using weighted gene co-expression network analysis (WGCNA), we generated 31 modules from the 10 clonal cell lines RNA-seq data ( Figure 4A , Figure 4 -source data 1). Several modules exhibited significant enrichment for period-correlated DE genes. Blue, lightgreen, green and darkred modules were enriched for positively correlated DE genes, while salmon, pink, red, and darkgreen modules were enriched for negatively correlated DE genes ( Figure 4B ). Ingenuity pathway analysis (IPA) revealed stress response signaling pathways and metabolic pathways were associated with the period-correlated DE genes, suggesting their important roles in circadian periodicity ( Figure 4C , Figure 4 -source data 1). IPA analysis of the correlated modules also revealed overlapping functional pathways. For example, the blue module is highly enriched for DE genes, and is also enriched for the EIF2 signaling pathway, which has been recently shown to regulate circadian period, consistent with the predicted elevated translational activity in LP group (Pathak et al., 2019;  Figure 4C ). To validate these results further, we used two different small molecules to activate the EIF2 signaling pathway in parent culture and observed significantly shortened period, consistent with what has been previously reported (Pathak et al., 2019;  Figure 4D ). In addition, the darkred module was enriched for the mTOR signaling pathway; the green and salmon modules were enriched for the protein ubiquitination pathway; and the pink module was enriched for NRF2-mediated oxidative stress response pathway (Figure 4 -source data 1). Interestingly, all three of these pathways have been shown to be functional in circadian periodicity, further confirming the functional importance of the co-expressed gene networks (Stojkovic et al., 2014; Ramanathan et al., 2018; Wible et al., 2018) . Further analysis of Protein-protein Interactions (PPI) revealed that co-expressed DE genes were also physically interconnected. For example, within the blue module there were several different tightly linked clusters, including those enriched for ribosomal RNA processing, protein ubiquitination, nucleotide and amino acid metabolism, and mRNA splicing, emphasizing the blue module as a transcriptional/translational related gene network ( Figure 4E ). Taken together, our results suggest that period heterogeneity is regulated by changes in large-scale functional gene co-expression networks.

To explore whether there was a genetic basis for heterogeneous gene expression, we performed whole-exome sequencing on SP clone#33 and LP clone#114. Interestingly, only four annotated genes carrying unique variants were identified (Supplementary file 1), but 2 of them are not expressed (Figure 2 -source data 1), and none of them have known circadian functions, suggesting that somatic mutations are unlikely to underlie the heterogeneous period distributions.

Cell-to-cell variability is also partially heritable via epigenetic modifications such as DNA methylation (Jaenisch and Bird, 2003; Jones, 2012) . To assess the contribution of DNA methylation in heterogeneous circadian periodicity, we used reduced representation bisulfite sequencing (RRBS) to explore DNA methylation profiles and their correlation with the period-correlated transcriptomes. Using 1,000 bp tiling windows genome-wide, we identified 16,520 significant differentially methylated regions (DMRs). Importantly, none of the core clock genes, even the few that were differentially expressed in the parental lines, had coding mutations or differential DNA methylation, except for a small DMR spanning~10 nucleotides located in exon 1 of Per1 (Table 2) . Of the DMRs found, 62% (10,212 DMRs) were up-regulated, whereas 38% (6,308 DMRs) were down-regulated in the SP group ( Figure 5A , Figure 5 -source data 1). 6055 genes were annotated as DMR-associated with DMRs falling in either the gene body or 5 kb upstream of the transcription start site (TSS), and of these, 1,315 DMR-associated genes overlapped with period-correlated DE genes ( Figure 5B ). Interestingly, for period-correlated DE genes associated with DMRs, in addition to negative correlations, we also observed positively correlated DMRs, indicating both repression and enhancement of functional (Figure 5C -D) as reported by others (Jones, 2012; Rinaldi et al., 2016; Yin et al., 2017b; Harris et al., 2018) .

The overall clustering pattern of the methylomes resembled that of the transcriptomes, indicating an important role for global DNA methylation in regulating the co-expressed genes ( Figure 5-figure supplement 1) . We examined the modules enriched for period-correlated DE genes and found that several hub genes were regulated by differential DNA methylation. For example, the hub gene of the blue module, Htatip2, which exhibited the same expression pattern of the module eigengene ( Figure 6A-B) , was hypermethylated at the promoter region and repressed in the SP group ( Figure 6C-D) . On the contrary, Parvb and Rftn1, two hub genes from negatively correlated modules, were hypermethylated and repressed in the LP group ( Figure 6A-D) . Except for these negative correlations, some genes with hypermethylation in the gene body or enhancer showed enhanced expression levels ( Figure 6-figure supplement 1) , supporting recent findings that DNA methylation in these regions may activate gene expression (Jones, 2012; Rinaldi et al., 2016; Yin et al., 2017b) . To validate the function of DMR-associated DE genes further, we also performed gene knockdown experiments in two different clonal cell lines. Knockdown of Htatip2 and Dusp18 in LP clone#114 significantly shortened period, whereas knockdown of Rftn1 in SP clone#128 significantly lengthened circadian period ( Figure 6E) , consistent with predictions that deficiency of Htatip2 shortens circadian period possibly by activating the AKT/mTOR signaling pathway (Yin et al., 2017a ; Ramanathan et al., 2018) , and that hypomethylation and upregulated expression of Dusp18 lengthens circadian period, possibly by inhibiting the SAPK/JNK signaling pathway (Wu et al., 2006; Chansard et al., 2007; Yoshitane et al., 2012) .

To assess the role of DNA methylation in circadian periodicity further, we manipulated global DNA methylation either by knocking down DNA methyltransferases or by applying small molecule inhibitors. Interestingly, deficiency of Dnmt1 significantly shortened period length, whereas knockdown of Dnmt3a slightly, but significantly, lengthened period ( Figure 7A -B). Dnmt1 and Dnmt3a knockdown in the ten clonal cell lines showed the same overall results, suggesting that DNA methylation affects circadian periodicity in the same way in all clones tested ( Figure 7C ). As pharmacological validation, administration of SGI-1027, which selectively induces degradation of the DNMT1 protein (Datta et al., 2009) , significantly shortened period, while administration of zebularine, which induces significant reduction of both DNMT1 and DNMT3A (Billam et al., 2010; You and Park, 2012) , lengthened period ( Figure 7D ). Drug administration in primary MEF cells with PER2::LUCsv and NIH3T3 cells carrying an E2-box-luc reporter also revealed similar results ( Figure 7E ). Taken together, these findings suggest that different DNA methyltransferases contribute to the regulation of circadian periodicity, likely via different mechanisms.

Using clonal cell analysis, we show that the heterogeneity of single-cell circadian periodicity is heritable and stable for at least 20 cell passages. The heritability of circadian period is consistent with an epigenetic mechanism, likely mediated by DNA methylation. By analyzing gene expression profiles of multiple clonal cell lines with different circadian periods, we identified groups of differentially expressed genes that were significantly correlated with period length. Although a few core clock genes were differentially expressed in parental cultures, there were no significant differences in these genes among subclones, suggesting they are not responsible for the period heterogeneity seen in these homogeneous cell populations. By comparing subclones, we narrowed down the common candidate gene list and further validated that 86% of the novel candidates regulated circadian period using gene knockdown assays. While some of these genes had effects on period length that were aligned with our predictions, others had effects counter to our expectations which were probably masked in the complex gene networks. Overall, our results are consistent with the hypothesis that period is determined by the ensemble interactions of many genes that can either shorten or lengthen period individually. Importantly, the vast majority of the DE genes identified here have never been reported as having effects on circadian period. Thus, we have provided a new pool of candidate genes functioning in circadian periodicity. We also provide evidence that the genome-wide DNA methylation landscape underlies much of the complex gene networks. Multiple hub genes of period-correlated modules were under the regulation of DNA methylation, showing remarkable coherence in DNA methylation, gene expression, and circadian phenotype. The similar clustering patterns of transcriptomes and methylomes further suggested an important role of DNA methylation in shaping circadian period heterogeneity through regulating large-scale gene networks. Previous studies have linked DNA methylation of core clock genes with different diseases (Joska et al., 2014; Peng et al., 2019) ; however, the results presented here have revealed how global DNA methylation can regulate circadian clock function via genomewide changes in gene expression. Our whole exome sequencing failed to detect significant coding mutations, further supporting the role of differential DNA methylation in establishing circadian heterogeneity. However, we cannot rule out that genetic variation in regulatory regions, or other epigenetic modifications could be involved. Additional experiments will help to understand better the full array of underlying mechanisms regulating circadian period. We observed both negatively and positively correlated DMRs in almost equal proportions, indicating both repression and activation of gene expression by DNA methylation and supporting the revised view of the functions of DNA methylation (Greenberg and Bourc'his, 2019) . In addition, we found that knockdown of DNMT1 and DNMT3A had opposite effects on circadian period. It is not surprising that DNMT1 knockdown alters period length, since it is the methyltransferase responsible for DNA methylation maintenance through mitotic inheritance (Jones, 2012) . However, as DNMT3A is responsible for de novo DNA methylation, it is less clear how its knockdown affects circadian period. One possibility is that DNMT3A is also involved in transcriptional activation associated with active enhancers (Rinaldi et al., 2016; Lyko, 2018) . Another possibility is that some genes might undergo dynamic demethylation and de novo methylation since both Tet2 and Tet3 are expressed at comparable levels to Dnmt3a in our cellular system (Oh et al., 2018; Oh et al., 2019 ; Figure 2 source data 1). Additional studies targeting DNMT1 and DNMT3A may help to explain the functions of different DNMTs in circadian regulation.

In conclusion, our findings have identified a novel pool of candidate genes involved in circadian period regulation, and have revealed the important role of DNA methylation underlying circadian period heterogeneity by bidirectionally regulating large-scale gene co-expression networks. Our study not only expands the knowledge about circadian clock regulation, but also may benefit epigenetic research by providing multiple candidate genes repressed or activated by DNA methylation. embryos. NIH3T3 cells stably expressing Per2 E-box (E2)-driven luciferase bioluminescence reporter were established by lentivirus transduction followed by blasticidin selection. Our cell line stocks have all tested negative for mycoplasma contamination. For authentication of cell lines, as described below, two clonal cell lines, #33 and #114, were sequenced by whole exome sequencing; and a total of 34 clones and subclones were assessed by RNA-seq and were found to be valid.

To measure luminescence rhythms from 35 mm culture dishes, confluent cells were synchronized with 100 nM dexamethasone for 2 hr, then changed to HEPES-buffered recording medium containing 2% FBS (Welsh et al., 2004) , and loaded into a LumiCycle luminometer (Actimetrics). The period was analyzed with LumiCycle Analysis program (Actimetrics). All LumiCycle period analysis results shown in this paper were averages of !3 experiments. Baseline-subtracted signals were exported to Excel to generate bioluminescence traces. For single-cell imaging, cells were changed to recording medium containing 2% B27% and 1% FBS without dexamethasone synchronization. An inverted microscope (Leica DM IRB) in a heated lucite chamber custom-engineered to fit around the microscope stage (Solent Scientific, UK) kept the cells at a constant 36˚C was mounted on an anti-vibration table (TMC) equipped with a 10X objective. A cooled CCD camera with backside illuminated E2V CCD 42-40, 2048 Â 2048 pixel, F-mount adapter, À100˚C cooling (Series 600, Spectral Instruments) was used to capture the luminescence signal at 30 min intervals, with 29.6 min exposure duration, for at least 12 days. 8 Â 8 binning was used to increase the signal-to-noise ratio. The bioluminescence signal of each single cell, outlined with a region of interest (ROI), was tracked using ImageJ (Schindelin et al., 2012; Rueden et al., 2017) with the Trackmate plugin (Tinevez et al., 2017) and analyzed as described previously (Li et al., 2020) .

For exome sequencing, two clonal cell lines #33 and #114 were sequenced representing short period and long period clones, respectively. Genomic DNA was purified using a ChargeSwitch gDNA Mini Tissue Kit (Invitrogen). Libraries were made using the SureSelect XT Reagent Kit (Agilent) following the manufacturer's instruction. All reads were mapped to mm10 genome assembly. We used Haplo-typeCaller and UnifiedGenotyper from GATK to call variants and the results were the union of both callers. SnpEff was used to annotate variants. Results were further filtered as follows: threshold GQ ! 20, total counts ! 8, and alternate frequency (defined as the ratio of alternates to total counts)!30%.

For RNA-seq, cells were collected at two time-points after synchronization: the first peak (T1) and the following trough (T2) based on LumiCycle recording. At each time-point, we collected 2 replicates for 10 clonal cell lines and 1 replicate for 24 subclones. RNA was isolated using TRIzol (Life technologies), and libraries were prepared as described previously (Takahashi et al., 2015) . Raw reads were tested for quality using FastQC. The resulting reads were mapped to mm10 annotation from UCSC using TopHat (Trapnell et al., 2009) . The output BAM file was then filtered for uniquely mapped reads using Samtools (Li et al., 2009) , and RPKM calculations were performed using ana-lyzeRepeats.pl of HOMER suite (Heinz et al., 2010) .

The average RPKM value for each gene was calculated separately for each of the six groups (SP, LP, SSP, LSP, SLP, LLP). To identify significant DE genes, the list was further filtered based on expression level. Only genes for which the maximum average RPKM value among six groups was greater than 0.5 were preserved. Differential gene expression analysis was carried out with DESeq2 (Love et al., 2014) and edgeR (Robinson et al., 2010) using a raw read counts matrix generated with featureCounts tool (Liao et al., 2014) . Genes with FDR < 0.05 were deemed significant. Results from both programs were combined to generate a final DE gene list. Pearson correlation coefficient between circadian period length and gene expression was calculated across all 88 samples (including replicates and different time-points) in Excel. P-value was adjusted using Benjamini-Hochberg (BH) method, and FDR < 0.05 was considered as significant. The overlaps between significant DE genes and period-correlated genes were defined as period-correlated DE genes. Multidimensional scaling (MDS) analysis with Euclidean distance was performed using edgeR. Ingenuity Pathway Analysis (Qiagen) was used to identify the pathways associated with period-correlated DE genes, using all expressed 22,786 genes (average RPKM >0) as a reference set.

For DNA methylation sequencing, cells were collected at the first peak (T1) after synchronization. Each clone included two replicates. DNA was purified using a PureLink Genomic DNA Mini Kit (Invitrogen). Libraries were made using the Premium Reduced Representation Bisulfite Sequencing (RRBS) Kit (Diagenode) following the manufacturer's instruction. Raw reads were tested for quality using FastQC and trimmed with Trim Galore. The trimmed reads were aligned to mm10 using Bismark (Krueger and Andrews, 2011) . The CpG reports from Bismark methylation extractor were then analyzed using methylKit (Akalin et al., 2012) . We used default settings to discard bases that had coverage below 10X and/or more than 99.9th percentile of coverage in each sample. Differentially methylated regions (DMR) were identified using a tiling window of 1,000 bp and a step size of 1,000 bp comparing SP and LP group. Clone#44 was excluded for DMR analysis because of the outlying clustering ( Figure 5-figure supplement 1 ). Overdispersion correction with Fisher's extract test was applied. P-value was adjusted with BH method. DMRs with FDR < 0.05 and methylation difference >25% were considered as significant. Genes with significant DMRs located either in the gene body or 5 kb upstream of the transcription start site (TSS) were considered as DMR-associated genes. Principal component analysis (PCA) was performed using methylKit. All sequencing was performed by the UTSW McDermott Sequencing Core Facility.

Weighted gene co-expression network analysis was performed using WGCNA package (Langfelder and Horvath, 2008) . Only genes for which the maximum average RPKM value among six groups was greater than 0.5 RPKM were used. A soft-threshold power was automatically calculated to achieve approximate scale-free topology (R 2 >0.85). Networks were constructed with blockwiseModules function with biweight midcorrelation (bicor). We used corType = bicor, networkType = signed, TOMtype = signed, TOMDenom = mean, maxBlockSize = 16000, mergingThresh = 0.10, minCoreKME = 0.4, minKMEtoStay = 0.5, reassignThreshold = 1e-10, deepSplit = 4, detectCutHeight = 0.999, minModuleSize = 100, power = 26. The modules were then determined using the dynamic tree-cutting algorithm. Deep split of 4 was used to split more aggressively the data and create more specific modules. Spearman's rank correlation was used to compute module eigengene -covariates associations. Gene set enrichment applied for module -period-correlated DE genes was performed using a Fisher's exact test in R with the following parameters: alternative = 'greater', conf.level = 0.99. The PPI network was generated using STRING without textmining, and the minimum required interaction score was 0.7 (Szklarczyk et al., 2019) .

Gene Knockdown Assay shRNA sequences were cloned into pLKO.1-TRC vector (gift from David Root, Addgene plasmid # 10878) (Moffat et al., 2006) . Scramble shRNA (5' -CCTAAGGTTAAG TCGCCCTCG-3') was used as control. Lentiviruses were produced using HEK293T cells as described previously (Huang et al., 2012) . Viruses were harvested twice after transfection, at 48 and 72 hr, to infect fibroblasts. Forty-eight hours after first infection, cells were synchronized and loaded for Lumi-Cycle analysis. RNA was extracted at the first peak after synchronization to check knockdown efficiency via qPCR. Average of three reference genes (Gapdh, Hprt and Ywhaz) served as internal control. See Supplementary files 2 and 3 for shRNA target sequences and primer sequences, respectively.

The EIF2 signaling pathway activator halofuginone (Sigma-Aldrich) was dissolved in DMSO as 10 mM stock and used at 50 nM. Tunicamycin (Sigma-Aldrich) was dissolved in DMSO as 5 mg/ml stock and used at 5 mg/ml. Cells were treated for 4 hr and 6 hr, respectively, before loading for LumiCycle analysis. DNMT inhibitor SGI-1027 (Sigma-Aldrich) was dissolved in DMSO as 200 mM stock and used at 10 mM. Zebularine (Sigma-Aldrich) was dissolved in water as 200 mM stock and used at 50 mM or 100 mM. The parent culture was continuously treated for up to 60 days and split when necessary. MEFs and NIH3T3 cells were treated for 3 days.

methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles

Effects of a novel DNA methyltransferase inhibitor zebularine on human breast Cancer cells

The causes and consequences of genetic heterogeneity in Cancer evolution

Advances in epigenetics link genetics to the environment and disease

c-Jun N-terminal kinase inhibitor SP600125 modulates the period of mammalian circadian rhythms

Identification of diverse modulators of central and peripheral circadian clocks by high-throughput chemical screening

A new class of quinoline-based DNA hypomethylating agents reactivates tumor suppressor genes by blocking DNA methyltransferase 1 activity and inducing its degradation

The diverse roles of DNA methylation in mammalian development and disease

A DNA methylation reader complex that enhances gene transcription

Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities

Exploring expression data: identification and analysis of coexpressed genes

Crystal structure of the heterodimeric CLOCK:BMAL1 transcriptional activator complex

BioVenn -a web application for the comparison and visualization of biological lists using area-proportional venn diagrams

Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals

Functions of DNA methylation: islands, start sites, gene bodies and beyond

Regulated DNA methylation and the circadian clock: implications in Cancer

Single-cell epigenomics: recording the past and predicting the future

Emergence of noise-induced oscillations in the central circadian pacemaker

The circadian clock and pathology of the ageing brain

Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications

WGCNA: an R package for weighted correlation network analysis

Persistent cell-autonomous circadian oscillations in fibroblasts revealed by six-week single-cell imaging of PER2::LUC bioluminescence

The sequence alignment/Map format and SAMtools

Noise-driven cellular heterogeneity in circadian periodicity

featureCounts: an efficient general purpose program for assigning sequence reads to genomic features

Cellular construction of a circadian clock: period determination in the suprachiasmatic nuclei

Multi-omic measurements of heterogeneity in HeLa cells across laboratories

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

The DNA methyltransferase family: a versatile toolkit for epigenetic regulation

A lentiviral RNAi library for human and mouse genes applied to an arrayed viral high-content screen

Central and peripheral circadian clocks in mammals

Circadian gene expression in individual fibroblasts: cell-autonomous and self-sustained oscillators pass time to daughter cells

Cytosine modifications exhibit circadian oscillations that are involved in epigenetic diversity and aging

Circadian oscillations of cytosine modification in humans contribute to epigenetic variability, aging, and complex disease

The eIF2a kinase GCN2 modulates period and rhythmicity of the circadian clock by translational control of Atf4

DNA methylation of five core circadian genes jointly contributes to glucose metabolism: a Gene-Set analysis in monozygotic twins

Nature, nurture, or chance: stochastic gene expression and its consequences

mTOR signaling regulates central and peripheral circadian clock function

Noise in gene expression: origins, consequences, and control

Dnmt3a and Dnmt3b associate with enhancers to regulate human epidermal stem cell homeostasis

edgeR: a bioconductor package for differential expression analysis of digital gene expression data

ImageJ2: imagej for the next generation of scientific image data

Fiji: an open-source platform for biological-image analysis

A central role for ubiquitination within a circadian clock protein modification code

STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets

The genetics of mammalian circadian order and disorder: implications for physiology and disease

ChIP-seq and RNA-seq methods to study circadian control of transcription in mammals

Transcriptional architecture of the mammalian circadian clock

TrackMate: An open and extensible platform for single-particle tracking

TopHat: discovering splice junctions with RNA-Seq

Qqman: an R package for visualizing GWAS results using Q-Q and Manhattan plots

Dynamic DNA methylation across diverse human cell lines and tissues

Individual neurons dissociated from rat suprachiasmatic nucleus express independently phased circadian firing rhythms

Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression

NRF2 regulates core and stabilizing circadian clock loops, coupling redox and timekeeping in Mus musculus

Ggplot2: Elegant Graphics for Data Analysis

Dplyr: A Grammar of Data Manipulation

Dual specificity phosphotase 18, interacting with SAPK, dephosphorylates SAPK and inhibits SAPK/JNK signal pathway in vivo

TIP30 regulates lipid metabolism in hepatocellular carcinoma by regulating SREBP1 through the akt/mTOR signaling pathway

Impact of cytosine methylation on DNA binding specificities of human transcription factors

Period2 3'-UTR and microRNA-24 regulate circadian rhythms by repressing PERIOD2 protein accumulation

JNK regulates the photic response of the mammalian circadian clock

Zebularine inhibits the growth of HeLa cervical Cancer cells via cell cycle arrest and caspase-dependent apoptosis

This research was supported by the Howard Hughes Medical Institute. All bioinformatics analyses were carried out on Stampede2 cluster of TACC at UT Austin. The authors would like to thank all Takahashi lab members, Dr. Carla B Green, and Dr. Shin Yamazaki for helpful discussions, and the McDermott Bioinformatics Lab at UT Southwestern Medical Center for their bioinformatics support. JST is an Investigator in the Howard Hughes Medical Institute.

Immortalized mouse ear fibroblast cells from male mice carrying PER2::LUCsv bioluminescence reporter were maintained in DMEM (Corning) supplemented with 10% fetal bovine serum (FBS). To generate clonal cell lines, cells were diluted and seeded at a density of~30 cells per 96-well plate with conditioned medium. Each well was monitored on a daily basis to make sure only single colonies were picked. 20 clonal cell lines were randomly selected and cultured continuously for 20 generations (3 days/generation) to verify stability of circadian period. Primary mouse embryonic fibroblast (MEF) cells carrying PER2::LUCsv bioluminescence reporter were isolated from 13.5 day mouse

Statistical analysis of single-cell imaging was performed with a Python code as described previously (Li et al., 2020) . Student's T-test and two tailed F-test were performed in Excel. P-values were adjusted using Benjamini-Hochberg (BH) method. Two-way ANOVA analysis with multiple comparisons via Tukey test was performed using GraphPad Prism. Heatmaps for single-cell imaging analysis and gene expression were generated using MeV based on z-score. GraphPad prism was used to generate heatmaps for T-test and F-test based on log transformed q-value. Volcano plot was generated in R using ggplot2 (Wickham, 2016) . Venn diagrams were generated using BioVenn (Hulsen et al., 2008) . Manhattan plots were generated in R using qqman (Turner, 2014) . Quadrant plots were generated using dplyr package in R.

Funding

Howard Hughes Medical Institute

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.