http://orcid.org/0000-0002-9145-5953
Figures
Abstract
There is considerable insight into pathways and genes associated with heat-stress conditions. Most genes involved in stress response have been identified using mutant screens or gene knockdowns. Yet, there is limited understanding of the temporal dynamics of global gene expression in stressful environments. Here, we studied global gene expression profiles during 12 hours of heat stress in the nematode C. elegans. Using a high-resolution time series of increasing stress exposures, we found a distinct shift in gene expression patterns between 3–4 hours into the stress response, separating an initially highly dynamic phase from a later relatively stagnant phase. This turning point in expression dynamics coincided with a phenotypic turning point, as shown by a strong decrease in movement, survival and, progeny count in the days following the stress. Both detectable at transcriptional and phenotypic level, this study pin-points a relatively small time frame during heat stress at which enough damage is accumulated, making it impossible to recover the next few days.
Citation: Jovic K, Sterken MG, Grilli J, Bevers RPJ, Rodriguez M, Riksen JAG, et al. (2017) Temporal dynamics of gene expression in heat-stressed Caenorhabditis elegans. PLoS ONE 12(12): e0189445. https://doi.org/10.1371/journal.pone.0189445
Editor: Denis Dupuy, INSERM U869, FRANCE
Received: August 24, 2017; Accepted: November 25, 2017; Published: December 11, 2017
Copyright: © 2017 Jovic et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The microarray datasets supporting this article have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5753.
Funding: LBS was funded by ERASysbio-plus ZonMW project GRAPPLE - Iterative modelling of gene regulatory interactions underlying stress, disease and ageing in C. elegans (project 90201066) and The Netherlands Organization for Scientific Research (project no. 823.01.001). KJ, JG and JAGR were funded by the Human Frontiers Science Program Organisation (HFSPO), grant RGP0028-2014.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Heat-stress results from an exposure to potentially harmful temperature conditions beyond the optimum range of an organism. An increase in the intracellular temperature can interfere with protein homeostasis, leading to an accumulation of misfolded proteins and protein aggregates [1]. Over the past few decades, detailed insights have been obtained about the molecular and genetic control of the cellular response to heat-stress. To avoid the detrimental effects of cytotoxic misfolded protein species and protein aggregates, multiple stress response systems have evolved as a first line of defence to maintain proteostasis, of which the highly-conserved heat-shock response (HSR) pathway is prominent [1,2].
The accumulation of misfolded proteins is also a hallmark of aging and age-related diseases, such as Alzheimer's and Parkinson's disease [3–5]. The connection between the processes involved in stress and aging is further substantiated by the fact that several components of the stress response pathways were found to function as regulators of lifespan [6–8]. For example, the evolutionary highly conserved transcription factor HSF-1 is a key component in the initiation of the HSR, as well as a regulator of lifespan [9]. Therefore, understanding how an organism perceives and handles heat-stress is fundamental for understanding the molecular mechanisms that underlie aging [4].
The nematode Caenorhabditis elegans is an established metazoan model for studying the effect of—and response to–heat stress in vivo [4,6,7,9]. One of the most widely studied responses in C. elegans is to acute heat stress, which can be easily applied by exposing the animal to temperatures between 33–37°C [9–11]. It was shown that C. elegans detects and responds to heat stress via transient receptor potential channels and a neuropeptide signaling pathway [12], and is capable of swiftly up-regulating a suite of protective proteins (mainly chaperones) to prevent protein denaturation and misfolding, a process which affects every aspect of the animals’ biology [4]. Short or mild stresses can be tolerated and can even protect individuals against future stresses [11].
The effects of heat stress on C. elegans are often quantified on a phenotypic level by recording complex traits such as survival rate, mobility, and reproduction [11,13,14]. Generally, the inflicted damage accumulates with increasing temperature and exposure time. For example, brood size decreases with moderate increases in temperature beyond the optimum [14,15], whereas a strong decrease in survival rates is only observed after prolonged exposures to heat stress [11,16]. At the level of the transcriptome, a heat shock induces a strong response. Genome wide gene expression analysis in C. elegans shows that a two hour exposure to 35°C affects genes associated with development, reproduction and metabolism [17]. Furthermore, an exposure of 30 minutes to 33°C already induced a massive global gene expression shift highly dependent on HSF-1, affecting genes associated with a wide range of functions such as cuticle structure, development, stress response, and metabolism [18].
Given the range of phenotypic effects under different heat stress regimes, it is to be expected that the transcriptional response during heat stress is highly dynamic. For example, the initial transcriptional response to a sub-lethal heat stress likely does not resemble the transcriptome after a heat stress exposure from which the organism cannot recover. So far, most knowledge of the transcriptional stress response in C. elegans was gained from studies using single time point exposures, providing us with a snap-shot of global gene expression [17,18]. Time-resolved global and gene specific expression studies have been performed in yeast and mammalian cells showing the transient nature of the expression changes at sub-lethal heat stress temperatures [19–22]. Overall, there is limited understanding of the temporal dynamics of global gene expression patterns under acute heat stress that lasts from sub-lethal to lethal exposure durations.
To gain more insight into the underlying gene-expression dynamics of the stress response, we have generated a high-resolution time-series of transcriptomic and phenotypic data of C. elegans exposed to heat stress conditions at 35°C for 0-12h. Transcriptomic analysis revealed a global shift in expression dynamics occurring between 3 and 4 hours into the heat exposure. The shift marks the end of an initially highly dynamic transcriptional response to heat stress that plateaus at longer exposures. On a phenotypic level, longer exposures (> 4h) were associated with much lower chances of recovery in the four days following the stress. Therefore, this phenotypic turning point follows shortly after the transcriptional turning point, and is marked by a strong decrease in movement, survival, and progeny count.
Results
Transcriptional variation during prolonged heat stress
We first assessed the impact of heat stress durations on genome-wide expression levels. Wild type Bristol N2 populations were exposed to heat stress conditions at 35°C for increasing exposure durations between 0.5–12 hours (Fig 1A). To find the main sources of variation during the transcriptional response to heat stress, we used principal component analysis (PCA). The first two principal components (PCs) captured 77% (1st 57%, 2nd 20%) of the total variation (Fig 1B). The first PC sorted the time points in chronological order, showing that variation in gene expression between samples was largely due to the increasing length of heat exposure. Furthermore, the distance between samples along the 1st PC was larger for early time points (0–3 hours) in comparison to longer exposure times (3–12 hours). The shift in sample distribution indicates a turning point in expression dynamics in which a large part of the overall changes in gene expression occurred early on in the stress response. Together, the 1st and 2nd PCs place this shift in expression dynamics at 3–4 hours of heat exposure.
(A) mRNA sampling schedule. Bristol N2 populations were grown at 20°C for 46 hours before the start of the heat-shock at 35°C. Clock-symbols indicate the time of sampling for subsequent transcriptome analysis of the dynamic stress response. Each time point (0, 0.5, 1, 2, 3, 4, 6, 8, and 12 hours) was sampled 3–5 times. (B) Principal component analysis of gene expression data averaged per time point. The first two components retain 77% of the variation in the data set, and placed the exposure duration (as indicated by the clock symbol) in chronological order.
Early transcriptional activation of heat shock proteins
Having identified a turning point in transcriptional dynamics, we further investigated the temporal dynamics of expression changes for a set of previously associated heat stress response genes. For C. elegans, the Gene Ontology database listed 72 genes with a role in the ‘response to heat stress’ (GO:0009408, WormBase version 257). Most of these genes show minor transcriptional changes in the course of the 12-hour heat exposure (S1 Fig). This is not surprising, since many components of the (heat) stress response are constituently expressed [23].
The fastest transcriptional response was found for five heat shock proteins, which are part of the heat shock response pathway activated by HSF-1 upon stress exposure. Within 30 minutes of stress exposure, hsp-16.1, hsp-16.2, hsp-16.41, hsp-16.48, and hsp-70 showed a ~16-fold increase in expression levels (S1 Fig). The expression levels of these genes peak 4 hours into the stress exposure, corresponding with the turning point identified in the PCA. Correlation analysis extracted two more genes from our data set that were not listed in the GO term ‘response to heat stress’, but presented with similar expression patterns to the above heat-shock proteins: F13E9.1 (ortholog of human NISCH) and F44E5.5 (member of the hsp-70 family). F44E5.5 was also detected in previous studies [18,24]. Together, this small set of seven genes presented the strongest, first and immediate reaction of the transcriptional response to stress.
Changes in gene expression reach a plateau
We further investigated the temporal dynamics of global transcriptome changes. About ~6200 (~30%) genes contributed significantly (q-value < 0.01) to the variation explained by the first two principal components. This sub-set was used as input for k-means clustering to extract common patterns in gene expression changes, identifying six distinct stress-response groups (Fig 2A and 2B). Cluster 1, 4, and 6 (representing ~3790 (60%) of the genes) contained genes downregulated during exposure to heat (Fig 2A), and cluster 2, 3, and 5 contained upregulated genes (~2450 (40%) of genes; Fig 2B). The largest changes were found in cluster 1 and 3 with an average 5-fold down- and 32-fold up-regulation, respectively.
Genes with similar patterns in expression (log2 ratio with the mean) were grouped by k-means clustering. Dark coloured bold lines present the average expression of the individual clusters; lighter corresponding colours present the expression of individual genes. Enrichment analysis of gene clusters was performed with DAVID 6.8. (A) Cluster 1, 4, and 6 showed a downward trend in gene expression during heat stress. (B) Cluster 2, 3, and 5 were upregulated in response to heat stress. (C) Enrichment of downregulated gene clusters. (D) Enrichment of upregulated gene clusters 2 and 5. Analysis of cluster 3 did not result in a significant enrichment.
Within these clusters, the initial in- or decrease in transcript levels started rapidly, between 0.5–1 hour after initiation of the heat stress exposure, and reached a plateau before 3 hours into the stress response. One exception is cluster 1, consisting of 420 genes that were down-regulated after 4 hours. Interestingly, this was the only pattern clearly distinguishing the later (>4h) time points. Overall, the transcriptional patterns reveal a global change in expression dynamics around 3–4 hours into the stress response, starting with a highly dynamic phase and ending on a plateau with minimal overall changes. To ensure that the observed change in expression dynamics is not due to an increased number of dead animals at the time of sampling, the survival was scored directly after removing the animals from the heat stress in a separate experiment. We found that >90% of animals initially survived an exposure up to 6 hours at 35°C (S2 Table). This shows that the observed transcriptional attenuation at 4 hours occurred before the animals are dead.
To explore the biological functions associated with the gene sets within the individual response clusters, a GO-enrichment analysis was performed (Fig 2C and 2D; S1 Table). Overall, the down-regulated clusters were enriched with structural constituents of the cuticle, particularly collagens (col, dpy, rol, sqt), as well as genes associated with transcription (nhr), metabolic processes, and locomotion (Fig 2C). In the upregulated clusters, genes involved in nucleosome assembly (his) were found to be overrepresented, as well as those regulating embryo and larval development (Fig 2D). Cluster 3, the smallest group (54 genes), had an immediate and strong reaction to the stress. This cluster could not be associated with an enrichment term. Half of the genes within this cluster have not previously been classified with any GO term.
Gene expression dynamics correlated with phenotypic changes
Through transcriptome analysis, we identified a turning point around 3–4 hours into the stress response, separating an initially highly dynamic phase from a later mostly stagnant phase. Next, we tested how these observed transcriptional patterns correlate with the effects of increasing heat stress durations on the phenotypic recovery of the animals. To measure the effects, we observed survival, progeny count, and movement in populations that were allowed to recover at 20°C after different heat stress durations (Fig 3). Since it has previously been shown that it can take three days after the exposure to a transient, severe heat-shock to observe the fatal effects in the survival scores of C. elegans [11], we recorded daily phenotypic observations over a four day recovery period following the stress.
After exposing N2 populations in the L4 stage to increasing heat-shock durations at 35°C, populations were maintained at 20°C. A total of ~65 individuals divided over 3 biological replicas per treatment group were observed. The following phenotypes were scored on the four consecutive days following the heat stress: (A) fraction alive, (B) fraction of worms with a healthy movement phenotype (i.e. sinusoidal, constant and unprovoked movement), (C) the average number of viable offspring produced per population with a cut-off point set at 1000 offspring. A detailed overview of the tested differences between groups can be found in S3 Table.
The heat exposure durations resulted in three phenotypically distinct groups (for details, see S3 Table). First, for survival, the animals exposed to heat for up to two hours show high survival chances equal to the control (Fig 3A). An intermediary group was formed by animals exposed to heat for 3–4 hours with about 80% surviving the first day, which steadily declined to ~60% survival by day 4. It is of note that the exposure duration of this group coincides with the critical time point (3-4h) identified in the transcriptomic data. In the third group, with heat-exposures over 6 hours, survival chances were already drastically reduced after the first day (<20%).
Analogous to the 3 distinct survivorship groups found for short-, intermediate- and long-term stress exposures, this division was also present in the fraction of nematodes regaining a healthy movement during the recovery period, as well as regaining a normal number of progeny (Fig 3B and 3C). The movement in populations exposed to a short heat stress (< 3 hours) did not differ from that of control populations. While the heat stress initially causes slightly lower numbers of progeny, the reproduction peeked 2–3 days after the heat stress together with the control population. For intermediate exposure durations (3–4 hours), 60–70% of animals displayed a normal movement, yet reproduction was further reduced and delayed in these populations. For longer heat exposures (>4 hours), the few surviving individuals commonly presented abnormal, slow, and sporadic movement.
Overall, the transcriptional patterns during heat stress changed dramatically around 3–4 hours which coincided with a phenotypic change, as shown by the drastic decrease in movement, viable offspring, and recovery chances in the days after.
Heat stress disrupts major developmental processes
To investigate the correlation between gene expression and the different phenotypes, we first looked at how normal developmental processes progressed under heat stress conditions. Snoek et al. have dissected the temporal patterns of global transcript levels of C. elegans spanning the entire 4th larval stage [25], which corresponds to the time frame used in this study. We analysed the heat stress expression patterns in relation to developmental gene expression. First, we selected gene clusters strongly upregulated during L4 development at 20°C (see Materials and methods for details; Fig 4A). These genes showed little change in heat stress conditions at 35°C. Likewise, genes with a strong transcriptional response to heat stress (cluster 3) displayed few expression differences during development. Next, we selected gene clusters with a strong decrease in expression levels (Fig 4B). While most of the transcriptional patterns differed between development and heat stress conditions, a relatively small number of genes (i.e. 82 genes) were present in both groups. An enrichment analysis of these genes found a strong overrepresentation associated with the cuticle structure and locomotion.
The log2 transformed gene expression is indicated by the colour scale. Developmental gene expression data obtained from [25]. Time was measured beginning 46 hours post age-synchronization. The order of genes within each heat map was retrieved through hierarchical clustering, and is therefore not the same between the upper and lower panels. (A) Gene clusters with a strong up-regulation during development (left panel; cluster 1 in [25]) or during heat stress (right panel; cluster 3). (B) Gene clusters with a strong down-regulation during development (left panel; cluster 5 in [25]) or during heat stress (right panel; cluster 1). Venn diagrams presents the number of genes within each group.
Together, these results showed that heat stress disrupted the major transcriptional changes normally occurring during L4 development, indicative for the heat stress induced developmental delay. Furthermore, it shows that the animal almost fully shifts its transcriptional program to deal with the acute heat stress conditions.
To identify the genes involved in the sharp decrease in survival after the turning point in expression dynamics, we compared gene expression levels of samples taken between 2–4 hours into the heat-shock with the samples taken at the last three time points (6, 8 and 12 hours). We found 262 upregulated and 667 downregulated genes (q-value < 0.0001; S1 Table and S2 Fig). Enrichment analysis revealed an overrepresentation of genes involved with cuticle development and metabolic processes in the late heat stress down-regulated group. Genes involved in reproduction, development, and locomotion were enriched in the late heat stress up-regulated group, possibly showing the continuation of the normal developmental processes after the initial slow down. However, it should be noted that upregulation occurred with very low effect sizes, which could indicate a hampering of the transcriptional processes.
Discussion
By analysing a series of stress exposure times in C. elegans we detected a shift in gene expression patterns between 3–4 hours into the stress response, separating an initially highly dynamic phase from a later mostly stagnant phase. Survival, progeny count, and movement revealed that exposure to a heat stress lasting longer than 4 hours resulted in irreversible damage. Overall, the heat stress response could be divided into three distinct phases: i) an early highly dynamic phase up to 2 hours of stress exposure (including a very early upregulation of heat shock proteins), ii) an intermediate phase in which the transcriptional response attenuates presenting a turning point in dynamics (3–4 hours), iii) a late phase with gradual transcriptomic changes (6–12 hours). Phenotypically, each phase corresponds with distinct trends in the ability to recover from the stress, the ability to recover normal movement, and to produce viable offspring in the four days of recovery following the stress.
Gene expression regulation under stress conditions is strictly controlled, its kinetics are rapid and very often it is reversible. This allows for extremely rapid adaptation of cells and tissues in response to general stress, in particular heat stress, and for returning to a baseline level [23]. We analysed the phenotypic recovery from these rapid adaptive changes occurring during stress and found that already a relatively early response to heat stress abruptly changes development.
Early dynamic response to heat-stress disrupts development
During the early phase of heat stress (< 3 hours), the transcriptional response is highly dynamic. About 400 genes (cluster 2 and 5) are highly upregulated. Comparison with transcriptional patterns normally occurring during development has shown that this gene-set specifically reacts in the response to stress. Furthermore, genes highly active during development show low transcriptional changes during stress conditions. These results indicate that the animal almost fully switches its transcriptional focus on counteracting the adverse effects of heat stress. In this context, it is not surprising that heat shock proteins are the first set of genes to show a strong and rapid increase in transcript levels. Shortly after, histones and genes associated with the nucleosome assembly are highly enriched in upregulated gene clusters. Nucleosome remodelling has previously been shown to be an important part of the stress response, e.g. by allowing access to transcription sites of stress responsive genes [23,26,27]. In C. elegans, depletion of a nucleosome remodelling complex leads to a higher thermal sensitivity [20]. Packaging of DNA into nucleosomes could be an additional protective mechanism during the stress response.
After short exposures to stress, the animals recovered a healthy movement phenotype and started reproducing, indicating that the protective mechanisms put in place by the early transcriptional heat shock response are sufficient in this time frame. However, the disruption of normal transcriptional development could be one of the causes for the observed delay in reproduction. A study in which a two hour 35°C heat-shock was compared to two hour recovery of that heat-shock showed that the transcriptional patterns in the recovery population had still not returned to normal [17]. Also, a delay in reproduction has previously been shown in heat-shocked pre-gravid adult C. elegans exposed to temperatures between 30–32°C [15]. Arresting reproduction ensured limited damage to reproductive compartments during stress conditions. We found this delay on a transcriptional level in the early heat stress response, as development and reproduction related genes did not show their normal up regulation.
Attenuation of dynamic response
At medium exposures (3–4 hours), the transcriptional stress response attenuates. Initially, most animals (>90%) survive this exposure duration. However, in the four days following the stress, a noticeable ~40% decrease in survival and an increased occurrence of animals with an abnormal movement phenotype is observed. This shows that a change in transcriptional dynamics during stress could be a first indicator of a change in recoverability.
The attenuation of the heat-shock response has mostly been studied in several cell lines [20,21]. An integral part is the activation and subsequent suppression of the HSF-1 transcription factor activity through a negative feedback loop, which is partially mediated by those chaperones that are transcriptionally induced by HSF-1, such as HSP-70 [28]. The attenuation of the heat shock response is believed to serve a protective function, as cell lines with defects in the process display lower growth rates and reduced fitness [21,23]. In C. elegans, it was shown that a gain-of-function mutation in a negative regulator of the heat-shock response (HSB-1) results in severe effects on survival after heat stress [20]. In our data set, transcripts of chaperones induced by HSF-1 strongly increased within the first 30 minutes of the stress response (S1 Fig). The increase slows down until peak levels are reached at 4 hours into the stress response, followed by a small decrease and complete attenuation. It is unclear if the observed global transcriptional slowing down is due to an actively regulated process, such as the HSF-1 feed-back loop, or due to a passive process, such as the accumulation of damage to key cellular processes through enhanced protein aggregation, or limitations in protein synthesis.
Another factor to consider in the observed transcriptional dynamics and the corresponding decrease in recovery is the developmental stage at which the experiment was performed. It is well known that the age of the animals influences the stress response; e.g. the inducibility of heat shock response genes abruptly declines in very early adulthood [29], and phenotypically even different L4 sub-stages have different survival rates after heat-stress [10]. During unperturbed development of L4 larvae at 20°C, the C. elegans transcriptome is highly dynamic, marked by a pronounced shift at 50 hours, which overlaps in time with the turning point observed during stress [25]. Besides the possible heat-induced cell damage, passing this point of attenuation might result in irreversible changes to the pronounced developmental dynamics leading to a the strong decrease in recovery chances.
At long exposures, recovery chances are drastically reduced. While most transcript levels have reached a plateau, a distinct exception is the pronounced decrease in expression of a set of genes highly enriched with collagen related genes. Collagens are key components of the nematode cuticle, which is critical for protection and locomotion [30]. During development, the transcription of cuticle collagens is tightly regulated between any of the four molts [30,31]. Comparison with normal development (S3 Fig), shows a disruption of these patterns and a general downregulation of all cuticle genes. More recently, gene expression studies in C. elegans have shown that collagen genes are highly expressed in short heat stress exposure and during oxidative stress [18,32]. The strongly reduced survival chances after longer exposure might be caused by the later reduction of transcript levels of these cuticle genes in our experiment.
Overall, our study links a strong shift in transcriptional dynamics upon exposure to heat stress with an inability to recover from the stress response. The inability to recover was reflected in a decrease in worm activity, progeny count, and survival in the days after. Therefore, we suggest this critical shift in the dynamics of gene expression marks a point-of-no return ultimately leading to death.
Materials and methods
Nematode culturing
Hermaphrodites of the Caenorhabditis elegans strain Bristol N2 were used for all experiments and kept under standard culturing conditions at 20°C on Nematode Growth Medium (NGM) seeded with Escherichia coli strain OP50 as food source [33]. For the experiments, starved populations were placed onto fresh NGM dishes seeded with E. coli OP50 by transferring a piece of agar and subsequently grown at 20°C for 3–4 days until sufficient gravid adults had developed. Age-synchronized populations were obtained by bleaching according to standard protocols using a hypochlorite solution [33]. After bleaching, eggs were transferred to fresh 9 cm NGM dishes and maintained at 20°C.
Heat stress treatment
The heat stress treatments were performed in an incubator set to 35°C. N2 populations were exposed to the heat stress treatment starting 46 hours after age-synchronization during the L4 stage. We selected the L4 stage because nematodes in this stage exhibit a strong response to heat-stress [11]. The response declines in adult worms [29]. Samples were taken at several time points during the stress period: 0.5, 1, 2, 3, 4, 6, 8, or 12 hours. In total, 3–5 samples were collected for each time point. As preparation for the transcriptome analysis, the populations were washed off the plate with M9 buffer, collected in Eppendorf tubes and flash-frozen in liquid nitrogen and stored at -80°C until further use. For phenotypic observations, the N2 populations were transferred back to pre-heat stress maintenance conditions at 20°C.
Phenotyping
The selected traits (movement, survival, and progeny count) were observed using a stereomicroscope at approximately 24, 48, 72, and 96 hours post heat stress. To allow for accurate scoring of all individual animals, the population size per dish was kept at a maximum of 25 animals at the start of the experiment. In total, 3 dishes per heat stress duration were scored, which amounts to a total of approximately 60 animals per treatment. Animals were transferred to fresh NGM dishes every day during the reproductive phase using a platinum wire. Bagging and suicidal animals were censored.
Movement and survival.
Movement was scored based on classification systems that have previously been described in association with aging studies, where it acts as a measure of the biological age [34,35]. These systems were combined and adapted to score the impact of the heat stress. Healthy nematodes are actively moving in a sinusoidal pattern (Hosono: type I; Herndon: Class A). As a result of the heat stress, a proportion of the animals deviated from the healthy phenotype in varying degrees such as visibly lower levels of activity, low responsiveness to touch with a platinum wire and/or an irregular shape of movement (for example due to a partially paralysed tail). This is corresponding to Class B and C of Herndon or Type II and III of Hosono). Worms were scored as dead, when no head movement was observed after 3 touches with a platinum wire.
Progeny count.
It has previously been shown that C. elegans can lay non-viable eggs after heat stress [11]. For this reason, the progeny count was measured, defined as the absolute number of living offspring per population. We counted the progeny one day after transferring the adults of the experimental populations to fresh dishes, at which time viable eggs have hatched. For populations with a high level of reproduction, the total number of live offspring was estimated based on the count of a quarter of the dish.
Transcriptome profile
RNA isolation.
RNA was isolated from the flash frozen samples using the Maxwell® 16 AS2000 instrument with a Maxwell® 16 LEV simplyRNA Tissue Kit (both Promega Corporation, Madison, WI, USA). The mRNA was isolated according to protocol with a modified lysis step. Here, 200 μl homogenization buffer, 200 μl lysis buffer and 10 μl of a 20 mg/ml stock solution of proteinase K were added to each sample. The samples were then incubated for 10 minutes at 65°C and 1000 rpm in a Thermomixer (Eppendorf, Hamburg, Germany). After cooling on ice for 1 minute, the samples were pipetted into the cartridges resuming with the standard protocol.
Sample preparation and scanning.
For cDNA synthesis, labelling and the hybridization reaction, the ‘Two-Color Microarray-Based Gene Expression Analysis; Low Input Quick Amp Labeling’—protocol, version 6.0 from Agilent (Agilent Technologies, Santa Clara, CA, USA) was followed. The C. elegans (V2) Gene Expression Microarray 4X44K slides manufactured by Agilent were used. The microarrays were scanned with an Agilent High Resolution C Scanner using the settings as recommended. Data was extracted with the Agilent Feature Extraction Software (version 10.5) following the manufacturers’ guidelines.
Data pre-processing.
Data was analysed using the ‘R’ statistical programming software (version 3.3.2 64-bit). For normalization, the Limma package was used with the recommended settings for Agilent [36]. Normalization within and between arrays was done with the Loess and Quantile method, respectively [37]. The obtained normalized intensities were log2 transformed and outliers were removed. Batch effects within the data set were calculated with a linear model and removed as previously described [25]. For further analysis, the expression values of biological replicas were averaged. To analyse temporal expression dynamics independent from absolute expression values, the individual intensities measured at each time point for each gene were rescaled to the average expression in time per gene. The obtained values were log2 transformed and are further referred to as the log2 ratio with the mean.
Data accessibility.
The microarray datasets supporting this article have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5753.
Data analysis.
Principal component analysis (PCA) was performed on the log2 ratio with the mean to explore the source of underlying variation in gene expression. PCA scores of the first and second component were used to select genes with a significant contribution to the variation in expression dynamics. Selection was based on a significance level p < 0.01 in a linear model relating expression values with PCA scores for the separate components. To explore overall trends in gene expression dynamics, gene clusters were extracted by k-means with 200 iterations on 10 different starting sets. Six clusters were sufficient to visualise distinct patterns in gene expression changes.
Differently expressed genes between time-points were deducted by a linear model using the log2 expression of individual samples (3–5 samples per heat stress duration). In cases where multiple time points were compared, they were grouped into one factor, e.g. Group 1 (2, 3, 4 hours) vs Group 2 (6, 8, 12 hours). A high significance level of q < 0.0001 was chosen.
To extract genes with similar expression patterns to heat-shock proteins, we used spearman correlation analysis on the log2 ratio with the mean averaged for hsp-16.1, hsp-16.2, and hsp-16.41. Genes were selected with a log2 change >2 within the first 30 minutes of heat exposure.
The phenotypic data was analysed for differences between groups. For survival, log-rank test was used for pairwise comparison between treatments. For movement, the distribution between healthy and non-healthy movement (absolute numbers; total over three replicas) was compared pairwise between treatments separately for each day using a chi-squared test. The average number of live offspring per adult on day 1 post stress was compared to the control (i.e. 0 hour heat stress) with a two-sided t-test. This was done separately for each treatment conditions. All p-values were corrected for multiple testing by the Benjamini-Hochberg method.
Enrichment analysis.
To explore the biological functions associated with selected gene sets, we used the functional annotation tool provided by DAVID 6.8 [38,39]. For the enrichment analysis (functional annotation chart), settings were limited to Gene Ontology terms (GOTERM_BP_DIRECT, GOTERM_MF_DIRECT, GOTERM_CC_DIRECT).
Developmental data.
List of genes within developmental cluster 1 and 5 (strongly up- and downregulated, respectively) were obtained from Snoek et al. [25]. The normalized developmental expression data set was retrieved from WormQTL [40,41]. From the developmental time series, a subset of samples were selected that correspond to the heat stress time series (i.e.: of the developmental time series 46h, 47h, 48h, 49h, 50h, 52h, 54h, and 58h corresponding with the time points in the heat stress time series 0h, 1h, 2h, 3h, 4h, 6h, 8h, and 12h, respectively). Expression data of replicates was averaged. To compare expression dynamics between the time series obtained during development and in heat stress conditions, we selected the expression data of subsets of genes with strong expression patterns during development (cluster 1 and 5, Snoek et al., 2014) and heat stress (cluster 1 and 3, Fig 2, S1 Table). Heatmaps (R; package: gplots) of the log2 ratio with the mean are used to visualize the comparison of the expression dynamics during development and in heat stress conditions for each subset of genes.
Supporting information
S1 Fig. Temporal expression patterns of heat stress responsive genes during 12 hours of heat stress (35°C).
Genes were selected based on the information provided by the Gene Ontology database for the GO term ‘response to heat stress’ (GO:0009408, WormBase version 257). Expression levels of individual genes are presented as the Log2 ratio with the mean.
https://doi.org/10.1371/journal.pone.0189445.s001
(PDF)
S2 Fig. Volcano plot showing the difference per gene in log2 transformed gene expression levels between medium exposure durations (2, 3, and 4 hours) and long exposure durations (6, 8, and 12 hours).
P-values were inferred from a linear model comparing the two groups, and corrected for multiple testing by the Benjamini-Hochberg method. The red line indicates the selected significance level resulting in the selection of 262 upregulated (positive effect size) and 667 downregulated (negative effect size) genes.
https://doi.org/10.1371/journal.pone.0189445.s002
(PDF)
S3 Fig. Heatmap of temporal expression (log2 ratio with the mean) of structural constituents of the cuticle during development and under heat stress conditions.
Developmental gene expression data was obtained from Snoek et al. 2014. Genes were selected based on the information provided by the Gene Ontology database for the GO term ‘structural constituents of the cuticle’ (GO:0042302, WormBase version 257).
https://doi.org/10.1371/journal.pone.0189445.s003
(PDF)
S1 Table. Gene lists used for GO enrichment analysis, and detailed output of the enrichment analysis performed with the functional annotation tool provided by DAVID 6.8.
https://doi.org/10.1371/journal.pone.0189445.s004
(XLSX)
S2 Table. Survival after heat stress.
Approximately 100 L4-stage animals per heat stress duration (1h, 2h, 3h, 4h, 6h, 8h, and 12h) were scored directly after removal from 35°C, and after overnight recovery. This experiment was performed separately from the phenotypic recovery experiment, due to the chance of damaging the animals in an effort to provoke movement directly after the stress exposure.
https://doi.org/10.1371/journal.pone.0189445.s005
(XLSX)
S3 Table. Output of statistical test on phenotype observations.
Each pairwise comparison is given. The p-values reported for survival are based on a log-rank test, for movement on a chi-squared test, and for number of offspring on a t-test. The tables report the Benjamini-Hochberg corrected p-values.
https://doi.org/10.1371/journal.pone.0189445.s006
(XLSX)
References
- 1. Tyedmers J, Mogk A, Bukau B. Cellular strategies for controlling protein aggregation. Nat Rev Mol Cell Biol. Nature Publishing Group; 2010;11: 777–788. pmid:20944667
- 2. Lindquist S. The Heat-Shock Response. Ann Rev Biochem. 1986;55: 1151–91. pmid:2427013
- 3. David DC, Ollikainen N, Trinidad JC, Cary MP, Burlingame AL, Kenyon C. Widespread protein aggregation as an inherent part of aging in C. elegans. PLoS Biol. Public Library of Science; 2010;8: 47–48. pmid:20711477
- 4. Rodriguez M, Snoek LB, De Bono M, Kammenga JE. Worms under stress: C. elegans stress response and its relevance to complex human disease and aging. Trends Genet. Elsevier Ltd; 2013;29: 367–374. pmid:23428113
- 5. Ben-Zvi A, Miller EA, Morimoto RI. Collapse of proteostasis represents an early molecular event in Caenorhabditis elegans aging. Proc Natl Acad Sci U S A. National Academy of Sciences; 2009;106: 14914–14919. pmid:19706382
- 6. Kourtis N, Tavernarakis N. Cellular stress response pathways and ageing: intricate molecular relationships. EMBO J. Nature Publishing Group; 2011;30: 2520–31. pmid:21587205
- 7. Zhou KI, Pincus Z, Slack FJ. Longevity and stress in Caenorhabditis elegans. Aging (Albany NY). 2011;3: 733–753.
- 8. Lithgow GJ, Walker GA. Stress resistance as a determinate of C. elegans lifespan. Mech Ageing Dev. 2002;123: 765–771. pmid:11869734
- 9. Morley JF, Morimoto RI. Regulation of longevity in Caenorhabditis elegans by heat shock factor and molecular chaperones. Mol Biol Cell. American Society for Cell Biology; 2004;15: 657–64. pmid:14668486
- 10. Zevian SC, Yanowitz JL. Methodological considerations for heat shock of the nematode Caenorhabditis elegans. Methods. 2014;68: 450–457. pmid:24780523
- 11. Rodriguez M, Snoek LB, Riksen JAG, Bevers RP, Kammenga JE. Genetic variation for stress-response hormesis in C. elegans lifespan. Exp Gerontol. Elsevier Inc.; 2012;47: 581–587. pmid:22613270
- 12. Glauser DA, Chen WC, Agin R, Macinnis BL, Hellman AB, Garrity PA, et al. Heat avoidance is regulated by transient receptor potential (TRP) channels and a neuropeptide signaling pathway in Caenorhabditis elegans. Genetics. Genetics Society of America; 2011;188: 91–103. pmid:21368276
- 13. Cypser JR, Wu D, Park SK, Ishii T, Tedesco PM, Mendenhall AR, et al. Predicting longevity in C. elegans: Fertility, mobility and gene expression. Mech Ageing Dev. Elsevier Ireland Ltd; 2013;134: 291–297. pmid:23416266
- 14. McMullen PD, Aprison EZ, Winter PB, Amaral LAN, Morimoto RI, Ruvinsky I. Macro-level modeling of the response of c. elegans reproduction to chronic heat stress. Sporns O, editor. PLoS Comput Biol. Public Library of Science; 2012;8: e1002338. pmid:22291584
- 15. Aprison EZ, Ruvinsky I. Balanced trade-offs between alternative strategies shape the response of C. elegans reproduction to chronic heat stress. PLoS One. Public Library of Science; 2014;9: e105513. pmid:25165831
- 16. Stroustrup N, Ulmschneider BE, Nash ZM, López-Moyado IF, Apfeld J, Fontana W. The Caenorhabditis elegans Lifespan Machine. Nat Methods. 2013;10: 665–70. pmid:23666410
- 17. Snoek LB, Sterken MG, Bevers RPJ, Volkers RJM, Van’t Hof A, Brenchley R, et al. Contribution of trans regulatory eQTL to cryptic genetic variation in C. elegans. BMC Genomics. BioMed Central; 2017;18: 500. pmid:28662696
- 18. Brunquell J, Morris S, Lu Y, Cheng F, Westerheide SD. The genome-wide role of HSF-1 in the regulation of gene expression in Caenorhabditis elegans. BMC Genomics. BioMed Central; 2016;17: 559. pmid:27496166
- 19. Gasch AP, Spellman PT, Kao CM, Carmel-Harel O, Eisen MB, Storz G, et al. Genomic Expression Programs in the Response of Yeast Cells to Environmental Changes TL—11. Mol Biol Cell. 2000;11 VN-r: 4241–4257.
- 20. Satyal SH, Chen D, Fox SG, Kramer JM, Morimoto RI. Negative regulation of the heat shock transcriptional response by HSBP1. Genes Dev. Cold Spring Harbor Laboratory Press; 1998;12: 1962–1974.
- 21. Abravaya K, Phillips B, Morimoto RI. Attenuation of the heat shock response in HeLa cells is mediated by the release of bound heat shock transcription factor and is modulated by changes in growth and in heat shock temperatures. Genes Dev. 1991;5: 2117–2127. pmid:1936996
- 22. Truttmann MC, Zheng X, Hanke L, Damon JR, Grootveld M, Krakowiak J, et al. Unrestrained AMPylation targets cytosolic chaperones and activates the heat shock response. Proc Natl Acad Sci. 2017;114: E152–E160. pmid:28031489
- 23. de Nadal E, Ammerer G, Posas F. Controlling gene expression in response to stress. Nat Rev Genet. Nature Publishing Group; 2011;12: 833–45. pmid:22048664
- 24. Guhathakurta D, Palomar L, Stormo GD, Tedesco P, Johnson TE, Walker DW, et al. Identification of a novel cis-regulatory element involved in the heat shock response in Caenorhabditis elegans using microarray gene expression and computational methods. Genome Res. 2002;12: 701–712. pmid:11997337
- 25. Snoek LB, Sterken MG, Volkers RJM, Klatter M, Bosman KJ, Bevers RPJ, et al. A rapid and massive gene expression shift marking adolescent transition in C. elegans. Sci Rep. 2014;4: 3912. pmid:24468752
- 26. Guertin MJ, Petesch SJ, Zobeck KL, Min IM, Lis JT. Drosophila heat shock system as a general model to investigate transcriptional regulation. Cold Spring Harb Symp Quant Biol. 2010;75: 1–9. pmid:21467139
- 27. Shivaswamy S, Bhinge A, Zhao Y, Jones S, Hirst M, Iyer VR. Dynamic remodeling of individual nucleosomes across a eukaryotic genome in response to transcriptional perturbation. Rando OJ, editor. PLoS Biol. Public Library of Science; 2008;6: 0618–0630. pmid:18351804
- 28. Morimoto RI. Regulation of the heat shock transcriptional response: Cross talk between a family of heat shock factors, molecular chaperones, and negative regulators. Genes Dev. Cold Spring Harbor Laboratory Press; 1998;12: 3788–3796.
- 29. Labbadia J, Morimoto RI. Repression of the Heat Shock Response Is a Programmed Event at the Onset of Reproduction. Mol Cell. Elsevier Inc.; 2015;59: 639–650. pmid:26212459
- 30. Page A, Johnstone IL. The cuticle. WormBook. 2007; 1–15. pmid:18050497
- 31. Johnstone IL, Barry JD. Temporal reiteration of a precise gene expression pattern during nematode development. EMBO J. European Molecular Biology Organization; 1996;15: 3633–9.
- 32. Shin H, Lee H, Fejes AP, Baillie DL, Koo HS, Jones SJM. Gene expression profiling of oxidative stress response of C. elegans aging defective AMPK mutants using massively parallel transcriptome sequencing. BMC Res Notes. 2011;4: 34. pmid:21303547
- 33. Brenner S. The genetics of Caenorhabditis elegans. Genetics. 1974;77: 71–94. pmid:4366476
- 34. Hosono R, Sato Y, Aizawa SI, Mitsui Y. Age-dependent changes in mobility and separation of the nematode Caenorhabditis elegans. Exp Gerontol. Pergamon Press Ltd; 1980;15: 285–289.
- 35. Herndon LA, Schmeissner PJ, Dudaronek JM, Brown PA, Listner KM, Sakano Y, et al. Stochastic and genetic factors influence tissue-specific decline in ageing C. elegans. Nature. 2002;419: 808–814. pmid:12397350
- 36. Zahurak M, Parmigiani G, Yu W, Scharpf RB, Berman D, Schaeffer E, et al. Pre-processing Agilent microarray data. BMC Bioinformatics. 2007;8: 142. pmid:17472750
- 37. Smyth GK, Speed T. Normalization of cDNA microarray data. Methods. 2003;31: 265–273. pmid:14597310
- 38. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. Nature Publishing Group; 2008;4: 44–57. pmid:19131956
- 39. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. Oxford University Press; 2009;37: 1–13. pmid:19033363
- 40. Snoek LB, Van Der Velde KJ, Arends D, Li Y, Beyer A, Elvin M, et al. WormQTL-public archive and analysis web portal for natural variation data in Caenorhabditis spp. Nucleic Acids Res. 2013;41: 1–6.
- 41. Van Der Velde KJ, De Haan M, Zych K, Arends D, Snoek LB, Kammenga JE, et al. WormQTLHD—A web database for linking human disease to natural variation data in C. Elegans. Nucleic Acids Res. Oxford University Press; 2014;42: 1–8. pmid:24217915