Circadian deep sequencing reveals stress-response genes that adopt robust rhythmic expression during aging

To identify age-dependent changes in diurnal gene expression, we performed RNA-seq with RNA from heads of 5- or 55-day-old() females collected every 4 h for one daily cycle in two biological replicates in light-dark (LD) 12:12 h. Gene expression, quantified in units of Fragments Per Kilobase of transcript per Million mapped reads (FPKM) using the Tuxedo suite(), showed high correlations between biological replicates in both young and old flies (;). Using ARSER (ref.), we identified 2,036 genes to be rhythmic in young flies, including 67% of genes previously deemed rhythmic at the messenger RNA (mRNA) level in heads of male flies; this overlap is substantial considering differences in sex and periodicity detection methods (Methods). ARSER classified 1,887 genes to be rhythmic in old, including 922 that were also rhythmic in young (). w w 1118 ^{18 19} Supplementary Table 1 Supplementary Fig. 1 Supplementary Table 2 ^{20 3} Supplementary Data 1

Although many CCGs maintained similar rhythmicity patterns throughout aging, a substantial number of genes showed marked changes in phase, amplitude or statistical rhythmic status. For example, 297 genes rhythmic in both young and old flies exhibited a ≥2 h shift in their time of peak expression (). Further, 48 genes were highly rhythmic in young flies but arrhythmic in old, while 38 genes were highly rhythmic in old and arrhythmic in young flies (, Methods). Examples of other age-related expression changes are shown in. RNA-seq expression plots for all FlyBase genes and isoforms are available at:. Supplementary Fig. 2 Fig. 1a Fig. 1b–e http://hendrixlab.cgrb.oregonstate.edu/youngAndOldExpression.html↗

In accordance with the diversity of alterations in CCG expression, expression of some core clock genes also showed significant changes. Althoughexpression remained highly rhythmic throughout aging, it showed a decline in peak expression, in agreement with previous qRT-PCR studies in heads of male and female flies.levels were highly rhythmic and not significantly reduced by day 55 (); however, qRT-PCR with RNA from 75-day-old flies corroborated prior reports of a decline inlevels in very old age(). In contrast with prior qPCR data, RNA-seq revealedpeak expression to increase with age, and this was confirmed by qPCR (). However, PER protein levels decreased significantly by day 55 (;), consistent with previous reports. Thus, the opposite age-dependent changes inandRNA levels suggest that the core circadian mechanism known in young flies may be altered during aging. tim Clk Clk per per tim ^{21 22} Fig. 2a ²¹ Fig. 2b ^{21 22} Fig. 2b Fig. 2c,d Supplementary Fig. 3 ^{21 22}

The most unexpected outcome of our analysis was the identification of a number of genes that adoptedrhythmic expression or robustly increased amplitudes in old flies. We named this gene subset ‘late-life cyclers' or LLCs. To computationally identify LLCs, we developed a scoring metric called the differential rhythmicity score (). Theaccounts equally for both the difference in the rhythmicity score, defined as the negative log of the ARSERvalues for the respective age groups, and the differential robustness, or the log fold change in the max–min expression (see Methods). We identified genes with significantvalues at an FDR of 0.05 (,,). This method successfully identified genes with strong loss or gain in expression rhythmicity with age; however, some LLC-like genes of potential interest that do not make the cutoff are shown in. Another set of genes with rhythms abolished during aging are listed in. A complete list of genes ranked bycan be found in. de novo S S P S S DR DR DR DR Fig. 3a Supplementary Fig. 4a Supplementary Table 3 ²³ Supplementary Fig. 4b Supplementary Fig. 4c Supplementary Data 2

summarizes properties of the top 25 LLCs. 21 have predicted orthologs in humans according to the DRSC Integrative Ortholog Prediction Tool (DIOPT (ref.)). Functional analysis using DAVID v6.7 (ref.) revealed the terms ‘heat shock' and ‘stress response' to be enriched among the top LLCs (FDR 0.05). Published microarray data sets reported 16 of the 25 LLCs to be upregulated under oxidative stress (OS) when tested at one unspecified time of day().shows superimposed RNA-seq expression plots for LLCs that were upregulated by OS in at least two previous single-time-point studies. Strikingly, all of these genes peaked within two hours of night onset, suggesting that they may be governed by common regulatory mechanisms. The strong enrichment of LLCs peaking in the late day/early night starkly contrasts with the global phases of genes rhythmic in young or old flies, which peaked predominantly in the early day/late night (). Supplementary Table 4 ^{24 25 26 27 28 29} Supplementary Table 4 Figure 3b Fig. 3c

We focused on the five LLCs with the most dramatic gains of rhythmicity for further investigations: small heat shock protein; fibroblast growth factor ortholog; lactate dehydrogenase (); Hsp40-like; and, which bears homology to the mammalian histidine-rich glycoprotein () according to DIOPT. Independent qPCR experiments in flies aged to 5, 35, 55 and 75 days confirmed increased expression of these genes in 55-day old and showed even further increase in heads of very old (75 days) females (). Middle-aged females already showed mild increase in the expression of these genes (). Altogether, these data reveal LLC mRNA levels exhibit strong age-dependence. These effects are not sex dependent as LLCs showed similar expression changes in heads of 5-day versus 55-day males (). Hsp22 bnl ImpL3 CG7130 CG15784 HRG Fig. 3d Supplementary Fig. 5 Fig. 3e

We hypothesized that oxidative stress (OS), which increases during aging, might play a role in the rhythmic activation of LLCs in old flies. To test this, we exposed youngflies to continuous hyperoxia (HO; 100% O). As the maximum lifespan under these conditions was 5–6 days, we collected flies at 4 h intervals on the fourth day after HO onset (Methods,). Remarkably, HO induced robust rhythmic expression of the tested LLCs in these young flies (), similar to those seen in old flies (). Furthermore, HO reproduced in young flies the opposing changes inandexpression observed during aging (,). These results suggest that oxidative stress contributes to rhythmic LLC upregulation in aging flies. ^{14 30} Fig. 4a Fig. 4b Fig. 3d Fig. 2b Supplementary Fig. 6a w per tim 2

Because constant hyperoxia induced rhythmic rather than constitutive LLC transcription, we tested whether CLK is involved in LLC regulation by measuring their expression in youngmutants collected on the 4th day in HO. Remarkably, this exogenous OS failed to significantly upregulate,andinflies relative tocontrols in normoxia (). Even for heat-shock LLCsand, the upregulation in HO-treatedflies was less significant than in HO-treatedflies. These results support an essential role for CLK in promoting rhythmic LLC upregulation during OS. Clk ImpL3 bnl CG15784 Clk w Hsp22 CG7130 Clk w out out out Fig. 4c

Interestingly, althoughexpression inmutants was arrhythmic as expected,expression in these HO-treated mutants remained weakly rhythmic, similar to several LLCs (). This suggests that CLK-independent mechanisms also contribute to rhythmic upregulation ofand some LLCs during OS. tim Clk per per out Supplementary Fig. 6b

We performedtranscript assembly using StringTieand Cuffmergeand found 154 unannotated genes with multi-exonic transcripts (). Twenty two of these were rhythmic in old flies (), and five exhibited LLC-like behaviour (). Among these five, only one (hereafter ‘') is conserved across several insects according to the 27 insect alignment and associated phastCons analysis from UCSC Genome Bioinformatics(). Whilepartially overlapped both a TE and a piRNA cluster annotation in the sense strand, the other four fully overlapped transposable elements (TEs) in antisense and mature Piwi-interacting RNAs (piRNAs) in the sense strand, implicating these four as primary piRNA transcripts(). de novo crescendo crescendo ^{31 19} Supplementary Data 3 Supplementary Table 4 Fig. 5a ³² Fig. 5b ^{33 34} Supplementary Fig. 7

We measuredexpression by qPCR in heads of 5-, 55- and 75-day old females, and observed an exponential increase with age, with the greatest increase occurring between days 55 and 75 (). As it oscillated in the same phase as other tested LLCs, we measuredlevels in HO-exposed flies; however, it was not significantly upregulated (), suggesting thatis not induced directly by HO but may be stimulated in response to other age-associated changes. crescendo crescendo crescendo Fig. 5c Supplementary Fig. 8

Our round-the-clock data afforded a high-confidence measure of age-dependent changes in average expression level for all genes, by treating individual samples from different time points as replicates. Among genes with an FPKM>1 in young or old flies, we found 1,504 genes to be significantly downregulated by the age effect (FDR 0.01) and 1,307 genes to be upregulated (;). Of these, the 583 genes upregulated and the 676 genes downregulated by >50% during aging are shown in. We found the ‘housekeeping gene'(actin) in the strongly upregulated subset (), indicating that it is a poor endogenous control for age-dependent qPCR experiments in. Notably, 33.5% of differentially expressed genes (FDR 0.01, fold-change ≥1.5) were robustly rhythmic in young flies, old flies, or both. Fig. 6a Supplementary Data 4 Fig. 6b Supplementary Fig. 9 Act5C Drosophila

Functional analysis with DAVID revealed several enriched annotation clusters for differentially expressed genes. Gene ontology terms related to immune response, glutathione metabolism and response to cellular and genotoxic stress were enriched among upregulated genes; terms associated with neural function, locomotory behaviour, ion homoeostasis and response to entrainment cues were enriched among downregulated genes (;). Fig. 6 Supplementary Data 5

were raised on a standard yeast (35 g l), cornmeal (50 g l) and molasses (5%) diet at 25±1 °C, under a light-dark (LD) 12:12h regimen. Mated flies were kept in groups of 50 males or 50 females in 300 ml round bottom polypropylene ventilated bottles (Genesee Scientific, San Diego, CA). Diet was changed three times a week without anaesthesia. The following genotypes were used in this study:(control) and(ref.). To induce oxidative stress by hyperoxia (HO), 5-day-old flies of each genotype were placed in clear, airtight chambers with 100% Oflow-through at atmospheric pressure, starting at ZT0. Control flies remained in normoxia (NO) next to the chambers. Two biological replicates of 25–50 flies each in HO or NO were collected simultaneously every 4 h for one 24 h cycle starting at 72 h after HO onset in LD or DD as shown in. Drosophila melanogaster w Clk −1 −1 1118 out ³⁸ Fig. 4a 2

One biological replicate of 50 flies was collected every 4 h for two 24 h cycles starting at ZT0 (lights on) for qRT-PCR. Each sample of 25–50 fly heads was homogenized in TRIzol Reagent (Thermo Fisher, Waltham, MA) using a Kontes handheld motor and pestle, and RNA was extracted according to the manufacturer's instructions. Samples were treated with rDNase I (Takara, Japan), followed by a phenol/chloforom extraction and ethanol/sodium acetate precipitation. Complementary DNA (cDNA) was synthesized from 1 μg of RNA using the iScript cDNA synthesis kit (Bio-Rad, Hercules, CA), or the Maxima First Strand cDNA Synthesis Kit for RT-qPCR (Thermo Fisher, Waltham, MA). Real-time PCR was performed with Power SYBR Green PCR Master Mix (Thermo Fisher, Waltham, MA) on a StepOne Plus Real-Time machine (Applied Biosystems, Foster City, California). Primers were obtained from Integrated DNA Technology (Coralville, Iowa), and all primer sets were verified to have >90% efficiency. Primer sequences are given in. Data were analysed using the 2method, using() as the endogenous control for normalization.was selected based on its low variance between time points and ages according to RNA-seq and qPCR (). Data for ZT 24 of cycle 1 at day 5 is repeated for ZT 0 of cycle 2 at day 5 for the experiment in females aged to 5, 55 and 75 days. Supplementary Table 6 Supplementary Fig. 9 −ΔΔCT Decapping protein 2 DCP2 DCP2

Heads of 5 and 55 day-old females (three biorepeats of 20 heads per timepoint for each age) were homogenized in Laemmli buffer, sonicated, boiled at 100 °C for 5 min and centrifuged at 12,000at 4 °C. A constant ratio of the buffer (7 μl per head) was used to ensure equal protein loading and separation on NuPAGE 4–12% gradient acrylamide gel (Life Technologies). Proteins were transferred to the 0.45 μm polyvinylidene fluoride (PVDF) Immobilon-FL membrane (Millipore Billerica, MA) and stained for 5 min with REVERT Total Protein Stain Kit (Li-Cor Biosciences). After staining, membranes were scanned in the 800 channel on the Odyssey Infrared Imaging system to quantify the total protein. The staining was reversed using the same kit, and the membranes were incubated in 1 × TBST (10 mM Tris, 0.15 M NaCL, 0.1% Tween-20, pH 7.5)+5% dry milk for 2 h, then overnight at 4 °C with 1:15,000 anti-PER (gift from Dr J. Price)in blocking buffer. Membranes were treated for 2 h with 1:20,000 goat anti-rabbit IRDye800 (catalogue # 926–32211, LI-COR Biosciences, Lincoln, NE) diluted in Odyssey Blocking Buffer. After washes, PER signal was quantified relative to total protein using the LI-COR Odyssey Infrared Image Studio software according to the manufacturer's instructions. Uncropped version of the PER western blot is present in. g ³⁹ Supplementary Fig. 3

Following poly(A) selection of RNA samples purified as above, strand-specific cDNA libraries were prepared using Illumina's TruSeq stranded mRNA kit (100 bp paired-end) following the manufacturer's directions, then sequenced on Illumina HiSeq 2000 with parallel samples from young and old heads multiplexed in the same lane.

Raw RNA-seq reads were filtered to exclude reads with mean quality scores <30 and to trim 3′ ends with mean quality <30, using the program skewer. Filtered reads were aligned to thegenome (BDGP release 6.06/dm6) using TopHat version 2.0.14 with a max intron length of 10,000. Aligned reads from individual time points were submitted to StringTie v1.2.0 for novel transcript assembly, with minimum junction coverage of 2; release 6.06 of the FlyBase genome annotation was used to guide the assembly process. The resulting list of novel transcripts was refined with Cuffmerge from the Cufflinks package (v2.2.1). Cuffdiff was used to compute abundance of genes and transcripts in units of Fragments Per Kilobase of transcript per Million mapped reads (FPKM) for novel genes and isoforms, separately from those annotated in the FlyBase genome annotation. We note that expression of someassembled transcripts, including the four LLC-like putative primary piRNA transcripts (), could not be reliably quantified because most of their reads aligned with high per cent identity to multiple genomic loci according to Bowtie 2 (ref.). However, each of these four assembled transcripts aligned only once with 100% identity according to the UCSC Genome Bioinformatics BLAT tool. ⁴⁰ Fig. 5 ⁴¹ Drosophila melanogaster de novo

Because of their short length resulting in unstable FPKM calculations, pre-microRNA hairpins were excluded from all subsequent analyses.

All genome browser tracks showing RNA-seq reads were generated in the Integrative Genomics Viewer (IGV). ⁴²

Rhythmic transcripts with 24-h periodicity were identified using ARSER (ref.), JTK_CYCLEand empirical JTK_CYCLE(;). For all three programs, input data was formatted as a series of two daily cycles. Thevalue used as a significance threshold for empirical JTK_CYCLE was the empiricalvalue. Genes reported as rhythmic have a median expression ≥1 FPKM, a max/min fold-change ≥1.5 and avalue≤0.05 in accordance with published thresholds for rhythmicity detection. Although all three methods showed substantial agreement (), ARSER showed the strongest proportional overlap between genes rhythmic in old versus young flies (). Thus, to be conservative when identifying age-dependent changes in rhythmicity, we used ARSER output for subsequent analyses. For basic comparisons, we define ‘highly rhythmic' asvalue≤0.01, and arrhythmic as>0.5. For comparison of overlap between our set of genes rhythmic in young females and genes reported by others as rhythmic in heads of young males, we evaluated the per cent overlap with genes that mapped to IDs in Flybase release 6.06. ^{20 43 44} Supplementary Table 1 Supplementary Fig. 2a ^{5 45} Supplementary Fig. 10 Supplementary Table 7 ³ P P P P P

We assigned anto each gene having an ARSERvalue<1 and median expression ≥1 FPKM in young or old flies, and nonzero expression in at least one time point in young and old flies. We then assignedvalues to these normally distributedvalues and subsequently computed their false discovery rates (FDRs) using the BH procedure to adjust for multiple hypothesis testing (). Among the resulting set of significantly (FDR 0.05) differentially rhythmic genes, we defined our top LLCs as those rhythmic in old flies (max/min fold-change≥1.5 and ARSERvalue≤0.05) with at least one isoform satisfying the rhythmicity criteria imposed at the gene level in old flies (). S P P S P DR DR Supplementary Fig. 4a and b ²³ Fig. 3a

When identifying genes with the strongest improvements in rhythmicity with age, we found that many genes either showed trivial expression levels and low peak/trough fold change in young but robust amplitudes in old, or showed enhanced precision in the periodicity of their expression with age. A few, including, fell into both of these categories. We sought to define a metric that would incorporate these two patterns of differential rhythmicity, and to this end we developed a differential rhythmicity score calculated as follows: ImpL3

is the sum of two Z-scores divided by; the factor ofin the denominator ensures that it obeys a standard normal distribution. The first term,, is a Z-score computed for the age-dependent change in periodicity,, using thevalues from ARSER in young and old. The second term, is a Z-score for the differential robustness Δ, the log fold change in the effective amplitude, given by S Z P Z R DR P R

whereis the effective amplitude, max FPKM–min FPKM. Each of these Z-scores showed a normal distribution, as did the combined score(). Avalue was computed for eachusing a Gaussian distribution based on the fit to the empirical distribution ofvalues, and the BH procedure was used to compute FDRs. S P S S DR DR DR Supplementary Fig. 4a

The custom Perl scripts that we used to implement this analysis can be found at http://hendrixlab.cgrb.oregonstate.edu/LLCs.html↗

For the bar plots in, phases from ARSER were rounded to the nearest integer and mapped to bins of size (4−2, 4+2) for=0, 1, ..., 5, where each integer 4is a time point sampled for RNA-seq. Fig. 3c n n n n

To identify age-induced differential expression independent of time of day, aligned reads for each individual data point (single fly cohort at single-time-point) were treated as replicates segregated by age into either the day 5 or day 55 group. Cuffdiff was run to compute average expression over samples for a given age, and to identify significantly differentially expressed genes. Using all the samples for a given age in a batch analysis using Cuffdiff properly identifies genes that are consistently differentially expressed due to age, regardless of time of day. The DAVID Functional Annotation Tool () was used to identify enriched functional annotations among significantly upregulated or downregulated genes having a minimum FPKM of 1 in young or old and a fold change ≥1.5. In addition, 525 out of the 582 upregulated genes and 572 of the 676 downregulated genes mapped to DAVID IDs. For annotation we included biological process and molecular function gene ontology (GO) terms, as well as KEGG pathways. For functional annotation clustering we used default classification stringency (Medium). The top 10 annotation clusters are summarized next to the heat map in. Because DAVID does not provide names for the annotation clusters, we assigned names that globally represented the terms present in each cluster. Complete results of GO analyses are presented in. https://david.ncifcrf.gov/summary.jsp↗ Fig. 6a Supplementary Data 5

Data for RNA-seq and processed files have been deposited to NCBI Gene Expression Omnibus (GEO) under the accession number. RNA-seq expression plots for all FlyBase genes and isoforms are available at:. GSE81100 http://hendrixlab.cgrb.oregonstate.edu/youngAndOldExpression.html↗

All other data supporting the findings of this study are included in the manuscript and its supplementary files or are available from the corresponding authors on request.