Dawn- and dusk-phased circadian transcription rhythms coordinate anabolic and catabolic functions in Neurospora

A total of 345 rhythmically transcribed genes (RNAPII-S2P rhythm, P <0.05) did not show corresponding in phase transcript abundance rhythms (Group 2) (see Additional file 5: Table S4). Heat-map analysis suggested that the transcript levels of these genes increase with a substantial delay after the transcription rhythm (RNAPII-S2P profile). Hence, these RNAs could have a long half-life resulting in a delayed phase and a blunted abundance rhythm [7]. Moreover, RNAPII-S2P amplitude of the oscillations of these genes was lower compared to genes identified by both methods (Figure 2B).

Finally, 170 genes with significant RNA abundance rhythms did not exhibit rhythmic RNAPII-S2P occupancy profiles (Group 3) (see Additional file 6: Table S5). This class of genes is potentially interesting since the RNA abundance rhythms might be based on hitherto unknown post-transcriptional mechanisms that are controlled by the circadian clock. However, the average RNAPII-S2P read coverage of these genes was low in comparison to genes that have significant profiles in RNAPII-S2P occupancy and RNA abundance (Figure 2C) (P <10⁻⁴). Thus, we cannot rigorously exclude transcription based expression rhythms of these genes. Together, the data suggest that expression of ccgs in Neurospora is predominantly associated with rhythmic transcription. Transcription independent generation of circadian transcript abundance rhythms, for example, on the level of rhythmic RNA turnover, may thus not be a major pathway used by the circadian clock.

A very recent RNA-seq analysis of three replicate time-courses by Hurley et al. [31] identified 872 genes with circadian RNA abundances rhythms. Of these 872 genes, 697 were expressed under our growth conditions and 327 of them were assigned as rhythmic in our study (see Additional file 7: Figure S2A, Additional file 4: Table S3). The remaining 370 genes were expressed at low levels under our conditions and oscillated with low amplitude in our and the Hurley et al. study (see Additional file 7: Figure S2B-D). A comparison of RNA abundance rhythms and RNAPII-S2P profiles of the rhythmic genes identified by both studies suggests that these genes are controlled on the level of transcription (see Additional file 7: Figure S2E).

Hurley et al. found among 187 luciferase reporters two genes (ccg1 and ccg9) with rhythmic RNA but no luciferase rhythm of the corresponding reporter gene (compare Additional file 8: Figure S6 and Additional file 9: Table S10 in Hurley et al.), suggesting a posttranscriptional regulation. We analyzed the RNAPII-S2P profiles of these genes by ChIP-seq and ChIP-PCR and found that both genes are rhythmically transcribed (see Additional file 7: Figure S2F and G).

We then analyzed the expression phases of ccgs that might be indirectly induced by the CSP1 repressor (303 genes) or indirectly repressed by the WCC activator (191 genes) (see Additional file 12: Figure S3). These ccgs were also expressed in either a morning- or evening-specific manner. The 191 ccgs that were upregulated in ∆wc2 were preferentially dusk-phased (Figure 3A and Additional file 12: Figure S3A) (see Additional file 10: Table S6) while the 303 ccgs upregulated in a CSP1 overexpressing strain (csp1_OE) were mainly dawn-phased (Figure 3B and Additional file 12: Figure S3B) (see Additional file 11: Table S7). Together, WCC and CSP1 appear to regulate, directly or indirectly, at least 1,137 of the 1,407 genes (approximately 80%) that were rhythmically transcribed in wt. Indeed, we found that 457 of 1,407 rhythmically expressed genes have CSP1 binding sites in their upstream regions suggesting a direct regulation by CSP1 (Additional file 11: Table S7). Hence, WCC and CSP1 are major determinants of clock-controlled transcription and circadian phase.

We have previously shown that CSP1 inhibits wc1 transcription and thereby regulates WCC expression levels [34]. Therefore, we also analyzed the expression of frq, which is a direct target of WCC. In ∆csp1 the levels of frq RNA were slightly elevated and the expression phase was slightly advanced (Figure 4D), consistent with the shorter period rhythm observed in ∆csp1 in high glucose conditions [34].

Interestingly, we identified 1,867 genes that were rhythmically expressed in ∆csp1 (P <0.05) but not in wt (see Additional file 14: Figure S4B and Additional file 15: Table S9). The expression rhythms of these ∆csp1-specific ccgs also were clustered in two phases. About 2/3 of the ∆csp1-specific ccgs were expressed with a dusk- and 1/3 with a dawn-specific phase. Intriguingly, the dusk- and dawn-phased ∆csp1-specific ccgs also clustered in two groups in wt, displaying rather complex temporal expression profiles (see Additional file 14: Figure S4B). The data suggest that expression of these genes is clock-controlled. However, under the conditions analyzed, that is, 2% glucose-containing medium, their circadian regulation in wt appears to be blunted by the glucose-dependent activity of CSP1. Examples of these ∆csp1-specific genes are shown in Additional file 14: Figure S4C.

Together, the data suggest that CSP1 regulates the phase and amplitude of both dusk-phased and dawn-phased rhythmic genes. The over-expression of CSP1 indicates that dusk-phased genes are repressed by CSP1 in the subjective morning whereas dawn-phased genes are indirectly activated via unknown pathways. Rhythmic expression of only a fraction of dusk-phased genes is affected in a ∆csp1 strain, suggesting that other transcription factors (TFs) contribute to evening-specific gene expression. The increased accumulation of WCC in ∆csp1 [34] could additionally affect expression of dawn- and dusk-phased ccgs, suggesting a further mechanism by which CSP1 may affect circadian gene expression.

Examples of dawn-phased genes involved in cell rescue and defense are shown in Figure 6C. The dawn-phased expression of heat-shock factors may contribute to the adaptation of Neurospora to higher temperature expected during the day. In order to prevent desiccation during the day fungi (as well as bacteria and plants) synthesize trehalose. Neurospora can synthesize trehalose via a one-step pathway by the trehalose synthase encoded by the morning-specific clock-controlled gene-9 (ccg-9) [45] and via a two-step reaction requiring α, α-trehalose phosphate synthase and trehalose phosphatase (see Additional file 16: Figure S5B). The corresponding genes are rhythmically expressed with a morning-specific phase (Figure 6C, Additional file 1: Table S1 and Additional file 3: Table S2). It is also noteworthy that the Neurospora lysozyme gene (lyz) was rhythmically expressed with peak levels at dawn (Figure 6C). Although the function of lysozymes in fungi is not well characterized, a potential anti-microbial activity might preferentially be required during the day, when bacterial growth is supported by elevated temperature.

A recent study revealed a role of the circadian clock in the coordination of ribosome biogenesis in mammals [46]. In Neurospora, the genes rpc-19, rpb-6 and rpa-12, which encode subunits of RNA polymerase I, were rhythmically expressed with a dusk-phased peak (Figure 6D). Moreover, rpf-2 and dbp-8, genes involved in rRNA maturation and ribosome assembly, were rhythmically expressed in an evening-specific manner. The data suggest that the circadian clock of Neurospora might affect ribosome biogenesis and protein translation.

Genes involved in DNA replication and cell division were also enriched among the dusk-phased circadian genes (Figure 6E). For example, mcm-3, mcm4 and mcm-5 encoding subunits of DNA replication licensing factor, required to initiate DNA replication in eukaryotes [47], were rhythmically expressed with a peak around dusk. The other subunits, mcm-2, 6, and 7, did not qualify as rhythmic genes by our criteria but appeared to be expressed with low-amplitude dusk-phased rhythms (see Additional file 8: Figure S6A). In addition rpa-1 and rpa-2, encoding the subunits of hetero-trimeric replication protein A, also showed circadian expression with a similar peak time with mcm genes. (Figure 6E and Additional file 8: Figure S6B). Furthermore, rfc-2, which encodes a subunit of the hetero-pentameric clamp loader complex, was expressed with a dusk-phased rhythm. The genes encoding the remaining subunits of the clamp loader were potentially also expressed in an evening specific manner (see Additional file 8: Figure S6C). Moreover, the genes encoding the core histones hH2A, hH2B, hH3 and hH4 were expressed at higher levels during dusk (Figure 6E, Additional file 8: Figure S6D). Finally, the gene encoding for the reverse transcriptase subunit of telomerase was robustly rhythmic with an evening specific phase (Figure 6E). The data strongly suggest that genes required for various aspects of cell division and growth are expressed in an evening-specific manner under the control of the circadian clock of Neurospora.

Together, the bi-phasic clustering of ccgs and the pronounced separation of the corresponding gene functions suggest a global and comprehensive temporal coordination of gene expression by the circadian clock to support rhythmic growth of Neurospora.

Our data suggest that clusters of functionally related genes are expressed in a morning- or evening-specific manner under the control of the Neurospora clock.

Neurospora strains wt (FGSC#2489), ∆csp1 (FGSC#11348) and ∆wc2 (FGSC#11124) were acquired from Fungal Genetics Stock Center (FGSC, Manhattan, KS, USA). Vogel's medium (1 X) supplemented with 2% glucose, 0.5% L-arginine and 10 ng/ml biotin was used as a standard growth medium. For the circadian time course experiments, 200 ml medium in 500 ml flasks were inoculated with approximately 10⁶ conidia and grown in light for six hours for conidial germination. Afterwards cultures were entrained in 11 hour/11 hour light/dark cycles for 2 days and released to constant dark for 22 hours. Cultures were transferred in a staggered manner, so harvesting of the cultures was performed within 10 hours.

Cultures were grown in Vogel’s medium with 2% Avicel as the carbon source for the cellulase assay. During harvesting, cultures were washed with three culture volumes of cold water to minimize extracellular cellulase contamination. Due to the adherence of Neurospora to Avicel we can measure extracellular cellulase that is trapped in the Avicel/Neurospora mesh in addition to intracellular cellulase in the secretory pathway. Extraction buffer (50 mM Tris–HCl pH 7.5, 150 mM NaCl, 1.5 mM MgCl₂, 0.1% NP-40, 1 mM phenylmethylsulfonyl fluoride (PMSF), 1 μg/ml leupeptin, 1 μg/ml pepstatin A) was mixed with 500 μl ground mycelia and incubated on ice for approximately 30 minutes with frequent vortexing. The cell homogenate was centrifuged at 20,000 g at 4°C for 30 minutes. Protein concentration was determined by NanoDrop 1000 spectrometer. Cellulase activity was calculated by using Azo-CM-Cellulose (Megazyme, Bray, Ireland) according to the manufacturer's protocol with minor modifications. Briefly, 120 μg total protein extract diluted in 200 μl 100 mM NaOAc, pH 4.6 was mixed with 200 μl Azo-CM-Cellulose and incubated at 40°C for one hour with constant shaking. A 1 ml precipitation solution (300 mM NaOAc, 20 mM ZnAc, pH 5, 75% EtOH) was added to the protein-Azo-CM-Cellulose mixture and centrifuged at 10,000 g at room temperature for 10 minutes. The color of the supernatant was measured at 590 nm with a Jenway 6320D spectrophotometer. To generate the standard curve for cellulase activity (see Additional file 17: Figure S7C), cellulase from Aspergillus niger (1U/μg) (Sigma-Aldrich Chemie Gmbh Munich, Germany) was used.

RNA was extracted with peqGOLD TriFAST (peqLab, Erlangen, Germany) according to the manufacturer’s protocol. RNA was dissolved in 70 μl nuclease free water with 80u Ribolock RNAse inhibitor (ThermoScientific, Waltham, MA US). For the cDNA preparation Maxima First Strand cDNA Synthesis Kit (ThermoScientific, Waltham, MA US) was used. Transcript levels were analyzed by quantitative real-time PCR in 96-well plates with the StepOnePlus Real-Time PCR System (Life Technologies, NY, USA) using TaqMan Gene Expression Master Mix (Life Technologies, NY, USA). Primers and probes are listed in Additional file: Table S11. rRNA was used for normalization. 18

ChIP was performed as described previously [32] by using specific antibody for serine-2 phosphorylated RNAPII C-terminal tail. Polyclonal anti-rabbit RNAPII Ser2-P was raised against the peptide (pS)PTSPSY(pS)PTSPSC. Primers and probes used for ChIP-PCRs are listed in Additional file 18: Table S11. actin gene (ncu04173) was used for normalization.

cDNA was prepared by using NEBNext® Ultra RNA Prep kit with NEBNext® Multiplex oligos according to the manufacturer’s instructions. ChIP DNA libraries were prepared with NEBNext® ChIP-Seq Library Prep Reagent Set for Illumina® with NEBNext® Multiplex oligos. A 2100 Bioanalyzer was used to check the size and the quality of the libraries. Un-paired sequencing with 50 bp reads was performed with a HiSeq 2000 at GeneCore EMBL Heidelberg for RNA-seq and by the BGI, Hong Kong, for ChIP-seq. Individual sequence reads for each run are available in the Sequence Read Archive (SRA) database under the study name PRJNA248256. Accession numbers for experiments and number of sequence reads are listed in Additional file: Table S12. 19

Raw sequence reads were mapped to the Neurospora crassa genome (NC10) using Bowtie [61], where parameters were set to allow maximum three mismatches and suppress alignments which mapped to more than one location. Gene expression was quantified by the number of reads falling into the annotated exons. For analysis of the RNAPII-S2P ChIP-seq the reads that fall to the 500 bp window upstream of gene end position were counted. Normalization was carried out using the size factor formula as described [62].

In order to prevent a read contamination from neighboring genes, an ‘Influence factor’ was calculated as shown below:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Influence\kern0.5em =\kern0.5em \frac{G/L{}_G}{\left(G\kern0.5em -\kern0.5em O\right)/\left(L{}_G\kern0.5em -\kern0.5em L{}_O\right)} $$\end{document}Influence=G/LGG−O/LG−LOwhere G is the read count mapped to the gene, L_G is the length of the gene, O is the read count mapped to the overlapping region and L_O is the length of the overlap. The gene with an Influence factor larger than 2 was filtered from the RNA-seq and ChIP-seq analysis. The genes with low RNA levels (lower 20% of all annotated genes) were excluded from further analysis. In order to assign the cut-off for RNAPII-S2P ChIP-seq, a background read was calculated by counting the reads that fall into 500 bp windows from non-coding regions of the genome. Genes with a lower RNAPII-S2P signal compared to the median of the background reads were excluded from further analysis. Normalized RNA-seq and RNAPII-S2P reads are shown in Additional file 1: Table S1 and Additional file 13: Table S8 (∆csp1 RNA-seq).

The read counts mapped to the genes were assumed to follow a negative binomial (NB) distribution:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {G}_i\approx NB\left({\upmu}_i,{\upsigma}_i\right) $$\end{document}Gi≈NBμi,σiwhere μ_i is the mean and σ_i is the variance. σ_i was estimated based on the mean and variance correlation. A ‘locfit’ R package was used to fit the relationship between mean and variance. Adapting from the Robinson and Smyth Exact Test [63], the two-sided P value can be calculated by using the formula:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\kern0.5em =\kern0.5em \frac{\Sigma_{f\left(a,b\right)\le f\left({G}_{treat},\kern0.5em {G}_{control}\right)}f\left(a,b\right)}{\Sigma f\left(a,b\right)} $$\end{document}p=Σfa,b≤fGtreat,Gcontrolfa,bΣfa,bwhere \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ a+b\kern0.5em =\kern0.5em {G}_i^{treat}+{G}_i^{treat},a,b\kern0.5em \in \kern0.5em 0\dots \left({G}_i^{treat}+{G}_i^{treat}\right),{G}_i^{treat} $$\end{document}a+b=Gitreat+Gitreat,a,b∈0…Gitreat+Gitreat,Gitreat is the reads count of i_th gene in treatment condition, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {G}_i^{treat} $$\end{document}Gitreat is the reads count of i_th gene in control condition, f (a, b) is the product of f (a) and f (b), which can be computed using ‘dnbinom’ of the R package.

The ARSER [36] program was used to identify rhythmic RNA and RNAPII-S2P profiles. Period length was set between 18 to 26 hours. In order to detect rhythmic expression of genes that have very high expression levels in light (0 time point), we ran ARSER twice with and without the 0 time point. All genes that have ARSER detected RNA and/or RNAPII-S2P profiles were selected for further analysis. The ARSER output is shown in Additional file 3: Table S2. Independent phase and amplitude determination was performed by fitting the data to a sine wave using a negative binomial generalized linear model. The ‘glm.nb/function from MASS package of R was used and the period length was set as 22 hours.

In total, 500 simulation data were generated for RNA and RNAPII-S2P by randomly shuffling the 12 time points by using the ‘sample’ function of R. The shuffled data were then analyzed by ARSER twice with and without 0 time point. It was observed that the occurrence of rhythmic genes follows an exponential decay function against the amplitude; exponential decay functions were fitted to both the original and the simulated data. The FDR was then computed by using the following formula: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ FDR(Amplitude)\kern0.5em =\kern0.5em \frac{f_{bg(Amplitude)}}{f_{ob(Amplitude)}\kern0.5em +\kern0.5em {f}_{bg(Amplitude)}}\times \kern0.5em Weight(Amplitude) $$\end{document} F D R Weight Amplitude Amplitude = × f b g Amplitude f o b Amplitude f b g Amplitude +

Where, f_{bg(Amplitude)}, f_{ob(Amplitude)} are the estimated occurrence of rhythmic genes against the amplitude using the fitted exponential decay function, and Weight (Amplitude) is introduced to remove the bias for the smaller number of occurrence when the amplitude is high. FDR according to Coverage is computed by using the same formula.