Environmental and Evolutionary Drivers of the Modular Gene Regulatory Network Underlying Phenotypic Plasticity for Stress Resistance in the Nematode Caenorhabditis remanei

Abstract

Organisms can cope with stressful environments via a combination of phenotypic plasticity at the individual level and adaptation at the population level. Changes in gene expression can play an important role in both. Significant advances in our understanding of gene regulatory plasticity and evolution have come from comparative studies in the field and laboratory. Experimental evolution provides another powerful path by which to learn about how differential regulation of genes and pathways contributes to both acclimation and adaptation. Here we present results from one such study using the nematode Caenorhabditis remanei. We selected one set of lines to withstand heat stress and another oxidative stress. We then compared transcriptional responses to acute heat stress of both and an unselected control to the ancestral population using a weighted gene coexpression network analysis, finding that the transcriptional response is primarily dominated by a plastic response that is shared in the ancestor and all evolved populations. In addition, we identified several modules that respond to artificial selection by (1) changing the baseline level of expression, (2) altering the magnitude of the plastic response, or (3) a combination of the two. Our findings therefore reveal that while patterns of transcriptional response can be perturbed with short bouts of intense selection, the overall ancestral structure of transcriptional plasticity is largely maintained over time.

When faced with novel and stressful environmental conditions, individuals may acclimate to survive and populations can adapt to flourish (Hoffmann and Hercus 2000). Phenotypic plasticity is one mechanism of acclimation (Bradshaw 1965). Like other complex phenotypes, plasticity has a genetic basis and can evolve in response to selection (West-Eberhard 2003; Moczek et al. 2011). The adaptive response of a population to new, stressful conditions may therefore involve both acclimation and adaptation via plasticity evolution (Via and Lande 1985; Gomulkiewicz and Kirkpatrick 1992; Gavrilets and Scheiner 1993; Lande 2009; 2014). Furthermore, novel environments in the wild may present several stresses simultaneously that can lead to a correlated response in means of different phenotypes (e.g., Grant and Grant 1995; Fischer et al. 2007), and a shared genetic basis may lead to covariance of phenotypic plasticity of different traits across environments (e.g., Czesak et al. 2006; Stinchcombe et al. 2010).

Phenotypic plasticity has been studied in the laboratory and the field for more than a century (Baldwin 1896b; 1896a; Clausen et al. 1940; Waddington 1953; 1956; Schmitt et al. 1995; Bennett and Lenski 1997; DeWitt 1998; Nussey et al. 2005; Cheviron et al. 2013) and has been shown to be adaptive in many different systems (e.g., Dudley and Schmitt 1996; Agrawal 1998; Aubret et al. 2004; Charmantier et al. 2008). More recently, molecular analyses of gene expression have begun to identify systems that underlie phenotypic plasticity and its evolution (e.g., Gasch et al. 2000; Swindell et al. 2007; Hodgins-Davis and Townsend 2009; Badisco et al. 2011; Schunter et al. 2014; Alvarez et al. 2015). Gene expression levels can be characterized as a norm of reaction for a given genotype across environments (Figure 1). If populations of organisms continue to experience the new environment long enough for genetic changes to accumulate, then gene reaction norms can change both in basal expression levels and in the plasticity of expression itself (Figure 1).

Despite the advances in understanding gene regulatory differences that underlie phenotypic plasticity that have come from comparative analyses (Alvarez et al. 2015), the relative balance between rapid plastic responses in gene regulation and longer term changes in baseline expression in response to a novel environment are largely unknown. Experimental evolution provides a powerful approach to provide this understanding because proximal exposure to both novel and ancestral environments can be separated from long-term adaptation to the novel environment in a robust and reproducible fashion (Yampolsky et al. 2012; Sikkink et al. 2015; Huang and Agrawal 2016). When whole genome approaches (e.g., RNA-seq) are utilized, correlated changes in reaction norms of numerous genes provide a systems level way to tackle this question (Eisen et al. 1998; Wolfe et al. 2005; Promislow 2005; Barchuk et al. 2007; Rose et al. 2016).

We evolved phenotypic plasticity in populations of the nematode C. remanei in the laboratory by selecting for resistance to heat stress and oxidative stress (Sikkink et al. 2014b; Sikkink et al. 2015). In this paper we use a differential gene expression approach via RNA-seq to determine the structure and evolution of gene coexpression networks. Our goal was to address the relative balance between rapid plastic responses in gene regulation and longer term changes in baseline expression of a trait in response. In addition, we asked if adaptation to the heat and oxidative stresses involved the same or independent co-regulated modules. Our findings reveal that while patterns of transcriptional response can be perturbed with short bouts of intense selection, and that these molecular changes are at least partially independent across stressors, the overall ancestral structure of transcriptional plasticity is largely maintained over time.

Materials And Methods

Experimental evolution of C. remanei

We used the experimentally evolved populations of C. remanei that have previously been described (Sikkink et al. 2014b; Sikkink et al. 2014a; Sikkink et al. 2015). Briefly, 26 isofemale strains of C. remanei were isolated from terrestrial isopods (Family Oniscidea) collected from Koffler Scientific Reserve at Jokers Hill, King City, Toronto, Ontario. These strains were crossed in a controlled fashion to promote equal genetic contributions from all strains. The resulting genetically heterogeneous population (PX443) was the ancestral population for the experimental evolution.

A subset of the ancestral population was used for transcriptional profiling. In addition to the ancestor, three experimentally evolved populations were sampled for RNA-sequencing. All selection lines had been evolved at 20° as described previously (Sikkink et al. 2014b; Sikkink et al. 2015). One representative control population, one heat-selected population, and one oxidative-selected population were used. The heat-selected line was generated by exposing age-synchronized L1 larval worms to a 36.8° heat shock approximately every second generation. The oxidative-selected was similarly treated with a 1mM solution of hydrogen peroxide. The control populations received a mock selection treatment, from which worms were selected at random to continue the selected line at a similar census size. All lines were frozen after every two selection events. The final experimentally evolved populations used for the transcriptomics had experienced a total of 10 acute selection events and five freeze-thaw cycles.

Expression data

We collected L1 tissue from the ancestral, control, heat-selected, and oxidative-selected populations to use for transcriptional profiling (Supporting Information, Figure S1). All lines except the oxidative-selected population have been previously described (Sikkink et al. 2014b; Sikkink et al. 2014a). Briefly, we thawed frozen stocks of worms from each population. Except in the oxidative-selected population, 6 replicates per treatment were collected from a minimum of two independently thawed populations from each line. For the oxidative-selected line, all replicates were collected from a single thawed population of worms. Worms were raised at 20° until the population was large enough to collect enough individuals for RNA isolation. Age-synchronized L1 larvae were raised for 20 hr in liquid medium at either 20° or 30° (Figure S1). Prior to tissue collection, larval worms were passed through a 20-μm Nitex screen to remove unhatched eggs and dead adults. Total RNA was isolated from approximately 100,000 pooled individuals using standard TRIzol methods. Sequencing libraries were prepared according to the protocols as previously described (Sikkink et al. 2014b; Sikkink et al. 2014a). Samples were sequenced from a single end, to a length of 100 nucleotides in six lanes on an Illumina HiSeq 2000 at the University of Oregon Genomics Core Facility.

Initial processing of RNA-seq data

Initial quality filtering of raw sequence reads was performed using the process_shortreads component of the software Stacks (Catchen et al. 2011; 2013). Reads were discarded if they failed Illumina purity filters, contained uncalled bases, or if sample identity could not be determined due to sequencing errors in the barcode sequence. Reads with ambiguous barcodes were recovered if they had fewer than two mismatches from a known barcode. Using the alignment software GSNAP (Wu and Watanabe 2005; Wu and Nacu 2010), we aligned all reads that passed the quality filters to a reference genome for C. remanei assembled from strain PX356 (Fierst et al. 2015). We then used the htseq-count tool from the Python package HTSeq (Anders et al. 2015) to count all reads unambiguously aligning to gene models.

For all expression analyses, we first normalized the gene counts from all samples to account for differences in library size using the scaling procedure implemented in the DESeq2 package (Anders and Huber 2010; Love et al. 2014) in R version 3.4.0 (R Development Core Team 2017). The expression dataset was next filtered to exclude the 40% of genes with the lowest variance across treatments. Independent filtering of genes with very low variance in expression across treatments generally improves power in subsequent analyses (Bourgon et al. 2010; Anders et al. 2013). The expression values of the remaining genes were transformed using the variance-stabilizing transformation within DESeq2 for further analysis.

Finally, we used surrogate variable analysis (SVA) to reduce signal from batch effects and other unknown sources of variance (Leek and Storey 2007). SVA uses the residual expression matrix after accounting for the variables of interest—in this case temperature and population—to identify latent variables within the expression matrix. Using the R package sva (Leek and Storey 2007), we identified five latent variables encapsulating the unknown effects. To remove the effects of these variables from the gene expression data, we used limma (Ritchie et al. 2015) to build a regression model including only the latent variables obtained from SVA. We retained the residual expression matrix for all downstream analyses.

Multivariate analysis of transcriptional variation

We used non-metric multidimensional scaling (nMDS), which is an unsupervised ordination method that enables highly-dimensional data to be projected onto a few axes for visualization. For RNA-seq data, nMDS may be preferable as an ordination method, because it does not assume linear relationships within the data, enabling nMDS algorithms to robustly extract complex patterns from gene expression data (Taguchi and Oono 2005). One drawback of this nonparametric approach is, however, that the scores for variables mapped onto ordination axes can not be easily interpreted (in contrast to principal component scores, for example), and other methods may be required to identify genes contributing to differences between groups.

To carry out the nMDS ordination, a dissimilarity matrix was calculated for the filtered dataset of SVA residuals using Bray-Curtis dissimilarities (Bray and Curtis 1957). Using other distance metrics did not substantially alter the ordination plot. Data transformation, ordination, and scaling were performed in five dimensions using the vegan package (Oksanen et al. 2013). We tested for significant differences among populations and treatments using a permutational analysis of variance performed on the Bray-Curtis dissimilarity matrix. Population, treatment, and the interaction term were included as effects in the model, and 1000 permutations were run.

de novo Network Analysis

We used weighted gene coexpression network analysis (WGCNA), implemented in the R package WGCNA (Langfelder and Horvath 2008) to identify sets of genes (modules) that are highly correlated in their expression patterns. The initial network was constructed from all replicate samples for all treatments (n = 48) using a signed adjacency matrix with power = 5 to construct the topological overlap matrix. Hierarchical clustering of the topological overlap matrix was performed using the hclust function of flashClust (method = “average”) (Langfelder and Horvath 2012). Initial module assignments were made using the dynamicTreeCut algorithm (Langfelder et al. 2008) with the following options: cutHeight = 0.905, deepSplit = 2, minClustSize = 30, pamRespectsHybrid = FALSE.

We used a resampling approach to determine the probability that each gene was assigned to the appropriate module. To do this, we selected four of the six replicates for each treatment group at random to create a new subsampled dataset with 32 samples each. A total of 100 resampled datasets were created in this way. Using the same parameters as in the full dataset, we reconstructed the gene coexpression network for each of the resampled datasets. We then assessed whether each gene belonged in a given module by the following criteria. (1) If a module within a resampled network comprised at least 10% of a module in the full network, the genes in the resampled module were considered a significant group within the original module. Thus, each module from the full network could consist of several well-supported modules from the resampled data, but the network topology is biased toward the modules created from the full dataset. (2) Every gene from the resampled modules that was included in such a group in at least 70% of the resampled networks was determined to be strongly supported as a member of that module in the original network. All genes that did not meet these criteria were removed to the “unassigned” bin. After poorly supported genes were removed from each module, we merged highly correlated modules together based on the correlations among module eigengenes. Modules with highly correlated eigengenes (r > 0.9) were merged into single modules. Finally, the genes from any remaining modules that contained fewer than the minimum cluster size of 30 genes were moved to the “unassigned” bin.

Functional Annotation of Evolving Modules

The eigengene for each module, defined as the first principal component of the expression of all the genes in the module, was calculated to represent the general pattern of expression seen within each module. With the samples from the ancestral population, we first performed a t-test for significant differences in eigengene expression across the two environments to identify modules contributing to the ancestral response to thermal stress. We then performed an analysis of variance on module eigengenes to test for effects of population, temperature, and population-by-temperature interactions on the overall module expression. For modules that had a significant effect of population, we used Tukey HSD to identify pairwise differences between lines. All statistical tests on the eigengene expression were carried out in the R statistical environment (R Development Core Team 2017).

To functionally annotate the genes expressed in our data, we searched the nr protein database using BLAST+ (v2.2.28) (Camacho et al. 2009) to find significant (E-value > 10⁻³) matches to each gene. We used Blast2GO v3.0 (Conesa et al. 2005; Conesa and Götz 2008) to map significant BLASTx results to gene ontology terms and to compute a Fisher’s Exact Test to test for significant over-representation of GO terms in each module.

We examined each module for enrichment of known regulatory targets of 23 transcription factors for which binding data are available for the related nematode C. elegans. Binding targets for all transcription factors except for the FOXO transcription factor DAF-16 were obtained from the C. elegans modENCODE project (Niu et al. 2011). These targets were all identified from chromatin immunoprecipitation sequencing (ChIP-seq) of transgenic strains tagged with a dual GFP:3xFLAG tag and immunoprecipitated using anti-GFP antibodies (Niu et al. 2011). Putative target genes bound by DAF-16 have been previously identified using two different approaches: ChIP using anti-DAF-16 antibodies (Oh et al. 2006) and DNA adenine methyltransferase identification (DamID; Schuster et al. 2010). Target genes could be included multiple gene sets if they are bound by more than one transcription factor.

C. remanei homologs for each of the C. elegans transcription factor targets were determined based on the annotations that have been curated in the WS220 release of WormBase (Harris et al. 2009). Homologous genes identified by any method were included as possible transcription factor targets in C. remanei. In cases where multiple C. remanei genes were matched to a single gene in C. elegans, all possible homologous genes were included in the gene set, since no information was available to determine whether transcription factor binding was preserved preferentially in either possible homolog.

Modules were tested for significant enrichment of target genes bound by each transcription factor using a one-tailed Fisher’s exact test. In addition, we tested for enrichment of the C. remanei heat shock proteins previously identified (Sikkink et al. 2014b).

Data Availability

Sequence data for the ancestral, control, and heat populations were previously deposited in the NCBI Gene Expression Omnibus (GEO) database as part of series GSE56510 with accession numbers GSM1362987-1363022. All of this data, as well as additional RNA-seq data for the oxidative population, were realigned and reannotated and deposited in GEO under accession number GSE115496. Supplemental material are available at Figshare: https://doi.org/10.25387/g3.7361183.

Results

Transcriptional regulation divergence occurs both across temperatures and between evolved populations

We first sought to determine whether samples from the various populations or temperatures could be differentiated based on global patterns of gene expression. To do this, we used non-metric multidimensional scaling (nMDS), a powerful ordination method that does not assume linear relationships among variables (Taguchi and Oono 2005). On the first axis, we observed distinct separation between the two temperature treatments (Figure 2A). To assess the significance of the observed differences between temperature treatments, we used a permutational analysis of variance (PERMANOVA) on the dissimilarity matrix. This analysis confirmed that temperature had a highly significant effect on global gene expression (F_1,40 = 6.41, P = 0.001).

The four populations also differed significantly from one another (PERMANOVA; F_3,40 = 3.17, P = 0.001). The population differences accounted for variation observed on Axis 2 and Axis 3 in the nMDS analysis (Figure 2B). Both the oxidative and heat selected lines diverged from the ancestral population on Axis 2, but in opposite directions (Figure 2D). In contrast, the control population separated from the ancestral and selected populations on Axis 3 (Figure 2E). This pattern of divergence from the ancestor suggests that the three different selection regimes (two stresses and lab adaptation) lead to unique changes in transcriptome regulation in these populations. In addition, the response to the temperature treatment was strongly dependent on the population (PERMANOVA, population-by-temperature interaction: F_3,40 = 1.77, P = 0.002). The consequences of the interaction effect are most apparent on Axis 4 (Figure 2F), where the control-selected population responds to temperature in the opposite direction compared to the remaining lines, and Axis 5 (Figure 2G), on which the heat population responds in the opposite direction.

Network modules are differentially associated with line- and temperature-specific variation in expression

Because nMDS is a non-metric method, the contribution of specific genes (or suites of genes) to divergence on each axis is not readily interpretable. To address this limitation, we used weighted gene co-expression network analysis (Zhang and Horvath 2005; Langfelder and Horvath 2008) to identify modules, which are sets of genes with strongly correlated expression patterns that are more loosely connected to other such modules. We sought to identify modules that were important in the differential regulation of stress resistance in our evolved populations of C. remanei, because members of a gene module often share a common function (Eisen et al. 1998; Wolfe et al. 2005), and highly correlated gene sets may share transcriptional regulators (Allocco et al. 2004, but see also Marco et al. 2009).

Network analysis identified 13 co-expressed modules containing a total of 5,622 genes (Table 1). An additional 9,212 genes were not consistently assigned to any module after resampling and were designated as “Unassigned”. For each module, we calculated the eigengene (Langfelder and Horvath 2008), defined as the first principal component of the module. An eigengene’s expression explains the largest proportion of variance for the genes within the module and is therefore representative of the expression of the combined set of correlated genes within the module (Figure S2). For all assigned (numbered) modules, the eigengene explains more than 40% of the expression variance within the module, with several explaining 50 to 70% (Table 1). The second principal component explains no more than 5.6% in any assigned module (Figure S2), confirming that the eigengene is in fact an appropriate representation of expression for the genes within the assigned modules.

Eigengene analysis of modules identified in gene coexpression network analysis.

Despite evolution the ancestral structures of gene expression reaction norms are largely retained

In the balance between environment and evolution in changes in patterns of transcriptional regulation (Figures 1 and 3), plasticity and/or the retention of ancestral plasticity via the evolution of baseline expression predominates. Significant temperature differences or temperature-by-population interactions were observed in 12 out of 13 modules (Table 1). In fact, the two largest modules, Module 1 and Module 2, show effects only attributable to temperature (“shared plasticity”). Module 1 is composed of genes that are upregulated in all populations under immediate heat stress, whereas Module 2 contains genes that are downregulated in response to heat (Figure 3). The largest class of modules retained the overall pattern of plasticity displayed across populations, but changed their baseline level of expression within each environment (“divergent baseline”, Figure 3, Table 1). Finally, four of the modules showed heterogeneity among the pattern of plasticity that was dependent upon their specific evolutionary history (“evolved and/or divergent plasticity”, Figure 3). In one case (Module 5), the oxidative stress population lost its apparent response to heat stress entirely, indicative of genetic assimilation.

To more specifically address evolutionary divergence among lines, we examined pairwise differences among them for modules that showed a significant population effect (Figure 3). Focusing specifically within the ancestral population, five modules (Modules 1, 2, 3, 4, and 7), which together contain 4,836 genes, had significant expression differences across environments (Table 2). Evolved lines that diverged from the ancestral population are of particular interest, as these could indicate a set of genes that are adaptive for stress resistance. Two modules, Module 6 and Module 9, differed in the heat-selected population only. These modules are expected to contain genes that are important for adaptation to heat stress. Similarly, five modules (3, 4, 5, 7, and 8) are significantly different from the ancestor only in the oxidative-selected population. Two additional modules, Module 10 and Module 13, have evolved in both stress-selected populations. In both cases, the responses are in opposite directions in each stress-selected population (Figure 3), similar to the observed differences on Axis 2 in the nMDS analysis (Figure 2B).

Effect of temperature on module expression in the ancestral population

Surprisingly, among the unassigned genes, we observed a significant effect of both temperature and population in the eigengene expression, supporting the observation that a pattern of “divergent baseline” predominates the overall structure of the evolution of gene expression across the genome (Table 1). Given the conservative approach we used to assign genes into modules following resampling, it is likely that genes contributing to the unassigned modules were falsely removed from a real module. However the proportion of variance explained by the eigengene for the unassigned module is relatively small (9.3%). Therefore, it is unlikely that the lack of module assignment for these genes substantially alters the overall network topology we have observed.

Gene expression modules are each enriched for functionally related genes

We next examined the functional relationships among genes in identified modules by looking for enrichment of Gene Ontology terms within each module, especially terms in the biological process ontology (Figure 4). Enrichment of molecular function and cellular component terms are shown in Figures S3 and S4, respectively. Modules 1, 4, and 5, which contain genes upregulated in response to temperature, were enriched for terms falling under the GO categories pertaining to embryonic development, reproduction, and cellular transport, as well as metabolic processes relating to nucleic acid metabolism and translation. Module 2, which encompasses the large module of genes down-regulated in response to heat, was enriched for genes involved in cell-cell signaling. Modules 3 and 7, which had lower expression in the oxidative-selected population, were enriched for immune system genes and genes involved with translation, respectively. Modules 9 and 11 were enriched for genes involved in metabolic processes acting on molecules other than nucleic acids.

Regulatory targets of stress-responsive transcription factors are enriched in network modules

In C. elegans, several transcription factors are known to be critical regulators of cellular responses to stress. However, these regulators may not be differentially expressed in response to stress themselves, but rather undergo protein modifications to activate them under certain conditions. For example, the FOXO transcription factor DAF-16 is a major target of the insulin/insulin-like growth factor signaling (IIS) pathway in worms and is responsible for mediating responses to heat and oxidative stress, among others (Honda and Honda 1999; Hsu et al. 2003). DAF-16 is normally localized to the cytoplasm, but in stress conditions, DAF-16 is activated and transported to the nucleus, where it regulates transcription of many target genes (Lin et al. 2001; Lee et al. 2001). We identified C. remanei homologs of known binding targets of 23 transcription factors with previously published transcription factor binding profiles (Oh et al. 2006; Schuster et al. 2010; Niu et al. 2011), and tested for significant enrichment in each of the network modules. We also examined enrichment of another stress-related candidate gene set, the heat shock protein families (hsps) previously examined in Sikkink et al. (Sikkink et al. 2014b).

We observed significant enrichment (FDR < 0.05) of regulatory targets for all but four of the available transcription factors (Figure 5). Many of the transcription factors are key regulators of embryonic development. Unsurprisingly, Module 1 is enriched for transcriptional targets of most of these genes, including three HOX transcription factors—LIN-39, MBA-5, and EGL-5—consistent with the functional annotation of Module 1’s role in development.

Several transcription factors that regulate stress responses also showed enrichment of their target genes in one or more of these modules. PHA-4, a developmental regulator necessary for formation of the pharynx (Mango et al. 1994; Horner et al. 1998), has also been implicated in regulating heat shock response through HSP90 (van Oosten-Hawle et al. 2013). Targets of PHA-4 were enriched in Module 1 as well as the oxidative-evolved Modules 4 and 5. Genes regulated by SKN-1, another target of IIS that is critical for oxidative stress resistance (An and Blackwell 2003), were also enriched in Modules 1 and 4. Modules 2 and 3 were both enriched for targets of PQM-1, a C2H2 zinc finger and leucine zipper-containing protein (Tawe et al. 1998). In C. elegans, PQM-1 is responsive to certain types of oxidative stress (Tawe et al. 1998), and is a key regulatory target of IIS, in addition to DAF-16 (Tepper et al. 2013). Module 2 was also enriched for DAF-16 targets.

Module 9, which was significantly divergent in the heat-selected population, was significantly enriched for heat shock proteins. This module was also enriched for targets of ELT-3, a GATA transcription factor that functions during hypodermal development in C. elegans (Gilleard et al. 1999) and may also function downstream of IIS to influence longevity (Budovskaya et al. 2008), pathogen resistance (Pujol et al. 2008), and osmotic stress response (Rohlfing et al. 2010). ALR-1, a homeodomain transcription factor involved in development of sensory and GABAergic motor neurons (Tucker et al. 2005), and EOR-1 which regulates RAS/RAF-mediated signaling during development (Rocheleau et al. 2002; Howard and Sundaram 2002) also bind to more genes than expected within this module.

Modules 7, 8, 10 and 13 were also significantly different between the ancestor and at least one evolved population; however, they did not show significant enrichment of target genes for the available transcription factors. However, ChIP binding data from C. elegans was not yet available for some key transcription factors involved in stress response, particularly HSF-1 and HIF-1. Heat shock proteins are known to be regulated by HSF-1 in response to heat stress (Wu 1995; Åkerfelt et al. 2010), therefore enrichment of hsps in Module 9 may indicate a role for HSF-1 in regulation of that module.

Discussion

For many organisms, phenotypic plasticity is a vital adaptation that has evolved to cope with environmental stress. After almost a century of work on plasticity, researchers are now beginning to dissect the genetic regulatory networks that underlie a plastic phenotypic response (Promislow 2005; Barchuk et al. 2007; Rose et al. 2016; Nocedal et al. 2017). Here, we add to this growing body of work by examing global changes in gene expression in the evolution of phenotypically plastic responses using experimental evolution of C. remanei in response to the two related, but distinct, evolutionary stresses of heat and oxidative shock.

Global patterns of gene expression describe evolutionary divergence

We observed clear differentiation attributable to the induction of a response to temperature (i.e., plasticity), as well as evolved differences between populations (Figure 2). Exposure to the inducing temperature resulted in very pronounced changes in the global patterns of gene expression in all populations, which were primarily encapsulated by the primary axis resulting from nMDS analysis (Figure 2A). The scale of the observed response to the temperature shift in all populations (Figure 2C) largely confirms our previous observations for the heat selected population. Increased resistance to acute heat stress is not conferred by changing the degree to which genes respond to changes in the thermal environment (Sikkink et al. 2014b). However, the more powerful multivariate statistical framework used here reveals that the response to temperature changes did in fact differ among the evolved populations on secondary axes of ordination. These interaction effects were most evident on nMDS Axes 4 and 5 (Figure 2F and 2G), and suggest that some changes in transcriptional plasticity may be responsible for changes in phenotypic plasticity as well.

We also detected changes in gene regulation attributable to the evolutionary history of each line. Notably, the three selected populations have diverged from the ancestor in different directions on nMDS Axis 2 and Axis 3 (Figure 1B). This pattern suggests that at least partially independent axes of gene regulation contribute to adaptation in each case. These findings are consistent with the observations we have previously made—that there is no genetic correlation between heat and oxidative resistance under the environmental conditions in which these populations evolved (Sikkink et al. 2015). In short, although one might reasonably hypothesize a correlated selective response to heat and oxidative stresses that acts through a generic stress response pathway, our data support the alternative hypothesis. Evolution results from changes in different GRNs, or least different modules within a GRN, for these two related stresses.

Gene expression evolves in a modular fashion

The pattern of expression differences that we observed in our data indicates a high degree of modularity. Despite a relatively small number of experimental treatments, we were able to identify 13 transcriptional modules with highly correlated patterns of expression. Furthermore, the eigengenes that describe expression patterns within each module are differentially associated with the experimental treatments.

To test our ability to draw meaningful inferences from RNA-seq data we examined a well-known pathway, heat shock proteins (hsps), which are molecular chaperones known to be a critical component of response to heat stress (Lindquist and Craig 1988). Therefore, we expected these genes to form one or more modules that covary strongly with temperature. One module (9) seems to fulfill this expectation by capturing many of the elements of the hsp response. This module was strongly enriched for the set of heat shock proteins (Figure 5), and was also significantly regulated by temperature (Table 1).

Significant expression differences attributable to line were also observed within this heat shock response module (Figure 3). On closer examination, the heat-selected population was the only selected line to show divergence from the ancestral population. The eigengene expression of this module reveals that the heat-selected population has higher overall expression of these genes relative to the other populations even at 20° (Figure 3). However, these genes were still upregulated in response to temperature. The shift in expression in this module is consistent with our previous observation that the heat selected line evolved resistance to heat stress by shifting the thermal threshold for induction of the heat stress response (Sikkink et al. 2014b).

Gene expression module evolution is independent between stresses

Most of the identified modules, with the exception of Modules 1 and 2, show significant evidence for an evolved response in either of the stress-adapted populations (Table 1). Interestingly, the response in the heat and oxidative selected lines occurs primarily in independent modules, consistent with the lack of phenotypic correlated responses we have observed previously in this system (Sikkink et al. 2015).

Two modules, Module 10 and Module 13, do show evidence of evolutionary change in regulation in both the heat- and oxidative-selected lines (Table 1). However, our data do not support the evolution of a generalized stress response pathway contributing to adaptation in both of the populations that were evolved under these different stressors. In both modules, selection for heat resistance results in an overall change in gene expression in one direction, while selection for oxidative resistance occurs in the opposite direction (Figure 2). Presumably, these modules would contain a portion of the overlap in pathways expected based on the pleiotropy in the known stress response networks in C. elegans (Sikkink et al. 2015). However, these two modules together contain only 104 total genes, about 5% of the genes assigned to evolved modules, and the effect of selection to one stressor is antagonistic to the preferred response for the other stress.

The overall independence we observed in the genetic network between the evolved lines may contribute to the decoupling of the phenotypic responses to these different stressors, as certain modules within the greater response network have different degrees of pleiotropy, and can potentially allow for fine-tuning of the response.

Adaptation to different stressors involves non-overlapping subsets of the ancestral stress response network

Six modules, together comprising 4836 genes (about 86% of the genes assigned to modules) had significant responses to temperature in the ancestral population (Table 2). Most of these genes were in Modules 1 and 2, which show no evidence for evolved changes in any of the selective populations (Table 1). This suggests that the majority of the ancestral heat shock response remains unchanged under diverse selective scenarios. Several of the modules that have evolved only in the oxidative selective environment (Modules 3, 4, and 7) also had a significant response to temperature in the ancestral population (Tables 1 and 2). A reasonable interpretation of this pattern is that these modules contributed to the core stress response in the ancestral population, along with the genes in Modules 1 and 2. However, in the oxidative-selected populations, different subsets of genes are selected upon, enabling resistance to the oxidative stress instead of heat stress.

In contrast, in the modules that respond to heat selection (Modules 6, 9, 10, and 13), there was no strong evidence for plasticity in the ancestral population (Table 2). The genes in these modules are therefore unlikely to be involved in the ancestral plasticity in response to the inducing temperature, but rather are other sets of genes from a different regulatory cascade that have been modified under particular conditions. It is surprising that it is the heat-selected line that utilized novel stress response components, while oxidative selection in large part modified the existing stress response pathway. Because the eigengene is a combined metric of expression for multiple correlated genes, a trivial interpretation is that some individual genes within this module do in fact respond significantly to temperature in the ancestor, but we have limited power to detect the change when considered as a group. If this is true, then selection for heat stress may use components of the ancestral heat response, although these modules are still largely independent of those invoked for the oxidative stress response.

A more interesting explanation of the pattern is that the core heat stress response in the ancestor is already maximized and therefore that further adaptation requires cooption of additional genes that are not part of the core stress pathway. These interpretations may not be mutually exclusive, and determining the relative contributions of the ancestral vs. novel components of the stress response pathway in adaptation to heat stress will require additional investigation into the roles of individual genes within each population.

Regulation and function of the evolved gene expression plasticity modules

Genes that are co-regulated by a common transcription factor are likely to have highly correlated expression (Marco et al. 2009) and therefore should be classified as part of the same module. Identifying the transcriptional regulators of each module can provide important insight into which pathways contribute to the evolution of plasticity. In this study, we tested for enrichment of known targets of 23 transcription factors within each of the identified gene modules (Figure 4). Most of these transcription factors have vital roles in regulating developmental processes, but a few also have well-characterized roles in mediating stress responses to a variety of different stress types. It is important to note that the sets of target genes included here were initially identified in C. elegans and may not be comprehensive. However, the functions of these transcription factors are likely to be highly conserved in C. remanei, and these targets represent the best current hypothesis for transcription factor binding during Caenorhabditis development.

Many of the tested transcription factors are enriched in Module 1 (Figure 4). Given that the tested factors are key regulators of development, it is not surprising this module is also functionally annotated as involved in growth, embryonic development, and reproduction—organismal processes which are highly dependent on temperature in the poikilothermic C. elegans and its relatives (Riddle 1997). Consistent with that role, the genes in Module 1 do not appear to have evolved differences in gene expression between the ancestor and any of the evolved populations, and likely represent a core set of highly conserved developmental programs that are difficult to alter on short evolutionary timescales.

Not all of the components of these pathways are so highly conserved however. Module 4, for example, is also enriched for regulatory targets of many of the same transcription factors as Module 1, with the notable exception of the transcription factors required for neuronal developmental (Figure 4). This module had higher overall expression in the oxidative-selected population compared to the ancestor. It is reasonable to think that tissue-specific differences, and the combinatorial control of gene expression, together might allow for the increased evolvability of transcription in this module, despite the overlapping functional role of these genes and those in Module 1.

Another intriguing result from this study is the enrichment of PQM-1 targets in Modules 2 and 3. PQM-1, like DAF-16, is a major target of the IIS pathway, and the two transcription factors appear to function in opposition to one another (Tepper et al. 2013). Both transcription factors appear to act together to control a portion of the ancestral stress response identified in Module 2. However, Module 3 also shows evidence for evolved differences in expression in the oxidative-selected line (Table 1), and is enriched for targets of PQM-1 (Figure 4). PQM-1 is known to respond to oxidative stress, although previous studies describe a response to the oxidative stressor paraquat (methyl viologen; Tawe et al. 1998)) rather than the hydrogen peroxide used in this study. Further study will be required to determine whether PQM-1 and its targets are in fact important contributors to evolution of the oxidative stress response.

Conclusion

Using a powerful experimental evolution approach in C. remanei nematodes, we identified the evolution of transcriptional modules in response to selection in two stress environments. In general, patterns of ancestral plasticity dominated the evolutionary response, with the predominant mode of adaptation occurring through changes in baseline gene expression within an environment rather than changes in plasticity per se. The evolutionary responses to the two different stressors each modified unique modules, consistent with previous observations of low pleiotropy between these two responses (Sikkink et al. 2015). Surprisingly, the response in the population selected under heat stress conditions did not modify components of the ancestral response to heat shock, while the oxidative selection line did. We identified a number of modules with significant evolutionary and plastic responses that were enriched for targets of key developmental and stress response transcription factors. The scale of the total number of genes involved is daunting, however. Indeed, a major conclusion from this work is that short-term adaptation, although it may be built upon existing systems that facilitate a plastic response to the environment, can nonetheless alter the overall pattern of regulation of the majority of the genes in the genome (Boyle et al. 2017). Complex responses to environmental stress remain complex even when adapting to simplified selective environments.

Acknowledgments

We thank the many volunteers who have assisted with maintaining the C. remanei selection lines. This work was supported by a GRFP and DDIG (1210922) from the National Science Foundation to KLS, a Ruth L. Kirschstein NRSA Postdoctoral Fellowship to RMR (AG032900), National Institutes of Health grants to PCP (AG022500, AG043988, GM096008) and WAC (RR032670), and the Ellison Medical Foundation fellowship to PCP.

Footnotes

Literature Cited

Author notes