Epigenetic inheritance of circadian period in clonal cells
Abstract
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.
Introduction
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.
Results
Heritable circadian periodicity in clonal cell lines
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–B; 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).
Transcriptomics identifies novel gene candidates determining period length
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, Figure 2—source data 1). 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 demonstrating the overall consistency of the gene knockdowns on circadian period (Figure 3G, Figure 3—source data 1). 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.
Global DNA methylation contributes to 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 gene expression by DNA methylation (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).
Opposite effects of different DNMTs on circadian period
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.
Discussion
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 genome-wide 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.
Materials and methods
Generation of clonal cell lines and cell culture
Request a detailed protocolImmortalized 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 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.
Bioluminescence imaging and data analysis
Request a detailed protocolTo 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).
Next generation sequencing and data analysis
Request a detailed protocolFor 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 SureSelectXT Reagent Kit (Agilent) following the manufacturer’s instruction. All reads were mapped to mm10 genome assembly. We used HaplotypeCaller 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 analyzeRepeats.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 (WGCNA)
Request a detailed protocolWeighted 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 (R2 >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’ -CCTAAGGTTAAGTCGCCCTCG- 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 LumiCycle 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.
Drug treatment
Request a detailed protocolThe 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 µg/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 µM. Zebularine (Sigma-Aldrich) was dissolved in water as 200 mM stock and used at 50 µM or 100 µM. The parent culture was continuously treated for up to 60 days and split when necessary. MEFs and NIH3T3 cells were treated for 3 days.
Quantification and statistical analysis
Request a detailed protocolStatistical 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.
Data availability
RNA Sequencing data have been deposited in GEO under accession codes: GSE132663 and GSE132665. Exome sequencing data have been deposited in SRA under accession number: PRJNA548837. All data generated or analyzed during this study are included in the manuscript and supporting files. Source data have been provided for Figures 2 and 4.
-
NCBI Gene Expression OmnibusID GSE132663. Transcriptional Profiling of Clonal Cell Lines with Different Circadian Period.
-
NCBI Gene Expression OmnibusID GSE132665. RRBS Profiling of Clonal Cell Lines with Different Circadian Period.
-
NCBI BioProjectID PRJNA548837. Exome-seq of mouse immortalized ear fibroblast clonal cell lines with different circadian periods.
References
-
Effects of a novel DNA methyltransferase inhibitor zebularine on human breast Cancer cellsBreast Cancer Research and Treatment 120:581–592.https://doi.org/10.1007/s10549-009-0420-3
-
The diverse roles of DNA methylation in mammalian development and diseaseNature Reviews Molecular Cell Biology 20:590–607.https://doi.org/10.1038/s41580-019-0159-6
-
Functions of DNA methylation: islands, start sites, gene bodies and beyondNature Reviews Genetics 13:484–492.https://doi.org/10.1038/nrg3230
-
The circadian clock and pathology of the ageing brainNature Reviews Neuroscience 13:325–335.https://doi.org/10.1038/nrn3208
-
The sequence alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
Multi-omic measurements of heterogeneity in HeLa cells across laboratoriesNature Biotechnology 37:314–322.https://doi.org/10.1038/s41587-019-0037-y
-
The DNA methyltransferase family: a versatile toolkit for epigenetic regulationNature Reviews Genetics 19:81–92.https://doi.org/10.1038/nrg.2017.80
-
Central and peripheral circadian clocks in mammalsAnnual Review of Neuroscience 35:445–462.https://doi.org/10.1146/annurev-neuro-060909-153128
-
Fiji: an open-source platform for biological-image analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019
-
A central role for ubiquitination within a circadian clock protein modification codeFrontiers in Molecular Neuroscience 7:69.https://doi.org/10.3389/fnmol.2014.00069
-
The genetics of mammalian circadian order and disorder: implications for physiology and diseaseNature Reviews Genetics 9:764–775.https://doi.org/10.1038/nrg2430
-
ChIP-seq and RNA-seq methods to study circadian control of transcription in mammalsMethods in Enzymology 551:285–321.https://doi.org/10.1016/bs.mie.2014.10.059
-
Transcriptional architecture of the mammalian circadian clockNature Reviews Genetics 18:164–179.https://doi.org/10.1038/nrg.2016.150
-
TopHat: discovering splice junctions with RNA-SeqBioinformatics 25:1105–1111.https://doi.org/10.1093/bioinformatics/btp120
-
Dynamic DNA methylation across diverse human cell lines and tissuesGenome Research 23:555–567.https://doi.org/10.1101/gr.147942.112
-
BookGgplot2: Elegant Graphics for Data AnalysisNew York: Springer-Verlag.https://doi.org/10.1007/978-0-387-98141-3
Article and author information
Author details
Funding
Howard Hughes Medical Institute
- Joseph S Takahashi
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
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.
Copyright
© 2020, Li et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,992
- views
-
- 496
- downloads
-
- 12
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Neuroscience
Specialized chemosensory signals elicit innate social behaviors in individuals of several vertebrate species, a process that is mediated via the accessory olfactory system (AOS). The AOS comprising the peripheral sensory vomeronasal organ has evolved elaborate molecular and cellular mechanisms to detect chemo signals. To gain insight into the cell types, developmental gene expression patterns, and functional differences amongst neurons, we performed single-cell transcriptomics of the mouse vomeronasal sensory epithelium. Our analysis reveals diverse cell types with gene expression patterns specific to each, which we made available as a searchable web resource accessed from https://www.scvnoexplorer.com. Pseudo-time developmental analysis indicates that neurons originating from common progenitors diverge in their gene expression during maturation with transient and persistent transcription factor expression at critical branch points. Comparative analysis across two of the major neuronal subtypes that express divergent GPCR families and the G-protein subunits Gnai2 or Gnao1, reveals significantly higher expression of endoplasmic reticulum (ER) associated genes within Gnao1 neurons. In addition, differences in ER content and prevalence of cubic membrane ER ultrastructure revealed by electron microscopy, indicate fundamental differences in ER function.
-
- Medicine
- Neuroscience
Pain after surgery causes significant suffering. Opioid analgesics cause severe side effects and accidental death. Therefore, there is an urgent need to develop non-opioid therapies for managing post-surgical pain. Local application of Clarix Flo (FLO), a human amniotic membrane (AM) product, attenuated established post-surgical pain hypersensitivity without exhibiting known side effects of opioid use in mice. This effect was achieved through direct inhibition of nociceptive dorsal root ganglion (DRG) neurons via CD44-dependent pathways. We further purified the major matrix component, the heavy chain-hyaluronic acid/pentraxin 3 (HC-HA/PTX3) from human AM that has greater purity and water solubility than FLO. HC-HA/PTX3 replicated FLO-induced neuronal and pain inhibition. Mechanistically, HC-HA/PTX3-induced cytoskeleton rearrangements to inhibit sodium current and high-voltage activated calcium current on nociceptive DRG neurons, suggesting it is a key bioactive component mediating pain relief. Collectively, our findings highlight the potential of naturally derived biologics from human birth tissues as an effective non-opioid treatment for post-surgical pain. Moreover, we unravel the underlying neuronal mechanisms of pain inhibition induced by FLO and HC-HA/PTX3.