The distributions of species and their responses to climate change are in part determined by their thermal tolerances. However, little is known about how thermal tolerance evolves. To test whether evolutionary extension of thermal limits is accomplished through enhanced cellular stress response (enhanced response), constitutively elevated expression of protective genes (genetic assimilation) or a shift from damage resistance to passive mechanisms of thermal stability (tolerance), we conducted an analysis of the reactionome: the reaction norm for all genes in an organism’s transcriptome measured across an experimental gradient. We characterized thermal reactionomes of two common ant species in the eastern U.S, the northern cool-climate Aphaenogaster picea and the southern warm-climate Aphaenogaster carolinensis, across 12 temperatures that spanned their entire thermal breadth.
Results
We found that at least 2 % of all genes changed expression with temperature. The majority of upregulation was specific to exposure to low temperatures. The cool-adapted A. picea induced expression of more genes in response to extreme temperatures than did A. carolinensis, consistent with the enhanced response hypothesis. In contrast, under high temperatures the warm-adapted A. carolinensis downregulated many of the genes upregulated in A. picea, and required more extreme temperatures to induce down-regulation in gene expression, consistent with the tolerance hypothesis. We found no evidence for a trade-off between constitutive and inducible gene expression as predicted by the genetic assimilation hypothesis.
Conclusions
These results suggest that increases in upper thermal limits may require an evolutionary shift in response mechanism away from damage repair toward tolerance and prevention.
Temperature regulates biological activity and shapes diversity from molecular to macroecological scales [1, 2]. Many species, especially small-bodied arthropods, live at temperatures close to their thermal limits and are at risk from current increases in temperature [3–5]. Thermal tolerance, the ability of individuals to maintain function and survive thermal extremes, depends on a complex interplay between the structural integrity of cellular components and activation of physiological response mechanisms to prevent and/or repair damage [6, 7]. Thermal defense strategies are shaped by the environmental regime organisms experience [8] and thermal limits vary considerably among species and populations [3, 4, 9, 10]. These differences in thermal tolerance are largely genetic [11, 12] with a highly polygenic basis [13–16]. Outside of candidate genes [13], little is known about the evolution of thermal tolerance or the link between short-term physiological acclimation and longer-term adaptation to novel temperature regimes. This information is critical for understanding the adaptive potential of species to future climates [17].
To address this gap of knowledge, we need information on the extent to which selection has acted upon the diversity and plasticity of genes involved in thermal tolerance [17, 18]. In recent years, whole-organism gene expression approaches (e.g. transcriptomics) using high-throughput RNA sequencing (RNAseq) technology have been widely applied to identify genes involved in thermal tolerance [19–22] and other traits. Such studies typically use an ANOVA-type experimental or sampling design, with only a few environmental levels, and often find only a few dozen to hundred genes with differential expression in different thermal regimes. However, temperature and other environmental factors vary continuously in nature. As a result, such categorical comparisons (e.g. high vs. low temperatures) are likely to miss key differences that are due not just to whether it is hot, but rather how hot it is. Continuous variation is better characterized with a reaction norm approach, which describes variation in the phenotype of a single genotype across an environmental gradient [23]. Reaction norms differ not only in mean values, but also in their shapes [10, 24], and differences in the shape of reaction norms are often larger than differences in mean values at both the species and the population level [24].
In this study, we extend the reaction norm approach to RNAseq analysis and introduce the reactionome, which we define as a characterization of the reaction norm for all genes in an organism’s transcriptome across an environmental gradient. Although temporal patterns of transcriptional activity (e.g. fast- vs. slow- responding genes) are also important components of an organism’s transcriptional response to environmental conditions [25], we focus here on the response of transcripts across conditions at the same time point.
We use the reactionome method to identify genes that are thermally responsive in two closely-related eastern North American ant species, Aphaenogaster carolinensis and A. picea [26, 27]. Aphaenogaster are some of the most common ants in eastern North America [28] where they are keystone seed dispersers [29–31]. Ants, and ecotherms in general, have little or no thermal safety margin [5] and thus are highly susceptible to climate change [4, 32], putting at risk important ecosystem services [33]. Growth chamber studies have demonstrated that reproduction of Aphaenogaster will be compromised by increased temperatures [34], while field studies [32] and simulations [35] indicate that ant species persistence will depend on combinations of physiology and species interactions. Aphaenogaster carolinensis experiences a higher mean annual temperature (MAT) (14.6 °C) and less seasonal temperature variation (temperature seasonality = 7678°) than does A. picea (MAT = 4.6 °C, seasonality = 10,008°; [36]) at their respective collection sites. In controlled laboratory experiments, these warm- and cold-climate species exhibit corresponding differences in their critical maximum (44.7 °C for A. carolinensis versus 41.3 °C for A. picea; see Methods) and minimum temperatures (6.1 °C for A. carolinensis versus −0.1 °C for A. picea). These differences between species in their thermal environments and physiological tolerances allowed us to investigate adaptation to both lower and upper thermal extremes in this system.
To characterize the thermal reactionome, we measured the reaction norm for each gene using a regression-based statistical approach to identify temperature-dependent patterns of change in gene expression. We used these response patterns to quantitatively test three mechanistic hypotheses of thermal adaptation. First, the enhanced response hypothesis [37–39] proposes that species extend their thermal limits through a stronger induced response to provide greater protection from more frequently encountered stressors. This hypothesis would predict that the cool-adapted A. picea would activate more genes, and induce them more strongly, in response to low temperatures than would the warm-adapted A. carolinensis, which would show greater induction in response to high temperatures.
Second, the tolerance hypothesis [9, 40] proposes that existing inducible stress responses become insufficient or prohibitively costly as environmental stressors increase in frequency, resulting in a shift away from induced response in favor of structural changes [41] or behavioral adaptations [5, 42]. This hypothesis predicts adaptation to stress should be associated with lower transcriptional responsiveness and less sensitivity to temperature perturbation, as well as a shift to an alternate suite of tolerance genes and pathways [43, 44].
Finally, the genetic assimilation hypothesis [45, 46] proposes that exposure to more extreme stressors selects for a shift from inducible to constitutive expression of stress-response genes. This hypothesis predicts that transcripts responsive to high temperatures in A. picea will have higher constitutive expression in A. carolinensis, whereas transcripts responsive to low temperatures in A. carolinensis will have higher constitutive expression in A. picea. Of course, these hypotheses are not mutually exclusive. By taking a reactionome approach, we can quantify if, and under what conditions, these mechanisms contribute to thermal adaptation.
To summarise, in this project we generated the transcriptomes of two closely-related temperate ant species, and quantified their gene expression across a wide range of thermal conditions. We then evaluated three non-mutually exclusive hypotheses (enhanced response, tolerance and genetic assimilation) of the evolution of thermal adaptation by comparing the number and expression patterns of transcripts between species in response to extreme low and extreme high temperatures. Finally, we used gene ontology information to determine which gene products and pathways are involved in thermal adaptation in the two species.
Results
Reaction norms of thermally-responsive transcripts
The combined Aphaenogaster transcriptome assembly contained 99,861 transcripts. About half of these (51,246) transcripts had a signficant BLAST hit, of which 50 % (25,797) had a top hit to Insecta and 37 % (18,854) had a top hit to Formicidae. We performed a BUSCO analysis [47] to assess the quality of the transcriptome assembly against the arthropod Benchmarking Universal Single-Copy Orthologs (BUSCOs). This analysis revealed that the transcriptome is largely complete, as we recovered 1,426 complete single-copy BUSCOs (62%) and an additional 435 fragmented BUSCOs (16%), which is in line with the results of Simao et al [47] for transcriptomes of other non-model species. Moreover, only 8% of the BUSCOs were found to be duplicated in the combined transcriptome, which indicates that the steps (see Methods) we took to collapse homologs between the two species in the transcriptome assembly were successful.
We quantified gene expression using the program Sailfish [48], and fitted polynomial regression models to the expression values of each transcript to identify those that had a linear or quadratic relationship. To account for multiple tests, we both applied a False Discovery Rate (FDR) correction, and performed a resampling analysis to determine the number of transcripts that would be expected to have a significant relationshp by chance alone. We retained the 2509 (2.5 % of total) true positive transcripts that exceeded null expectation from the resampling analysis for further analysis (Additional file 1: Table S1). Of these transcripts, 75 % (1553) had a non-linear relationship with temperature that would likely have been missed with a standard differential expression experiment (e.g. high vs. low temperature). The proportion of responsive transcripts is similar if we focus only on those transcript with a BLAST hit (725 significant transcripts out of 51,246). However, as with all de novo transcriptome assemblies, this assembly is fragmented due to partial contigs and alternative transcripts [49] so this estimate is likely a lower bound for the true proportion of transcripts that are thermally responsive.
We used the predicted transcript expression levels to partition transcripts for each species into five expression categories (Fig. 1) which were defined a priori to allow us to test predictions derived from three thermal adaptation hypotheses of relative response severity in the two species: High transcripts had greatest expression at temperatures > 31 °C, Low transcripts had greatest expression at temperatures < 10 °C, Intermediate transcripts had greatest expression between 10 to 30 °C, Bimodal transcripts had increased expression at both high and low temperatures, while NotResp transcripts were those that were not thermally responsive in the focal species but did respond in the other.
Fig. 1
Illustrations of the pattern in relation to temperature for each of the four expression categories; Bimodal, High, Intermediate and Low. The fifth category of Not Responsive (a flat line) is not shown
Expression response to thermal extremes differs between species
Although the total number of thermally-responsive transcripts did not differ between species (χ21 = 0.08, P = 0.77), the two species differed in the number of transcripts in each expression category (Table 1, χ24 = 302.896, P < 0.001). Aphaenogaster picea induced significantly more transcripts in response to both temperature extremes (Bimodal transcripts in Table 1; χ21 = 71.617, P < 0.001) than did A. carolinensis, which downregulated more transcripts under these conditions (Intermediate transcripts in Table 1; χ21 = 256.329, P < 0.001). Consistent with the enhanced response hypothesis, the cool-climate A. picea induced 273 (~50 %) more transcripts in response to low temperatures than the warm-climate A. carolinensis (Low transcripts in Table 1; χ21 = 71.227, P < 0.001). However, there was no difference among species in the number of transcripts upregulated at high temperatures (High transcripts in Table 1; χ21 = 0.53, P = 0.47).
Table 1
Table of the number of thermally responsive transcripts by expression type for A. carolinensis and A. picea
Low
Intermediate
High
Bimodal
NotResp
A. picea
1193
249
248
278
110
A. carolinensis
920
680
232
117
129
Low are transcripts with increased expression at low temperatures (<10 °C), Intermediate are transcripts with maximum expression between 10–30 °C, High are transcripts with increased expression at high temperatures (>31 °C), Bimodal are transcripts with increased expression at both low and high temperatures, while NotResp are transcripts that are not thermally responsive in one species but are in the other species
In addition, we also examined the specific patterns of shifts from one expression category to another between species. As transcripts may change expression between species due to neutral drift alone, we used the Stuart-Maxwell test of marginal homogeneity to test if the number of responsive transcripts in each expression category differed between the species when controlling for overall differences in the number of responsive transcripts. We found that the expression categories of individual transcripts between the two species were not randomly distributed (Stuart-Maxwell test of marginal homogeneity χ24 = 319, P < 0.001, Additional file 2: Figure S1). Specifically, the two species differed significantly in expression pattern, which captures differences in slope as well as category, for 1553 (62 %) of the thermally responsive transcripts. Most of these shifted from one to another expression category.
The enhanced response and tolerance hypotheses make opposing predictions concerning the overlap in response patterns between the two species (Fig. 2). The enhanced response hypothesis posits that temperature adaptation uses existing mechanisms for thermal resistance, which should result in significant overlap in response and fewer transcripts shifting expression categories than expected by chance (Fig. 2, left). In contrast, the tolerance hypothesis predicts that transcripts involved in active defense will become non-responsive or shift to other expression categories in the better-adapted species (Fig. 2, right). We tested these predictions by examining if the transcripts upregulated in response to the temperature extreme experienced less frequently by a species (cool temperatures for the warm-climate A. carolinensis, and warm temperatures for the cool-climate A. picea) displayed the same response profile in the species that more frequently experiences those conditions.
Fig. 2
Illustrations of the expected thermal response patterns in the two species under alternative mechanistic hypotheses of temperature adaptation. Although both temperature extremes were investigated in a similar way, for simplicity only the response to low temperatures is illustrated here. Each column indicates the distribution across all response categories in A. picea, which has a lower CTmin and is therefore better adapted to low temperatures, for the set of transcripts identified as cold-induced (either High or Bimodal categories) in the species with higher CTmin , A. carolinensis, relative to the null hypothesis of equal marginal frequencies. The dashed boxes highlight cells that would indicate matched responses in the two species, and the color of each cell (blue = excess, orange = deficit) represents the deviation of the observed from expected number of transcripts. The Enhanced Response hypothesis (left) proposes that the increase in cold tolerance in A. picea is achieved by amplifying existing molecular mechanisms, and thus there should be an excess of shared response types between species. In contrast, the Tolerance Hypothesis (right) predicts that A. picea is less reliant on induced responses to confer cold-tolerance than A. carolinensis, leading to an excess of shifts from induction in A. carolinensis to the Not Responsive or down-regulation categories in A. picea
Transcripts upregulated at low temperatures in A. carolinensis (Low and Bimodal transcripts) were significantly biased toward this same category and away from other expression categories in A. picea (Fig. 3a), suggesting shared response pathways as predicted by the enhanced response hypothesis. In contrast, transcripts upregulated in response to high temperatures in A. picea (High and Bimodal) shifted expression categories in A. carolinensis (Fig. 3b), primarily to the Intermediate category (Fig. 3b). These transcripts are less likely to be upregulated in any context, consistent with the tolerance hypothesis.
Fig. 3
Results of analysis of thermal response patterns in the two species. The color of each cell (blue = excess, orange = deficit) represents the deviation of the observed from the expected number of transcripts based on hypothetical equivalence of the marginal frequencies. The units are number of transcripts. For each temperature extreme, the species expected to be less well adapated to that extreme is displayed on the x-axis for the two response categories corresponding to upregulation (Bimodal and Low for the low temperatures, or Bimodal and High for high temperatures). The distribution of response categories for those transcripts in the better-adapted species is arrayed along the y-axis. The dashed boxes indicate the matched responses (e.g. High - High). a Low temperature extreme: there is an excess of shared Low and Bimodal expression types and a bias away from all other categories in A. picea, consistentwith the enhanced response hypothesis (Fig. 1). b High temperature extreme: in addition to an excess of matched categories, there is an excess of High and Bimodal transcripts in A. picea that are not upregulated in A. carolinensis (Intermediate and Not Responsive), partially consistent with the tolerance hypothesis. The complete set of matched observations is shown in Additional file 1: Figure S1. Expression types are defined in Table 1
Molecular processes suggest a generalized stress response mechanism
The gene set enrichment analysis revealed a number of gene groups enriched in each expression category (Additional file 3: Table S2). Across both species, there were 9 terms enriched in the Bimodal category, including terms involved in stress response (regulation of cellular response to stress, signal transduction by p53 class mediator), cell death (apoptotic signaling pathway) and cellular organization (e.g. protein complex localization). The 6 terms enriched in the Low category suggest that proteins undergo structural (e.g. protein acylation) and organizational (single-organism organelle organization) changes to tolerate colder temperatures, possibly to maintain membrane fluidity [50]. The High category included only a single enriched GO term, “nicotinamide metabolic process”, while the Intermediate category had 5 terms including DNA packaging and metabolic process terms.
A. carolinensis has greater inertia of expression change to increases in temperature than does A. picea
As an additional test of the tolerance hypothesis, we examined the critical temperature of gene induction in response to increasing and decreasing temperatures. We compared between species the mean temperatures of transcript upregulation, defined as the temperature at which the transcript showed the greatest positive change in expression. In support of the enhanced response but not the tolerance hypothesis, the temperature of induction at low temperatures was significantly higher for the cool-climate A. picea than for A. carolinensis (12.4 °C) than A. picea (13.1 °C; t1308 = −3.1, P < 0.002; Fig. 4a), though the temperature of induction did not differ between species for high temperatures (t567 = 0.8, P < 0.403).
Fig. 4
Histogram with smooth density estimate of temperature of maximum rate of change in expression for transcripts that have (a) increased expression at Low and High temperatures and (b) decreased expression at Low and High temperatures. Red bars and lines are for A. carolinensis while blue bars and lines are for A. picea
Similarly, for down-regulated (Intermediate) transcripts, we compared the mean temperatures of downregulation of transcript expression between species at both high (>20 °C) and low (<20 °C) temperatures. Consistent with the tolerance hypothesis, A. carolinensis had greater inertia of gene expression in response to increasing temperatures. The temperature of downregulation for Intermediate transcripts was 28.6 °C for A. carolinensis compared to 27.2 for A. picea (t294 = 3.8, P < 0.001). The difference between species was not significant with decreasing temperatures (t251 = 0.5, P = 0.584, Fig. 4b).
No evidence for genetic assimilation
We tested the genetic assimilation hypothesis by comparing the log ratios of relative inducibility to relative baseline expression at the rearing temperature (25 °C). If stress-response transcripts have shifted between species from inducible to constitutive expression, there should be a negative relationship between the two. We found no evidence of such a relationship for either temperature extreme: transcripts more upregulated at high temperatures in the cool-climate A. picea were not expressed at higher baseline levels in the warm-climate A. carolinensis (Fig. 5a). Similarly, transcripts more upregulated at low temperatures in A. carolinensis did not show higher baseline levels in A. picea (Fig. 5b). In fact, for both comparisons we found a weakly positive relationship between relative inducibility and baseline expression between the two species (β1 = 0.31, P < 0.001 and (β1 = 0.21, P < 0.001). In addition, the thermally responsive transcripts in A. carolinensis, regardless of expression pattern, had higher baseline expression than those in A. picea, including those with Intermediate expression profiles in both species (Wilcoxon V = 68842, P < 0.001). An important exception to this pattern is the set of transcripts that had High or Bimodal expression in A. picea but were not thermally responsive in A. carolinensis (top-row of Fig. 3b). These transcripts are less likely to be upregulated in any context, consistent with the tolerance hypothesis.
Fig. 5
Scatterplots of log ratios of relative inducibility to relative constitutive expression, defined as expression level at the common rearing temperature (25 °C) for (a) High transcripts in A. picea (P < 0.001, r2 = 0.07) and (b) Low transcripts in A. carolinensis (P < 0.001, r2 = 0.1). Blue lines and confidence intervals are from ordinary least squares regressions
Discussion
The potential for many species to persist in face of climate change will depend in part upon their thermal tolerances. However, for most species little is known about how plasticity or adaptive changes in gene expression underlie thermal tolerance. By using a reactionome approach, we were able to quantitatively describe plasticity in transcript expression across a thermal gradient, and identify putative changes in gene expression associated with shifts in thermal tolerance between the ant species Aphaenogaster picea and A. carolinensis. We found non-linear patterns of gene expression changes in response to temperature, with both quantitative and qualitative differences between species, consistent with different mechanisms of thermal adaptation to low and high temperature extremes.
Under the enhanced response hypothesis, stress-adapted species are hypothesized to induce a stronger and earlier response to extreme conditions. We found evidence for this hypothesis at low temperatures: although the lower thermal limit for A. picea is substantially lower than A. carolinensis, A. picea upregulated responsive transcripts at slightly less extreme temperatures (Fig. 4a). Moreover, the transcripts upregulated in A. picea included about about half (55 %) those upregulated in A. carolinensis as well as an additional set of 261 transcripts (Table 1), enriched for metabolism, organization and translation processes (Additional file 3: Table S2). Two non-mutually exclusive hypotheses may explain this pattern. First, surviving prolonged low temperatures, such as would be experienced during overwintering, generally requires advance production of specialized cryoprotectants [43] and a suite of preparatory physiological modifications [51]. The northern species A. picea may induce a greater response to survive the longer winter period. Alternatively, the response to low temperatures may reflect countergradient expression to counteract reduction in enzyme efficiency, and maintain activity as temperature declines [41]. This requirement may be under stronger selection in A. picea given the shorter growing season that would necessitate foraging under a broader range of temperatures.
In contrast to cold tolerance, the enhanced upper thermal limit in A. carolinensis is best explained by the tolerance hypothesis. High temperatures were associated with significantly fewer upregulated transcripts in A. carolinensis (Table 1), and a large proportion (25 %) of the transcripts upregulated at high temperatures in A. picea were either downregulated or expressed at negligible levels overall in A. carolinensis. These results suggest that mechanisms other than the heat shock response are acting to maintain protein stability in face of temperature increases. Such mechanisms may include novel constitutive defenses [19, 21, 22], enhanced proteome stability [52] or behavioral quiescence [5] to tolerate thermal stress. These differences are in line with expectations that A. carolinensis, with a growing season over twice the length of its northern congener, may be better able to afford to restrain from foraging in suboptimal conditions. Indeed, quiescence under stressful conditions by the red harvester ant Pogonomyrmex barbatus has been shown to increase colony fitness [42].
The one hypothesis that did not receive support was the genetic assimilation hypothesis (Fig. 5), which predicts that exposure to more frequent stressors will select for a shift from inducible to constitutive expression of stress-response transcripts. This contrasts with other recent studies on adaptation in field populations to thermal stress [21]. However, in a short-term selection experiment for heat tolerance, Sikkink et al. [46] also found no evidence for genetic assimilation at the expression level after ten generations of selection for heat tolerance in Caenorhabditis remanei, even though there was a substantial increase in heat tolerance. Both the genetic assimilation and tolerance routes to increasing thermal limits are functionally similar in that they emphasize damage prevention rather than repair. Whether a particular taxon evolves one strategy over another may be related to availability of alternative mechanisms as well as the the intensity, frequency and duration of temperature stress in a given environment.
Given the differences in the patterns of thermal responsiveness between species (Fig. 3), it is worth noting a number of similarities. In both species, there were 2 -- 3 times more transcripts upregulated at low than high temperatures (Table 1). The degree of upregulation at low temperatures is surprising given previous studies [53, 54] that found little transcriptional activity at low temperatures. However, these studies exposed organisms to a few extreme (−10 -- 0 °C) temperatures. At these extremes, we also found few upregulated transcripts (Fig. 3a), whereas the peak of low-temperature transcriptional activation occurred near 10 °C (Fig. 4). A potential explanation for this pattern is that increased gene expression functions to support elevated metabolism at moderately cold temperatures, as suggested by the metabolic cold adaptation hypothesis [55]. The observation that more transcripts were upregulated at low than high temperatures could also be due to stronger selection on upper than lower thermal limits, thereby reducing both genetic variation and gene expression plasticity at high temperatures [4, 56]. This explanation is consistent with the observation in Aphaenogaster rudis [57] and other ectotherms [10, 58] that critical maximum temperatures vary less among taxa than do critical minimum temperatures.
Critical maximum and minimum temperatures are hypothesized to be genetically correlated [10, 58], but this was not evident in terms of gene expression in this study. Only ~10 % of transcripts upregulated in response to temperature were bimodal, and for both activation and down-regulation, thresholds differed between species at only one temperature extreme (Fig. 4). This suggests that species do not face a fundamental trade-off between these two limits and may be able to shift upper and lower thermal limits independently to match requirements of more seasonally variable environments.
A major contribution of this study is the construction of a reactionome for gene expression data. Similar approaches have been used in other species [59, 60], but to our knowledge, none have applied a regression approach to identify a complete list of responsive transcripts across an environmental gradient. This approach revealed quantitative patterns of temperature response not captured in categorical comparisons. For example, the degree of upregulation at cool (~10 °C) but not extreme cold temperatures was missed in previous studies that focused on extreme cold limits, as discussed above. Further, a number of issues have hampered RNA-seq studies to date. Namely, lists of differentially expressed transcripts are prone to false positives [61], depend on the genetic background of the organism [62] and are prone to “storytelling” interpretations [63]. Our findings are robust to these issues as we focus on the average change in the shape of the reaction norms across many hundreds of responsive transcripts in each species. Although we use gene ontology information to interpret our results, the key findings about differential plasticity of expression between species do not depend on functional annotation.
Moreover, by characterizing responses across thousands of transcripts, the reactionome approach can help to distinguish selection from neutral drift in gene expression [64–66]. Although we cannot rule out drift as a source of variation for individual transcripts, we would not expect to see systematic differences in expression type categories or critical temperature thresholds as we do here (Fig. 3, Additional file 2: Figure S1). Thus, our method provides an example of how focusing on transcriptome-wide changes in gene expression—as opposed to identifying lists of differentially-expressed transcripts—can provide meaningful insight on the process of evolution. It should be noted, however, that although including non-linear relationships between expression and temperature captured a significantly larger range of biologically-relevant responses, it also led to a substantial increase in false positives. Empirical estimation of these rates via randomization tests, combined with robust sampling designs, can help to minimize this bias and focus results on biologically-meaningful gene sets.
A number of caveats do apply to our work. First, species may differ in gene expression along axes which we have not measured here, especially temporal patterns of gene expression [25], which could be studied in further work. Second, the de novo transcriptome assembly is highly fragmented, given that all sequenced ant genomes to date have only about 18,000 genes [67]. Although we took steps to remove contaminants and redundant transcripts, some likely remain, in addition to partially assembled transcripts. A genome assembly, in progress, will help to reduce fragmentation. Third, the quality of the annotation for a non-model system such as Aphaenogaster is not as good as it would be for model arthropods such as Drosophila and Apis. Finally, the mapping of changes in gene expression to organismal fitness is far from direct [68], and large differences in patterns of gene expression may have only small effects on fitness. In particular, functional protein levels cannot be expected to be fully linked to mRNA abundance due to post-transcriptional modification, regulation, mRNA fluctuations and protein stability [68].
Our results are congruent with evidence from other systems [21] that thermally-stressful habitats select for investment in tolerance, whereas organisms from less stressful environments rely on plastically-induced resistance. Although the heat-shock response is one of the most conserved across living organisms [39], it is energetically expensive, particularly under chronic stress conditions [69]. Under such circumstances, it may be advantageous to proactively prevent thermal damage even at the cost of reduced metabolic efficiency, either by maintaining a higher constitutive level of chaperone proteins [11] or by increasing the thermal stability of proteins at the expense of catalytic activity [70]. Thus, although in the short term increasing temperature stress leads to a quantitatively stronger induced response, adapting to such stress over evolutionary time appears to require a qualitative shift in mechanism of resistance that can alter not only the magnitude, but the sign of gene expression change in response to temperature. Whether such a shift would be possible in the compressed time frame of projected climate change, particularly for long-lived organisms such as ants, is likely to be critical in determining the capacity of populations to adapt to more frequent and long-lasting stressors.
Conclusions
In this work, we have brought reaction norms to the genomic era by characterizing the thermal reactionomes of two temperate ant species, Aphaenogaster picea and A. carolinensis. At least 2 % of their transcriptomes are thermally responsive. Our results indicate that these two ant species have different responses to thermal extremes. A. picea responds by increasing expression of transcripts related to metabolism, stress response and other protective molecules, whereas A. carolinensis decreases expression of transcripts related to metabolism and likely relies on other mechanisms for thermal tolerance. The thermal reactionomes of these two species provide key insights into the genetic basis of thermal tolerance, and a resource for the future study of ecological adaptation in ant species. Finally, the reactionome itself illustrates a new direction for characterizing acclimation and adaptation in a changing climate.
Methods
Samples
Ants of the genus Aphaenogaster are some of the most abundant in eastern North America [71], and species as well as populations within species differ in critical maximum and minimum temperatures [57]. Temperature is a potentially strong selective force for ground-nesting ant populations, which must tolerate seasonally freezing winters and hot summers. On shorter time scales, individual workers can experience extreme thermal environments when they leave the thermally buffered ant nest to forage for food [32].
In fall 2012, we collected a single colony of Aphaenogaster picea from Molly Bog, Vermont (University of Vermont Natural Areas; 44.508° N, −72.702° W) and a single colony of Aphaenogaster carolinensis, part of the A. rudis species complex [26], from Durham, North Carolina (36.037° N −78.874° W). These sites are centrally located within each species’ geographic range. Species identity was confirmed with morphological characters (Bernice DeMarco, Michigan State University). Colonies of both species were maintained in common conditions at 25 °C for 6 months prior to experimentation. Due to colony size limitations, we were unable to determine the critical thermal limits of these particular colonies. In summer 2013 we collected additional colonies of Aphaenogaster from Molly Bog, VT and North Carolina (Duke Forest, 36.036° N, 79.077° W). We tested the upper and lower critical thermal limits for five ants from each of these colonies using a ramp of 1 °C per minute, starting at 30 °C, and recorded the temperature at which the ants were no longer able to right themselves, following the protocol of Warren & Chick [57].
Common garden design
Ideally, genetically-based variation in gene expression profiles would be identified by comparing individuals completely reared under common-garden conditions to eliminate environmental variation experienced either as adults or during development. However, Aphaenogaster colonies are long-lived, cannot be bred under laboratory conditions, and do not achieve complete turnover of the workforce for at least a year or longer. Thus, as is commonly done with other long-lived organisms [21, 65], we exposed both colonies to common-garden rearing conditions for six months to fully acclimate adult workers to common temperatures. Over this time, roughly 1–2 cohorts of new workers are expected to join each colony (~1/3 of the total), such that the workers sampled for thermal traits and gene expression are likely to have included a mix of adult-acclimated and fully lab-reared individuals.
Unlike ANOVA-based experimental designs, which derive statistical power from replication within each experimental treatment level, regression designs have greater power when sampling additional values across the range of the continuous predictor variable [72]. Ideally, the treatments should be replicated at each level of the predictor variable [73]. However, even with no replication, the regression design is still more powerful than an ANOVA design with comparable replication, and provides an unbiased estimator of the slope [72]. For these reasons, we focused our sequencing efforts on maximizing the number of temperatures at which the transcriptome was profiled, rather than on replication at each temperature.
To limit differences in gene expression not related to the experimental treatment, on 12 different days we haphazardly collected three ants from each 2012 colony at the same time of day to minimize variation due to circadian oscillations. We measured response to temperature with a one-hour static temperature application, which is ecologically relevant for workers that leave the thermally-buffered nest and are immediately exposed to ambient temperatures while foraging [71]. Each day, the ants were placed in glass tubes immersed in a water bath maintained at one of 12 randomly-assigned temperatures (0° to 38.5 °C, in 3.5° increments) for one hour. The minimum and maximum temperatures were selected based on previous work showing that these temperatures are close to the critical minimum (~0 °C) and maximum (~43 °C) temperatures for Aphaenogaster [57], and these treatments did not cause mortality. At the end of the hour, the ants were flash frozen in liquid nitrogen and stored at −80 °C. Thus, our reactionome characterized early, but not late, responding genes. We extracted mRNA by homogenizing the three pooled ants in 500 uL of RNAzol buffer with zirconium silicate beads in a Bullet Blender (Next Advance; Averill Park, NY), followed by RNAzol extraction (Molecular Research Center Inc; Cincinnati, OH) and then an RNeasy micro extraction (Qiagen Inc; Valencia, CA) following the manufacturer’s instructions.
Sequencing, assembly and annotation
For each species, the 12 samples were barcoded and sequenced in a single lane of 2 × 100 bp paired-end reads on an Illumina HiSeq 1500 yielding 200 and 160 million reads for the A. picea and A. carolinensis samples respectively. Reads were filtered to remove Illumina adapter sequences and low quality bases using the program Trimmomatic [74].
We assembled the sequenced reads into the full set of mRNA transcripts, the transcriptome, for the combined data set from both species using the Trinity de novo transcriptome assembly program [75]. De novo transcriptome assembly is prone to falsely identifying alternative transcripts and identifying inaccurate transcripts that are chimeric (e.g. regions of two separate transcripts that assemble into a false, or chimeric, third transcript) [76]. We removed potentially false transcripts by first running the program CAP3 [77] to cluster sequences with greater than 90 % similarity and merge transcripts with overlaps longer than 100 bp and 98 % similar in length. Second, we ran the program uclust which clusters sequences completely contained within longer sequences at greater than 90 % similarity (see Additional file 4). We used liberal values (90 % similarity) to merge orthologous transcripts in the two species that may not have assembled together in the initial de novo transcriptome assembly. To identify contaminant sequences, we screened our full transcriptome using the program DeconSeq [78] with the provided bacteria, virus, archaen and human databases of contaminants.
The Trinity de novo transcriptome assembly for both species assembled together included 126,172 transcripts with a total length of 100 million bp. Filtering to remove redundant or chimeric reads resulted in an assembly with 105,536 transcripts. The total length was 63 million bp with an N50 length of 895 bp and a mean transcript size of 593 bp. Of the 105,536 filtered transcripts, 55,432 had hits to the NCBI-nr database. Of these, 38,711 transcripts mapped to GO terms, 1659 transcripts were identified to an enzyme and 18,935 transcripts mapped to a domain with >50 % coverage. We removed 5675 transcripts identified as known contaminants, leaving 99,861 clean transcripts.
We assessed the quality of the transcriptome using the BUSCO program [47] available from (http://busco.ezlab.org/). BUSCO asseses transcriptome completeness by measuring the number of near-universal single-copy orthologs using the Arthropod database from OrthoDB.
To determine the putative function of the transcripts, we used functional annotation of the transcriptome assembly using the web-based tool FastAnnotator [79] which annotates and classifies transcripts by Gene Ontology (GO) term assignment, enzyme identification and domain identification.
Identification of thermally-responsive transcripts
We quantified expression of each transcript using the program Sailfish [48] and used the bias-corrected transcripts per million (TPM) [80] as our measure of transcript expression. We included the contaminant transcripts identified by DeconSeq at the quantification stage to avoid incorrectly assigning reads to other transcripts, but removed these from further analyses. Because preliminary examination of the data (Additional file 4) indicated that the 7 °C samples may have been mis-labeled, we omitted these data from the analysis. The expression values were highly correlated between species at each temperature treatment (r2 > 0.98) indicating that assembling the transcriptome with data from both species was justified (Additional file 3).
To identify transcripts that had significant changes in expression across the thermal gradient, we fit to each transcript an ordinary least-squares polynomial regression model
Temperature and species were both fixed effects, with a quadratic term included for temperature. We used log(TPM + 1) as the response to control for skew in the expression data. For a continuous predictor such as temperature, this regression approach is preferred to an ANOVA approach as it can reveal non-linear responses such as hump-shaped or threshold effects [72]. This method is robust to over-dispersion because we expect errors in the read count distribution [81] to be independent with respect to temperature.
To evaluate the statistical significance of the patterns, we computed parametric P-values for each model and adjusted these P-values using the False Discovery Rate (FDR) approach of Benjamini and Hochberg [82]. As a more stringent filter for false positives, we then randomly re-assigned each transcript within a species to a different temperature, fit the polynomial models as above, and again calculated P-values and FDR. Ideally, these randomized data sets should not yield any significant associations. We repeated this resampling approach 100 times, and used the 95th quantile of false significant transcripts as the null expectation for retaining transcripts from the true data.
Of these overall significant transcripts, we identified thermally-responsive transcripts as the subset that had significant β2(temp), β3(temp2),β4(species * temp) or β5(species * temp) terms after step-wise model selection by AIC. For each thermally-responsive transcript, we predicted expression levels using the final linear model for each species across the tested thermal range. We used the predicted transcript expression levels to partition transcripts for each species into five expression categories: High transcripts had greatest expression at temperatures > 31 °C, Low transcripts had greatest expression at temperatures < 10 °C, Intermediate transcripts had greatest expression between 10 to 30 °C, Bimodal transcripts had increased expression at both high and low temperatures, while NotResp transcripts were those that were not thermally responsive in the focal species but did respond in the other. For the Bimodal group, we required that expression at both low and high temperatures was at least one standard deviation greater than the expression at the rearing temperature of 25 °C. Because expression category was defined by the temperature of maximal expression, both Low and High categories were biased toward transcripts up-regulated at that temperature extreme, but also likely included some transcripts down-regulated at the opposing extreme. The two categories which could unambiguously distinguish up- from down-regulation are Bimodal (up at both extremes) and Intermediate (down at both extremes).
Statistical analyses
We used χ2 tests to determine if the total number of responsive transcripts, and the number of transcripts in each expression category differed between species. To evaluate if shifts from one expression category to another between the two species were randomly distributed, we used the Stuart-Maxwell test of marginal homogeneity from the coin package [83] in R [84] which tests if the row and column marginal proportions are in equity.
To test whether the temperature at which thermally-responsive transcripts were activated differs between species, we identified the temperature at which there was the greatest change in expression for each transcript in each species, using only the transcripts with a significant species x temperature interaction. For upregulated transcripts, we grouped the High transcripts along with the high temperature end of the Bimodal transcripts, and did the same for Low transcripts. We then performed a t-test to determine if the mean temperature of transcript activation differed between the two species for each group. For downregulated transcripts (i.e. Intermediate), we identified the greatest change in expression for each transcript in response to both increasing (>20 °C) and decreasing (<20 °C) temperatures, and used a t-test to compare the mean temperature of down-regulation between species.
To test for a tradeoff between induciblity and constitutive baseline expression between species, we fit ordinary least squares regressions with the log ratio of relative constitutive expression as the response variable and the log ratio of relative inducibility as the predictor variable for High transcripts in A. picea and for Low transcripts in A. carolinensis. Constitutive expression was defined as predicted expression at 25 °C, whereas inducibility of each transcript was defined as ((maximum TPM - minimum TPM)/minimum TPM) x 100. In addition, we used a Mann–Whitney test to compare the baseline constitutive expression between species for all responsive transcripts.
Gene set enrichment analysis
To describe the molecular processes involved in thermal adaptation, we performed gene set enrichment analysis (GSEA) using the parentChild algorithm [85] from the package topGO [86] in R [84]. Briefly, this approach identifies GO terms that are overrepresented in the significant transcripts relative to all GO terms in the transcriptome, after accounting for dependencies among the GO terms.
All analyses were performed with R 3.2 [84] and are fully reproducible (Additional file 4).
The Illumina short-read sequence data supporting the results of this article are available in the NCBI Short Read Archive BioProject repository, PRJNA260626 http://www.ncbi.nlm.nih.gov/bioproject/PRJNA260626/.
The Trinity transcriptome assembly, FastAnnotator annotation file and Sailfish gene expression quantification files supporting the results of this article are available from the LTER data portal, datasets hf113-38, hf113-41, and hf113-42 (http://dx.doi.org/https://doi.org/10.6073/pasta/05ea6464df30efa2f1e2c7439366bf47).
Declarations
Acknowledgements
Support made possible by NSF DEB Award #1136644 and the Vermont Genetics Network through Grant Number 8P20GM103449 from the INBRE Program of the National Institute of General Medical Sciences (NIGMS) of the National Institutes of Health (NIH). Its contents are solely the responsibility of the authors and do not necessarily represent the official views of NIGMS or NIH. We thank Lauren Nichols, Bernice DeMarco and Clint Penick, for assistance, and four anonymous reviewers for their feedback.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
(1)
Department of Biology, University of Vermont
(2)
Department of Ecology and Evolutionary Biology, University of Tennessee
(3)
Vermont Genetics Network, University of Vermont
(4)
Department of Biological Sciences, North Carolina State University
(5)
Harvard Forest, Harvard University
(6)
Center for Macroecology, Evolution and Climate, University of Copenhagen
(7)
Data Scientist
References
Brown JH, Gillooly JF, Allen AP, Savage VM, West GB. Toward a metabolic theory of ecology. Ecology. 2004;85:1771–89.View ArticleGoogle Scholar
Deutsch CA, Tewksbury JJ, Huey RB, Sheldon KS, Ghalambor CK, Haak DC, et al. Impacts of climate warming on terrestrial ectotherms across latitude. Proc Natl Acad Sci. 2008;105:6668–72.View ArticlePubMedPubMed CentralGoogle Scholar
Kingsolver JG, Diamond SE, Buckley LB. Heat stress and the fitness consequences of climate change for terrestrial ectotherms. Funct Ecol. 2013;27:1415–23.View ArticleGoogle Scholar
Sunday JM, Bates AE, Kearney MR, Colwell RK, Dulvy NK, Longino JT, Huey RB. Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation. Proc Natl Acad Sci. 2014;111:5610–5615.View ArticlePubMedPubMed CentralGoogle Scholar
Huey RB, Kingsolver JG. Evolution of thermal sensitivity of ectotherm performance. Trends Ecol Evol. 1989;4:131–5.View ArticlePubMedGoogle Scholar
Richter K, Haslbeck M, Buchner J. The heat shock response: Life on the verge of death. Mol Cell. 2010;40:253–66.View ArticlePubMedGoogle Scholar
Angilletta MJ, Wilson RS, Navas CA, James RS. Tradeoffs and the evolution of thermal reaction norms. Trends Ecol Evol. 2003;18:234–40.View ArticleGoogle Scholar
Hoffmann AA, Chown SL, Clusella-Trullas S. Upper thermal limits in terrestrial ectotherms: How constrained are they? Funct Ecol. 2013;27:934–49.View ArticleGoogle Scholar
Krebs R, Loeschcke V. Estimating heritability in a threshold trait: Heat-shock tolerance in Drosophila buzzatii. Heredity. 1997;79:252–9.View ArticlePubMedGoogle Scholar
Kellermann V, Overgaard J, Hoffmann AA, Fløjgaard C, Svenning J-C, Loeschcke V. Upper thermal limits of Drosophila are linked to species distributions and strongly constrained phylogenetically. Proc Natl Acad Sci. 2012;109:16228–33.View ArticlePubMedPubMed CentralGoogle Scholar
Krebs RA, Feder ME, Lee J. Heritability of expression of the 70KD heat-shock protein in Drosophila melanogaster and its relevance to the evolution of thermotolerance. Evolution. 1998;52:841–7.View ArticleGoogle Scholar
Williams BR, Van Heerwaarden B, Dowling DK, Sgrò CM. A multivariate test of evolutionary constraints for thermal tolerance in Drosophila melanogaster. J Evol Biol. 2012;25:1415–26.View ArticlePubMedGoogle Scholar
Morgan TJ, Mackay TFC. Quantitative trait loci for thermotolerance phenotypes in drosophila melanogaster. Heredity. 2006;96:232–42.View ArticlePubMedGoogle Scholar
Takahashi KH, Okada Y, Teramura K. Genome-wide deficiency screen for the genomic regions responsible for heat resistance in Drosophila melanogaster. BMC Genet. 2011;12:57.View ArticlePubMedPubMed CentralGoogle Scholar
Hoffmann AA, Willi Y. Detecting genetic responses to environmental change. Nat Rev Genet. 2008;9:421–32.View ArticlePubMedGoogle Scholar
Somero GN. Comparative physiology: A “crystal ball” for predicting consequences of global change. Am J Physiol Regul Integr Comp Physiol. 2011;301:R1–14.View ArticlePubMedGoogle Scholar
Meyer E, Aglyamova GV, Matz MV. Profiling gene expression responses of coral larvae (Acropora millepora) to elevated temperature and settlement inducers using a novel RNA-seq procedure. Mol Ecol. 2011;20:3599–616.PubMedGoogle Scholar
Teets NM, Peyton JT, Colinet H, Renault D, Kelley JL, Kawarasaki Y, et al. Gene expression changes governing extreme dehydration tolerance in an antarctic insect. Proc Natl Acad Sci. 2012;109:20744–9.View ArticlePubMedPubMed CentralGoogle Scholar
Barshis DJ, Ladner JT, Oliver TA, Seneca FO, Traylor-Knowles N, Palumbi SR. Genomic basis for coral resilience to climate change. Proc Natl Acad Sci. 2013;110:1387–92.View ArticlePubMedPubMed CentralGoogle Scholar
O’Neil ST, Dzurisin JDK, Williams CM, Lobo NF, Higgins JK, Deines JM, et al. Gene expression in closely related species mirrors local adaptation: Consequences for responses to a warming world. Mol Ecol. 2014;23:2686–98.View ArticlePubMedGoogle Scholar
Gomulkiewicz R, Kirkpatrick M. Quantitative genetics and the evolution of reaction norms. Evolution. 1992;46:390–411.View ArticleGoogle Scholar
Murren CJ, Maclean HJ, Diamond SE, Steiner UK, Heskel MA, Handelsman CA, et al. Evolutionary change in continuous reaction norms. Am Nat. 2014;183:453–67.View ArticlePubMedGoogle Scholar
Sørensen JG, Nielsen MM, Kruhøffer M, Justesen J, Loeschcke V. Full genome gene expression analysis of the heat stress response in drosophila melanogaster. Cell Stress Chaperones. 2005;10:312–28.View ArticlePubMedPubMed CentralGoogle Scholar
Umphrey G. Morphometric discrimination among sibling species in the fulva - rudis - texana complex of the ant genus aphaenogaster. Can J Zool. 1996;74:528–59.View ArticleGoogle Scholar
De Marco B, Cognato A. A multiple-gene phylogeny reveals polyphyly among eastern North American Aphaenogaster species. Zoologica. 2015. doi:https://doi.org/10.1111/zsc.12168.
King JR, Warren RJ, Bradford MA. Social insects dominate eastern US temperate hardwood forest macroinvertebrate communities in warmer regions. PLoS ONE. 2013;8:e75843.View ArticlePubMedPubMed CentralGoogle Scholar
Ness JH, Morin DF, Giladi I. Uncommon specialization in a mutualism between a temperate herbaceous plant guild and an ant: Are aphaenogaster ants keystone mutualists? Oikos. 2009;118:1793–804.View ArticleGoogle Scholar
Zelikova TJ, Sanders NJ, Dunn RR. The mixed effects of experimental ant removal on seedling distribution, belowground invertebrates, and soil nutrients. Ecosphere. 2011;2:art63.View ArticleGoogle Scholar
Rodriguez-Cabal MA, Stuble KL, Guénard B, Dunn RR, Sanders NJ. Disruption of ant-seed dispersal mutualisms by the invasive asian needle ant (pachycondyla chinensis). Biol Invasions. 2012;14:557–65. View ArticleGoogle Scholar
Diamond SE, Nichols LM, McCoy N, Hirsch C, Pelini SL, Sanders NJ, Ellison AM, Gotelli NJ, Dunn RR. A physiological trait-based approach to predicting the responses of species to experimental climate warming. Ecology, 2012;93:2305–2312.View ArticlePubMedGoogle Scholar
Toro ID, Ribbons RR, Pelini SL. The little things that run the world revisited: A review of ant-mediated ecosystem services and disservices (Hymenoptera: Formicidae). Myrmecological News. 2012;17:133–146.Google Scholar
Pelini SL, Diamond SE, Maclean HJ, Ellison AM, Gotelli NJ, Sanders NJ, Dunn RR. Common garden experiments reveal uncommon responses across temperatures, locations, and species of ants. Ecology and Evolution. 2012;2:3009–15.View ArticlePubMedPubMed CentralGoogle Scholar
Sharon B, Stuble KL, Lessard J-P, Dunn RR, Adler FR, Sanders NJ. Predicting future coexistence in a North American ant community. Ecology and Evolution. 2014;4:1804–1819.View ArticleGoogle Scholar
Hijmans R, Cameron S, Parra J, Jones P, Jarvis A. Very high resolution interpolated climate surfaces of global land areas. Int J Climatol. 2005;25:1965–78.View ArticleGoogle Scholar
Hofmann GE, Somero GN. Interspecific variation in thermal denaturation of proteins in the congeneric mussels mytilus trossulus and m. galloprovincialis: Evidence from the heat-shock response and protein ubiquitination. Mar Biol. 1996;126:65–75.View ArticleGoogle Scholar
Feder ME, Hofmann GE. Heat-shock proteins, molecular chaperones, and the stress response: Evolutionary and ecological physiology. Annu Rev Physiol. 1999;61:243–82.View ArticlePubMedGoogle Scholar
Kültz D. Molecular and evolutionary basis of the cellular stress response. Annu Rev Physiol. 2005;67:225–57.View ArticlePubMedGoogle Scholar
Fields PA. Protein function at thermal extremes: Balancing stability and flexibility. Comp Biochem Physiol A Mol Integr Physiol. 2001;129:417–31.View ArticlePubMedGoogle Scholar
Lockwood BL, Somero GN. Functional determinants of temperature adaptation in enzymes of cold- versus warm-adapted mussels (genus mytilus). Mol Biol Evol. 2012;29:3061–70.View ArticlePubMedGoogle Scholar
Gordon DM. The rewards of restraint in the collective regulation of foraging by harvester ant colonies. Nature. 2013;498:91–3.View ArticlePubMedGoogle Scholar
Neelakanta G, Hudson AM, Sultana H, Cooley L, Fikrig E. Expression of ixodes scapularis antifreeze glycoprotein enhances cold tolerance in drosophila melanogaster. PLoS ONE. 2012;7:e33447.View ArticlePubMedPubMed CentralGoogle Scholar
Franssen SU, Bergmann N, Winters G, Klostermeier UC, Rosenstiel P, Bornberg-Bauer E, et al. Transcriptomic resilience to global warming in the seagrass zostera marina, a marine foundation species. Proc Natl Acad Sci. 2011;108:19276–81.View ArticlePubMedPubMed CentralGoogle Scholar
Waddington C. Genetic assimilation of an acquired character. Evolution. 1953;7:118–26.View ArticleGoogle Scholar
Sikkink KL, Reynolds RM, Ituarte CM, Cresko WA, Phillips PC. Rapid evolution of phenotypic plasticity and shifting thresholds of genetic assimilation in the nematode Caenorhabditis remanei. G3 (Bethesda). 2014;4:1103–12.View ArticleGoogle Scholar
Patro R, Mount SM, Kingsford C: Sailfish enables alignment-free isoform quantification from RNA-seq reads using lightweight algorithms. Nature Biotechnology 2014, 32:462–464Google Scholar
Vijay N, Poelstra JW, Künstner A, Wolf JBW. Challenges and strategies in transcriptome assembly and differential gene expression quantification. a comprehensive in silico assessment of RNA-seq experiments. Mol Ecol. 2013;22:620–34.View ArticlePubMedGoogle Scholar
Ohtsu T, Kimura MT, Katagiri C. How drosophila species acquire cold tolerance. Eur J Biochem. 1998;252:608–11.View ArticlePubMedGoogle Scholar
Teets NM, Peyton JT, Ragland GJ, Colinet H, Renault D, Hahn DA, et al. Combined transcriptomic and metabolomic approach uncovers molecular mechanisms of cold tolerance in a temperate flesh fly. Physiol Genomics. 2012;44:764–77.View ArticlePubMedGoogle Scholar
Vesala L, Salminen T, Laiho A, Hoikkala A, Kankare M. Cold tolerance and cold-induced modulation of gene expression in two drosophila virilis group species with different distributions: Cold-induced changes in gene expression. Insect Mol Biol. 2012;21:107–18.View ArticlePubMedGoogle Scholar
Addo-Bediako A, Chown SL, Gaston KJ. Metabolic cold adaptation in insects: A large-scale perspective. Funct Ecol. 2002;16:332–8.View ArticleGoogle Scholar
Kelly MW, Grosberg RK, Sanford E. Trade-offs, geography, and limits to thermal adaptation in a tide pool copepod. Am Nat. 2013;181:846–54.View ArticlePubMedGoogle Scholar
Warren RJ, Chick L. Upward ant distribution shift corresponds with minimum, not maximum, temperature tolerance. Glob Chang Biol. 2013;19:2082–8.View ArticlePubMedGoogle Scholar
Addo-Bediako A, Chown SL, Gaston KJ. Thermal tolerance, climatic variability and latitude. Proc R Soc B Biol Sci. 2000;267:739–45.View ArticleGoogle Scholar
Hodgins-Davis A, Adomas AB, Warringer J, Townsend JP. Abundant gene-by-environment interactions in gene expression reaction norms to copper within saccharomyces cerevisiae. Genome Biol Evol. 2012;4:1061–79.View ArticlePubMedPubMed CentralGoogle Scholar
Aubin-Horth N, Renn SCP. Genomic reaction norms: Using integrative biology to understand molecular mechanisms of phenotypic plasticity. Molecular Ecology. 2009;18:3763–3780.View ArticlePubMedGoogle Scholar
González E, Joly S: Impact of RNA-seq attributes on false positive rates in differential expression analysis of de novo assembled transcriptomes. BMC Research Notes. 2013;6:503.View ArticlePubMedPubMed CentralGoogle Scholar
Sarup P, Sørensen JG, Kristensen TN, Hoffmann AA, Loeschcke V, Paige KN, et al. Candidate genes detected in transcriptome studies are strongly dependent on genetic background. PLoS ONE. 2011;6:e15644.View ArticlePubMedPubMed CentralGoogle Scholar
Pavlidis P, Jensen JD, Stephan W, Stamatakis A. A critical assessment of storytelling: Gene ontology categories and the importance of validating genomic scans. Mol Biol Evol. 2012;29:3237–48.View ArticlePubMedGoogle Scholar
Khaitovich P, Weiss G, Lachmann M, Hellmann I, Enard W, Muetzel B, et al. A neutral model of transcriptome evolution. PLoS Biol. 2004;2:E132.View ArticlePubMedPubMed CentralGoogle Scholar
Ogasawara O, Okubo K. On theoretical models of gene expression evolution with random genetic drift and natural selection. PLoS ONE. 2009;4:e7943.View ArticlePubMedPubMed CentralGoogle Scholar
Gadau J, Helmkampf M, Nygaard S, Roux J, Simola DF, Smith CR, et al. The genomic impact of 100 million years of social evolution in seven ant species. Trends Genet. 2012;28:14–21.View ArticlePubMedPubMed CentralGoogle Scholar
Feder ME, Walser J-C. The biological limitations of transcriptomics in elucidating stress and stress responses. J Evol Biol. 2005;18:901–10.View ArticlePubMedGoogle Scholar
Hoekstra LA, Montooth KL. Inducing extra copies of the hsp70 gene in drosophila melanogaster increases energetic demand. BMC Evol Biol. 2013;13:1–11.View ArticleGoogle Scholar
Arnold FH, Wintrode PL, Miyazaki K, Gershenson A. How enzymes adapt: Lessons from directed evolution. Trends Biochem Sci. 2001;26:100–6.View ArticlePubMedGoogle Scholar
Lubertazzi D. The biology and natural history of aphaenogaster rudis. Psyche (Camb Mass). 2012;2012:1–11.View ArticleGoogle Scholar
Gotelli NJ, Ellison AM. A Primer of Ecological Statistics. 2nd ed. Sunderland: Sinauer Associates, Inc; 2012.Google Scholar
Cottingham KL, Lennon JT, Brown BL. Knowing when to draw the line: Designing more informative ecological experiments. Front Ecol Environ. 2005;3:145–52.View ArticleGoogle Scholar
Lohse M, Bolger AM, Nagel A, Fernie AR, Lunn JE, Stitt M, et al. RobiNA: A user-friendly, integrated software solution for RNA-seq-based transcriptomics. Nucleic Acids Res. 2012;40:W622–7.View ArticlePubMedPubMed CentralGoogle Scholar
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.View ArticlePubMedPubMed CentralGoogle Scholar
Schmieder R, Edwards R. Fast identification and removal of sequence contamination from genomic and metagenomic datasets. PLoS ONE. 2011;6:e17288.View ArticlePubMedPubMed CentralGoogle Scholar
Chen T-W, Gan R-CR WTH, Huang P-J, Lee C-Y, Chen Y-YM, Chen C-C, et al. FastAnnotator- an efficient transcript annotation web tool. BMC Genomics. 2012;13 Suppl 7:S9.View ArticleGoogle Scholar
Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012;131:281–5.View ArticlePubMedGoogle Scholar
Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.Google Scholar
Hothorn T, Hornik K, van de Wiel M, Zeileis A. Implementing a class of permutation tests: The coin package. J Stat Softw. 2008;28:1–23.View ArticleGoogle Scholar
R Core Team. R: A language and environment for statistical computing. 2013.Google Scholar
Grossmann S, Bauer S, Robinson PN, Vingron M. Improved detection of overrepresentation of gene-ontology annotations with parent–child analysis. Bioinformatics. 2007;23:3024–31.View ArticlePubMedGoogle Scholar
Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22:1600–7.View ArticlePubMedGoogle Scholar
