Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome

Abstract

Cold tolerance is a key determinant of insect distribution and abundance, and thermal acclimation can strongly influence organismal stress tolerance phenotypes, particularly in small ectotherms like Drosophila. However, there is limited understanding of the molecular and biochemical mechanisms that confer such impressive plasticity. Here, we use high-throughput mRNA sequencing (RNA-seq) and liquid chromatography - mass spectrometry (LC-MS) to compare the transcriptomes and metabolomes of D. melanogaster acclimated as adults to warm (rearing) (21.5 °C) or cold conditions (6 °C). Cold acclimation improved cold tolerance and led to extensive biological reorganization: almost one third of the transcriptome and nearly half of the metabolome were differentially regulated. There was overlap in the metabolic pathways identified via transcriptomics and metabolomics, with proline and glutathione metabolism being the most strongly-supported metabolic pathways associated with increased cold tolerance. We discuss several new targets in the study of insect cold tolerance (e.g. dopamine signaling and Na⁺ -driven transport), but many previously identified candidate genes and pathways (e.g. heat shock proteins, Ca²⁺ signaling, and ROS detoxification) were also identified in the present study, and our results are thus consistent with and extend the current understanding of the mechanisms of insect chilling tolerance.

Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome

OPEN

R

A P

Heath A. MacMillan,, Jose M. Knee,, Alice B. Dennis,,,*, Hiroko Udaka,, Katie E. Marshall,, Thomas J. S. Merritt & Brent J. Sinclair

Cold tolerance is a key determinant of insect distribution and abundance, and thermal acclimation Drosophila. However, there is limited understanding of the molecular and biochemical mechanisms that confer such impressive plasticity. Here, we use high-throughput mRNA sequencing (RNA-seq) and liquid chromatography mass spectrometry (LC-MS) to compare the transcriptomes and metabolomes of D. melanogaster acclimation improved cold tolerance and led to extensive biological reorganization: almost one third of metabolism being the most strongly-supported metabolic pathways associated with increased cold tolerance. We discuss several new targets in the study of insect cold tolerance (e.g. dopamine signaling and Na+ shock proteins, Ca+ our results are thus consistent with and extend the current understanding of the mechanisms of insect chilling tolerance.

Low temperature tolerance is a key determinant of insect distribution, because the physiological effects of temperature ultimately determine physical performance and reproductive success14. The extensive inter- and intra-specic variation in Drosophila thermal tolerance makes the genus a useful model for linking the eects of temperature on biochemistry and physiology to ecological patterns and processes58. Because of its cosmopolitan distribution in the wild, large degree of thermal tolerance plasticity810, and status as a model organism1113,D. melanogaster Meigen is at the forefront of these investigations. Like many insects, Drosophila enter a state of paralysis termed chill coma at their critical thermal minimum (CTmin)1417. As chill-susceptible insects, adult ies lose ion and water balance during cold exposure, which eventually leads to cell death, tissue damage, physical performance decits, and mortality1719. However, if substantial injury has not yet occurred, ies recover the ability to move when returned to favorable temperatures20. The response to cold in D. melanogaster is phenotypically plastic6,9, labile to articial selection21,22, and segregates genetically11. Drosophila populations from higher latitudes are generally more cold-tolerant (dened as improvement in any cold tolerance trait)23,24,26,27, and there is a

wealth of among-species variation in thermal tolerance within the Drosophila genus that persists when species are reared under common thermal conditions5,25. These evolved dierences in cold tolerance appear, like phenotypic

Department of Biology, University of Western Ontario, London, ON, Canada. Department of Chemistry and Biochemistry, Laurentian University, Sudbury, ON, Canada. Landcare Research, Auckland, New Zealand. Allan Wilson Centre for Molecular Ecology and Evolution, Auckland, New Zealand. Present address: Department of Biology, York University, Toronto, Canada. Present address: Advanced Medical Research Institute of Canada, Sudbury, Canada. Present address: Institute of Integrative Biology, ETH Zrich, Zrich, Switzerland. *Present address: EAWAG, Swiss Federal Institute of Aquatic Science and Technology, Dbendorf, Switzerland. Present address: Department of Zoology, Kyoto University, Kyoto, Japan. Present address: Department of Zoology, University of British Colombia, Vancouver, Canada. Correspondence and requests for materials should be addressed to H.A.M. (email: [email protected])

SCIENTIFIC REPORTS

1

www.nature.com/scientificreports/

plasticity, to be partially driven by modications in hemolymph ion and metabolite composition5,18,28, adaptive

shis in membrane composition that maintain membrane uidity at low temperatures29,30, and an ability to main

tain metabolic homeostasis in the face of chilling stress31.

Surprisingly, in spite of this wealth of phenotypic information, the molecular mechanisms underlying cold tolerance of D. melanogaster are not fully understood, and relatively few candidate genes or pathways have been linked directly to chilling tolerance. Quantitative trait loci (QTL) analyses of chill coma recovery time (CCRT) have only identified three QTLs associated with faster CCRT12,32,33, while microarray and RNA-seq studies have identied few genes consistently associated with cold tolerance or upregulated following cold exposure in D. melanogaster34,35. Although adult acclimation induces such large shis in thermal tolerance phenotypes, one proteomic study identied a handful of proteins (largely muscle- and reproduction-related) that diered from controls (25C) aer ve days acclimation at 11C36. Making direct biochemical links between these existing candidates and cold tolerance phenotypes has also proven difficult. For example, Frost appears to encode a disordered protein37, is consistently upregulated aer cold exposure34,3741, is present in a QTL associated with chill coma32, and shows geographic variation in amino acid sequence consistent with environmental variation42. However, RNAi knockdown of Frost expression does not consistently impair cold tolerance38,43. Thus, although there is substantial cold tolerance plasticity in D. melanogaster, this plasticity is associated with strikingly few candidate molecules.

It is possible that this mismatch between phenotype and candidate genes has arisen because most studies have focused on the response to a single brief sub-lethal cold exposure (e.g.34,35,41) or compared the basal transcriptome among ies with detectable, but oen small, phenotypic dierences (e.g.12). For example, several gene-expression studies have attempted to elucidate the mechanisms underlying the rapid cold-hardening response, which is now thought to be mediated by non-transcriptomic cell signaling44. By contrast, longer-term cold acclimation over a period of days or weeks dramatically improves cold tolerance in almost all insects, including Drosophila18,26,45, and

this robust response occurs over a timescale that is consistent with modications to gene expression and (consequently) metabolic pathways. We postulated that comprehensively identifying changes in the transcriptome and metabolome associated with these large acclimation-induced changes in cold tolerance phenotype would allow us to identify a well-supported suite of candidate genes and pathways, allowing us to strengthen or question existing hypotheses, and generate new hypotheses of the mechanisms underlying variation in cold tolerance.

Here, we compare D. melanogaster adults acclimated for six days at 21.5 and 6C, which leads to a 2.5C shi in CTmin in this population (from 3.4 0.2 to 0.9 0.1 C5). We compare the transcriptomic and metabolomic eects of acclimation temperature, and show that almost one third of the transcriptome and nearly half of the metabolome are dierentially regulated with thermal acclimation. Importantly, we found considerable overlap in the metabolic pathways identied via transcriptomics and metabolomics, which suggests that our candidate genes, metabolites, and metabolic pathways for chilling tolerance are particularly robust. In particular, our data suggest strong roles for glutathione and arginine and proline metabolism in cold acclimation, as these pathways were enriched in both the transcriptome and metabolome. A majority of our candidate genes and pathways are also consistent with our current understanding of the physiology underlying insect chilling tolerance. For example, immune-responsive genes and genes coding for heat shock proteins and Ca2+ binding proteins were dierentially expressed, and thus our results both support and extend the hypothesized mechanisms of insect chilling tolerance.

Results and Discussion

Thermal acclimation shifts the D. melanogaster transcriptome. We examined the impact of thermal acclimation on the transcriptome of D. melanogaster using high-throughput mRNA sequencing (RNA-seq). An average of 11.4 and 10.5 million high-quality 50bp reads (109M reads in total) were mapped to the D. melanogaster genome from each of ve biological replicates of cold- (6C) and warm-acclimated ies (i.e. control ies maintained at their rearing temperature of 21.5 C), respectively (Dataset 1). A total of 11,900 genes had reads mapped to them in every one of the ten samples. Acclimation treatment heavily inuenced the gene expression prole (Figs1a and 2); a total of 4,362 genes (c. 29% of coding genes) were dierentially expressed between the two acclimation temperatures based on a geometric algorithm and a false discovery rate (FDR)-corrected cuto of 0.05 (see methods for details). To focus our eorts on a smaller subset of (highly dierentially-expressed) genes we further subset this to genes with greater than a two-fold dierence in expression and Q<0.01 (P-value following FDR), resulting in 1,567 genes that diered in expression between the two acclimation treatments (c. 10% of the genome; hereaer referred to as dierentially expressed genes), with 649 genes upregulated and 918 downregulated following cold acclimation (Fig.1b, Dataset 2). These numbers are orders of magnitude greater than the number of cold-tolerance candidate genes previously identied following cold exposure in D. melanogaster34,35, and in most other insect species4649.

The majority (985; 63%) of dierentially-expressed genes mapped to at least one gene ontology (GO) term, and 108 GO terms were signicantly enriched in the dataset following FDR correction (Dataset 3). To simplify interpretation of this list, we used REVIGO to cluster GO terms across all dierentially-expressed genes, which resulted in 49 semantically-distinct GO terms that were signicantly enriched in response to acclimation (Dataset 3). We then separated genes by their direction of dierential expression between the two acclimation temperatures and ran the same GO term analysis with each subset, which resulted in 31 and 47 semantically-distinct GO terms with greater expression in cold-acclimated ies and ies maintained at 21.5 C, respectively (Fig.3). Complete lists of signicant GO terms with dierentially expressed genes associated with each term can also be found in the supplementary material (Dataset 3). Several terms for cellular components and processes identied here were also enriched in the previously-sequenced transcriptomes of either D. montana or D. virilis (both in dierent subgenera to D. melanogaster) following cold acclimation (for example, cytoskeletal components, ion transporters, calcium regulation and glutathione-s-transferase activity)48. Despite using similar acclimation conditions

SCIENTIFIC REPORTS

2

www.nature.com/scientificreports/

Figure 1. Overview of the eects of thermal acclimation on gene expression (a,b) and metabolite proles (c,d) of D. melanogaster. (a) Dendrogram of gene expression similarity among biological replicates basedon Jensen-Shannon distance of log10-transformed fragments per kilobase of exon per million mapped reads (FPKM). (b) Biplot of gene expression. (c) Metabolic proles of warm and cold acclimated ies visualizedwith principal component analysis (points represent scores of biological replicates). (d) Biplot of metabolite abundance in cold- and warm acclimated ies. Dierentially expressed genes (b); with 2-fold expression dierence between acclimation treatments and a Q-value< 0.01) or signicantly more abundant metabolites, in warm- and cold-acclimated ies are shown in red or blue, respectively. Genes or metabolites that did not dier in expression or abundance are shown in grey. Full lists of dierentially-expressed and abundant genes and metabolites are presented in Dataset 3 and Table S1, respectively.

to our study, Parker et al.48 associated far fewer candidate genes and GO terms with cold acclimation in these drosophilid species. This dierence may arise from analytical dierences between our approaches, or because the more cold-tolerant D. montana and D. virilis respond less to acclimation at this temperature. Regardless, although the GO terms associated with cold acclimation diered between D. montana and D. virilis, we observe overlap between each of these species and D. melanogaster, suggesting that candidate mechanisms for cold acclimation are generally conserved across the Drosophila phylogeny.

Several of the GO terms that were signicantly enriched in cold-acclimated ies related to the muscle contractile apparatus (Fig.3, Dataset 3), consistent with a previous proteomic survey36. Low temperatures strongly suppress maximal muscle force production in both sh and insects50,51, an eect that can be partially compensated for through acclimation51. Acclimation to low temperatures by poikilothermic sh involves reorganization

SCIENTIFIC REPORTS

3

www.nature.com/scientificreports/

Figure 2. Heat map of mRNA expression proles of male D. melanogaster acclimated to 6C (le) and 21.5C (right) as adults. Columns represent biological replicates (n= 5 per treatment) and rows represent individual genes. Abundance of mRNA is shown as log2-transformed fragments per kilobase of exon per million mapped reads (FPKM) scaled by row (gene). Transcripts with greater expression in a sample appear red and transcripts that are less abundant in a sample appear blue. Gene clustering is based on the similarity (Euclidean distance) of expression proles across samples. Dendrograms denote the overall similarity of gene expression proles by gene (y axis, le; limited to displaying nine large clusters shown as grey boxes), and by sample (x-axis, top). Only genes for which contigs were mapped in all ten samples are shown (n = 11,900 genes). Note that here the similarity among samples is based only on the data presented in the dendrogram, which thus diers from that presented in Fig.1.

of muscle ber type and abundance, a change that is associated with dierential expression of several proteins involved in muscle contraction, including the myosin light and heavy chains52,53. Our results suggest a similar reorganization of the muscle bers during adult cold acclimation in Drosophila, as the genes encoding myosin heavy chain (Mhc) and both myosin light chains 1 and 2 (Mlc1, Mlc2) were upregulated 3.5-fold (Dataset 2). In addition, cold-acclimated ies upregulated genes encoding troponin C and troponin T (including upheld, Tpn25D, Tpn41C, TpnC47D, and TpnC73F; Dataset 3), which facilitate Ca2+-mediated muscle contraction. This increased transcriptional investment in the muscle systems could mitigate the observed depolarization of muscle

SCIENTIFIC REPORTS

4

www.nature.com/scientificreports/

Figure 3. Gene Ontology (GO) terms enriched in genes upregulated with cold acclimation (blue) or downregulated with cold acclimation (i.e. higher expression in warm acclimation; red) in adult D. melanogaster. Bars denote the proportion of dierentially expressed genes (2-fold dierence in expression anda P-value< 0.01) relative to the total number of genes in the D. melanogaster genome mapped to each term. Asterisks denote terms with long titles that have been simplied for space (GO:0016705, GO:0016747). To avoid redundancy, signicant GO term lists were reduced using REVIGO (see methods), but full lists of signicant GO terms, including dierentially expressed genes that map to each term, are presented in Dataset 3).

cell resting potentials54,55 and decreased muscle force generation at low temperatures50. Warm-acclimated D. melanogaster lose the ability to y at relatively mild low temperatures (c. 1418C)56, and our results suggest that ies could modulate muscle force generation in response to low temperature exposure to maintain ight performance even if development at high temperatures has already xed wing morphology56. Calcium regulatory proteins largely determine the thermal sensitivity of the actin-myosin complex, and are also dierentially regulated during cold acclimation in sh57. Thus, these modications to Ca2+ binding and the actin-myosin machinery could be essential to maintenance of muscle contractility at low temperatures, and could be a critical determinant of reproductive success in natural settings. Diapausing mosquitos (Culex pipens) have greater abundance of polymerized actin at muscle ber intersections in the midgut58, and chilling disrupts cytoskeletal organization in primary embryonic cultures of Drosophila cells59. Thus, regulation of cytoskeletal function is likely to be an important component of cold acclimation in a variety of tissues beyond muscle.

Transcripts for several genes that have been previously associated with cold stress or cold tolerance in Drosophila were significantly upregulated in cold-acclimated flies. For example, Frost, smp-30, Starvin and hsr-omega transcripts were all upregulated during cold acclimation in the present study, and have all previously been associated with cold tolerance or cold exposure in D. melanogaster38,6064. Many genes encoding heat shock proteins were upregulated in cold-acclimated ies, including transcripts for HSP22, HSP23, HSP26, HSP67, and HSP70Bbb (Dataset 4). Many of these candidate genes are also upregulated in response to a brief cold stress (e.g. 0C for 2h)34 and/or cold acclimation65 in D. melanogaster. Expression of all of these candidate cold tolerance genes is highly enriched in the Malpighian tubules66, which supports the notion that the renal epithelia, and the Malpighian tubules in particular, are central to insect cold tolerance18,28,67,68, and targets the tubules as a tissue in which to better determine the physiological roles of these molecules. A tissue-focused approach may help to resolve conicting evidence about the role of these candidates in determining cold tolerance phenotypes, as whole-animal knockdowns oen have no eect on cold tolerance phenotypes (e.g.43,65). The upregulation of these genes during acclimation to a mild low temperature that does not induce chill coma supports their importance in acute thermal tolerance, but suggests that their upregulation following a cold stress reects the initiation of thermal acclimation, rather than an adaptive response to acute injury.

Cold acclimation also involved signicant upregulation of genes involved in chitin metabolism and cuticle binding (e.g. knk, verm, and serp), as well as those associated with immunity and lipid transport. Cuticle protein-related genes are also associated with cold tolerance and may be under positive selection in stick insects46.

Further, knockdown of chitin synthesis can disrupt the barrier functions of the midgut peritrophic matrix of our beetles69, and the barrier function of the gut epithelia appear to be of great importance to insect cold tolerance70. Several genes that were upregulated with cold acclimation in the present study have important roles in tracheal

SCIENTIFIC REPORTS

5

www.nature.com/scientificreports/

Figure 4. Heat map of all identiable metabolites in male D. melanogaster acclimated to 6C (le)and 21.5C (right) as adults. Rows represent individual metabolites, and columns represent biological replicates (n= 40 per treatment). Samples and metabolites are grouped according to similarity, denoted by the dendrograms on the x and y-axes.

system development and/or may inuence cuticular permeability (see Dataset 3). A relationship between cold exposure and immunity has also been suggested by numerous transcriptomic studies35,71,72, although the nature of this relationship remains unclear73, particularly since cold acclimation does not appear to improve realized immunity in crickets74. We observed widespread upregulation of both genes encoding antimicrobial peptides (Dro, Dipt, AttA, and AttC) and those associated with immune signaling (Thor, IM3, IM23, Tep4). Cold acclimation also involved upregulation of genes coding for all of the known components of the lipid-transporting larval serum protein complex (Lsp1, Lsp1, Lsp1, Lsp2, and CG8100; Fig.3, Dataset 3), suppression of which is associated with the early stages of diapause induction in the drosophilid Chymomyza costata75,76.

Impacts of thermal acclimation on the Drosophila metabolome. We detected and identied a broad suite of metabolites in ies from both acclimation groups (Fig.4). The overall metabolic proles of the acclimation groups were well-resolved in one axis of a principal component analysis (PC1), along which warm- and cold-acclimated ies were clearly resolved (Fig.1c). A total of 187 features (potential metabolites) were detected, 90 of which (48%) signicantly diered (P< 0.05) between ies acclimated to warm and cold conditions (Fig.1d). We positively identied 34 (18%) features and conrmed them with standards; 26 of these known metabolites signicantly diered in abundance between warm- and cold-acclimated ies (Table S1). Loadings of these conrmed metabolites along PC1 are presented in Fig. S1, and boxplots of individual metabolites are presented in Fig. S2.

Of the total pool of metabolite features, 22% (41/187) were only detected in one of the two groups of ies. We identied two of these features, which were both amino acids (valine and threonine), and both were only above detectable levels in cold-acclimated ies. Several other amino acids were detected and identied in both warm- and cold-acclimated ies and were signicantly more abundant in cold-acclimated ies, including argi-nine, asparagine, glutamine and proline (Figs1d and 5a, Fig. S2). Interestingly, several of these amino acids have been previously reported to decrease during rapid cold-hardening in adult D. melanogaster77, but proline and asparagine increase during long-term thermal acclimation of larvae78, as we observed here in adults. By contrast, adult ies reared under uctuating low temperature conditions increased valine (as observed here) but did not alter proline, glutamine, or threonine levels79. This variation in amino acid responses during dierent windows of development and reproductive maturity support the idea that rapid cold-hardening and long-term acclimation operate through dierent physiological mechanisms80, and suggest that developmental and adult acclimation have distinct metabolic outcomes.

SCIENTIFIC REPORTS

6

www.nature.com/scientificreports/

Figure 5. Results of KEGG pathway analysis of thermal acclimation in D. melanogaster. (a) Results of pathway analysis of identied metabolites in warm- and cold-acclimated ies. Metabolites mapped to pathways are represented along a row. Boxes surrounding metabolite names correspond to dierences in the relative abundance between cold and warm acclimated ies (up-regulated in cold in blue and down-regulated in cold in red). Metabolites shown in purple were only detected in cold-acclimated ies. Total: total number of metabolites in pathway, Hits: number of metabolites in pathway that signicantly dier in abundance between warm- and cold-acclimated ies. Inset: Statistical signicance and pathway impact, with numbers corresponding to the pathways listed in panel A. (b) Metabolic pathways in the KEGG database that were overrepresented in genes that were signicantly dierentially expressed between warm- and cold-acclimated ies. Pathways listed with black text were signicantly enriched when all dierentially expressed genes were included in the analysis. Those shown in blue or red were signicantly enriched when only genes with higher expression in cold- or warm- acclimated ies (respectively) were included in the analysis. Numbers shownin the green circles match to the same pathway shown in panel A, meaning these pathways were enriched in both analyses. Pathways shown in bold remained signicant following false discovery rate correction. (c) Venn diagram of enriched metabolic pathways in based on dierentially expressed genes and dierentially abundant metabolites. Numbers correspond to the number of pathways, before FDR correction in the case of the transcriptome, that were enriched in the metabolome (yellow), transcriptome (blue) or both (green).

SCIENTIFIC REPORTS

7

www.nature.com/scientificreports/

Adenosine mono- di- and triphosphate levels were also all higher in cold-acclimated ies (Figs1d and 5a), which suggests that an increase in the adenylate pool may be advantageous in the cold. Exposure to cold stress that causes injury and death has no eect on energy availability in the form of ATP in the fat body of rebugs81,

and slightly increases ATP availability in the cricket muscles and beetle larvae82,83. By contrast, chronic chilling decreases whole-body ATP levels in esh y larvae, a detrimental state that can be rescued by interrupting the cold exposure with pulses of high temperature84. These dierent eects of chilling on ATP availability may be related to whether chilling more strongly aects metabolic pathways of ATP production or consumption. For example, cold acclimation reduces the amount of active Na+/K+-ATPase protein (a major contributor to total ATP consumption rates) in D. melanogaster5, and thus increased ATP availability may be a passive consequence of reduced ATP consumption, rather than an adaptive mechanism of thermal acclimation. Alternatively, cold acclimation may increase ATP in anticipation of acute stress, potentially to provide a larger ATP pool to fuel repair processes or reestablishment of homeostasis. Increased ATP supply could, for example, improve low temperature survival by facilitating increased rates of ATP-dependent ion-transport in epithelia during or following cold. However, the specic functions of increased ATP availability and the size of the total adenylate pool in cold tolerance remains to be thoroughly examined.

Downregulation of genes and reductions in metabolite abundance with cold acclimation.

Although stress tolerance studies tend to focus on increases in transcript and metabolite abundance, transcriptional downregulation of key processes may be of equal importance in setting thermal tolerance phenotypes. We observed a greater number of transcripts and approximately 1.5 times more semantically-distinct GO terms enriched in ies maintained at their rearing temperature than those acclimated to cold (meaning these transcripts and GO terms were inversely related to cold tolerance; Fig.3). Transcripts that were less abundant with cold acclimation (i.e. enriched in warm-acclimated ies) were associated with juvenile hormone (JH) catalysis (driven by changes in juvenile hormone epoxide hydrolases; JHEH) as well as sodium-driven symporter activity, monosaccharide binding and peptidase (particularly trypsins) and hydrolase activity (Fig.3, Dataset 3). Signaling through JH suppresses ovarian development during short day conditions in D. melanogaster females85 (note that in the present study day length was the same in both acclimation treatments), and JH promotes courtship behavior in males (which were sampled for the present study)86. Furthermore, JH stimulates production of antifreeze proteins in larvae of the freeze avoidant beetle (Dendroides canadensis)87. To our knowledge, the role of JH in chilling tolerance of Drosophila has not been studied, but Jheh1 upregulation is associated with seasonal changes in the phenotype of D. suzukii88. On a shorter timescale, JH could be a candidate for coordinating physiological and behavioural adjustments in response to abiotic temperature, thus integrating multiple seasonal cues.

Several genes encoding proteins involved in multiple classes of sodium driven transport or exchange were also downregulated with cold acclimation (Fig.3, Dataset 3), including those implicated in sodium:amino acid symport (e.g. NAAT1, dmGlut) and sodium:inorganic phosphate symport (e.g. Picot, MFS10, MFS14). If these changes in gene expression reect a reduced reliance on sodium-driven symport mechanisms, such a change may be tied to previously-observed reductions in transcellular sodium gradients in cold-acclimated ies, a shi that is suggested to protect against ion balance disruption at low temperatures18. A transporter (Nha1) recently suggested to function as an H+-Cl cotransporter, rather than a Na+/H+ exchanger as originally suspected89, was also downregulated in cold-acclimated ies. This contrasts with the recent nding that knockdown of Nha1 considerably slows chill coma recovery time of D. melanogaster adults67. Roles for acid-base balance and the thermal sensitivity of hydrogen ion transport in chilling tolerance have received little experimental attention, but could be a fruitful avenue of future research.

A total of ten of the 34 positively-identied metabolites signicantly decreased in abundance following cold acclimation (i.e. were higher in ies maintained at 21.5C) including dopamine, trehalose and histidine (Table S1; Fig.5a). The role of dopamine signaling in thermal preference in D. melanogaster has been examined in two studies using knock-out mutants, but these led to opposite conclusions about whether increases in dopamine signaling increase cold preference or warm preference90,91. Another potential explanation for the decrease we observed in dopamine concentration is that dopamine secreted by epidermal cells is oxidized into melanin92,

and given that Drosophila invest in melanin production during development at low temperatures93, it may be that dopamine is being invested into melanin synthesis during cold acclimation94. Trehalose was also less abundant in cold-acclimated than warm-acclimated ies. This result is curious, given that trehalose is a well-described cryoand osmoprotectant, and increases in abundance in response to cold shock or rapid cold-hardening treatments in Drosophila adults77, as well as during long-term cold acclimation of Drosophila larvae78. In the present study, cold acclimation also reduced levels of histidine in D. melanogaster (Fig.5a), but selection for fast recovery from chill coma has been previously noted to increase histidine levels in the same species31.

Transcriptomic and metabolomic evidence for metabolic pathway regulation. We used MetaboAnalyst95 and goseq96 to nd metabolic pathways in the KEGG database that were signicantly altered in the metabolome and transcriptome, respectively. A total of 28 KEGG pathways were enriched with the metabolites in our analysis, and 25 of these pathways were signicantly modied during thermal acclimation (Fig.5a). A total of 283 of the 1,567 dierentially expressed genes (18%) from the transcriptomic analysis were represented in the KEGG database. Of these, 103 genes were more highly expressed in cold-acclimated ies, and 180 were more highly expressed in ies maintained at 21.5C. We mapped all dierentially-expressed genes, as well as just those enriched in one or the other acclimation group, to KEGG pathways, and present results from each of these three independent analyses in Dataset 4. To maximize identication of overlapping pathways (those enriched in both the transcriptomic and metabolomic data), we pooled signicant pathways from all three transcriptomic analyses (prior to FDR correction).

SCIENTIFIC REPORTS

8

www.nature.com/scientificreports/

In total, 25 metabolic pathways from the KEGG database were signicantly altered during thermal acclimation in each of the two omics approaches. From these, ten enriched pathways were shared across the two approaches (Fig.5b,c), and two of these ten pathways have the strongest support, given that they remained significant following false discovery rate correction of both datasets. These two pathways were 1) Arginine and Proline Metabolism and 2) Glutathione Metabolism. The observation that Arginine and Proline Metabolism genes and metabolites were so strongly enriched in both datasets ts with our current understanding of insect chilling tolerance, as proline has been repeatedly suggested to be preferentially used by the mitochondria as a metabolic fuel at low temperatures, and/or be involved in osmotic protection (at relatively low concentrations) and cryoprotection (at high concentrations) of insect tissues at low temperatures9799. Indeed, a proline-supplemented diet under winter-acclimating conditions confers partial freeze tolerance to D. melanogaster larvae98. In the present study, cold-acclimated ies had higher levels of both proline (1.3-fold) and arginine (1.8-fold) than warm-acclimated ies, and also contained a greater proportion of glutamate (1.4-fold), a precursor to arginine and proline (Table S1). Further, 21 of the 54 genes that map to this pathway (39%) were dierentially expressed during thermal acclimation (Fig.5b). Together, these ndings suggest that these polar and charged amino acids, and modulation of proline metabolism in particular, may be central to the acquisition of chilling tolerance. Curiously, however, selection for rapid recovery from chill coma downregulates proline metabolism, resulting in lower, rather than higher, levels in selected ies31. This discrepancy may arise because proline has dierent roles in chilling survival versus chill-coma recovery, or because selection for cold tolerance elicits trade-os with neighboring pathways such as those related to redox balance. Increases in intracellular or extracellular proline protect mammalian cells against oxidative stress, through preservation of the glutathione redox environment, meaning proline biosynthesis may instead be indicative of an oxidative stress response100.

Glutamate is not only a metabolic precursor of arginine and proline, but also of glutathione, and Glutathione Metabolism was the second KEGG-derived pathway that had strong support for dierential regulation in both the transcriptome and the metabolome. Interestingly, we found no dierence in the abundance of reduced glutathione (GSH), but rather a higher abundance of oxidized glutathione (GSSG) in cold-acclimated ies, relative to ies maintained at their rearing temperature (Fig.5a, Table S1). This apparent shi in redox state could imply that cold-acclimated ies are under greater oxidative stress at 6 C, as cold stress and chilling injury have been associated with oxidative damage and rightward shi in the GSH:GSSG ratio, particularly in chill susceptible and freeze avoiding insect species101103. In the transcriptome of D. melanogaster, we found 17/39 (29%) of genes mapping to the Glutathione Metabolism term were dierentially expressed, including eight genes that encode glutathione-S-transferases (GSTs) and one GST-like protein coding gene (Mgstl), all of which were upregulated with cold acclimation (Dataset 4). Cold acclimation induces increases in GST abundance and activity in the fat body and plasma of larval Lepidoptera104,105, and GSH:GSSG balance plays important roles in cold stress and cold acclimation signaling in plants106, suggesting that this observed shi might serve a signaling, as well as protective, role.

How does cold acclimation work? Cold acclimation is thought to improve insect cold tolerance through an improved ability to maintain homeostasis and cellular integrity at low temperatures107. The primary mechanisms proposed to underlie this phenotypic plasticity are 1) protection against cold-induced ion and water balance disruption18, 2) defense against oxidative stress102, 3) maintenance of metabolic balance31, 4) refolding or stabilizing proteins45, 5) modications to membrane uidity to maintain cellular integrity108, and 6) inhibition of apoptosis80.

Our results support most, but not all, hypothesized mechanisms of chill tolerance. While the whole-animal metabolome cannot oer substantial insight into adjustments to ion and water balance, a large number of genes involved in ion homeostasis were dierentially-expressed with cold acclimation. Our transcriptome and metabolome data suggest that glutathione metabolism may be important in avoiding or repairing oxidative damage at low temperatures, and the generally broad suite of changes we observed in metabolic enzyme transcripts and the metabolome with cold acclimation also supports a role for modications to protein, lipid, and carbohydrate metabolism networks in cold tolerance. The general upregulation of chaperone proteins we observed following cold acclimation supports a role for preemptive protection against misfolded proteins at stressful low temperatures. Although our transcriptomic data indicate altered lipid metabolism during thermal acclimation, we saw little evidence to specically suggest changes to membrane lipid metabolism specically. Because our metabolomics data was biased toward polar metabolites, we were unable to identify membrane components to address membrane-related mechanisms directly. Finally, we did not see strong evidence in the transcriptome (perhaps excepting shis in Ca2+ handling) to suggest changes to apoptotic signaling pathways during long-term adult thermal acclimation.

Our results beg one to ask why many of the molecular and biochemical changes noted here have not been previously associated with the low temperature biology of Drosophila. Previous studies in Drosophila have focused on acute responses to temperature (i.e. during or following minutes to hours in chill coma) rather than long-term acclimation to non- or less-stressful low temperatures. Thus, we hypothesize that these studies were exploring phenotypes driven by signaling or post-translational shis that would not be readily apparent in the transcriptome, and whose immediate signature in the metabolome may also be weak. Future studies should thus be directed at better understanding modications to thermal tolerance occurring over these longer time scales. As roles of specic epithelia and organs (such as the Malpighian tubules and gut epithelia), organ systems (such as the muscles and nerves), or body compartments (such as the hemocoel) in setting thermal tolerance pheno-types are now becoming apparent (e.g.7,28,54,67), we suggest omics approaches can, and should, now be used in a more hypothesis-driven and tissue-specic manner than they have in the past (and in the present study). Such an approach could facilitate rapid discovery of molecular and biochemical plasticity that can be causally, rather than coincidentally, linked to tissue-specic and whole-organism thermal tolerance phenotypes.

SCIENTIFIC REPORTS

9

www.nature.com/scientificreports/

Conclusions

Temperature aects all biochemical reactions in an organism, and thus low temperature acclimation undoubtedly involves a large suite of molecular, biochemical and physiological adjustments across interacting and co-dependent organs and tissues. It is this sub-organismal tuning of function ultimately determines whole organism survival and tness. Here, we generated a dataset that nally reects this expected complexity; we have increased the number of candidate genes for cold tolerance to nearly a third of the genome of D. melanogaster, and note that c. 50% of the measured metabolome is altered during the same acclimation treatment. By combining these two approaches, we have 1) provided strong support for some previously reported mechanisms of cold acclimation (e.g. investment in heat shock proteins and proline synthesis); 2) drawn attention to other mechanisms that are suggested but understudied in insects (e.g. protection against oxidative stress, Ca2+ dynamics, and plasticity of the actin and myosin machinery); 3) contributed to an ongoing discussion of consistently observed but mechanistically unclear molecular responses to cold stress (e.g. Frost, immune response pathways); and 4) provided new pathways of study in the mechanisms of thermal acclimation (e.g. sodium-driven symport and endocrine signaling).

Materials and Methods

Animal husbandry. The Drosophila melanogaster population used in this study was derived from 35 isofe-male lines caught in London and Niagara on the Lake, Ontario, Canada in 2007109 and was maintained at a constant 21.5 0.5C and at 50 5% relative humidity with a 13:11h (L:D) light cycle. Flies were raised from egg to adult in 35mL vials containing ~10mL of banana-based food medium. Adult ies were transferred to 3.7L plastic containers each generation for mass breeding and egg collection during which adult ies were provided with a petri dish of food topped with yeast paste (1.5:1 active yeast:water). Aer 24h, the food plate was replaced with a fresh plate of food, which was le for 14h before eggs were collected into fresh vials at a density of 75 eggs per vial.

Newly-eclosed adult males were collected under light CO2 anesthesia on the day of their emergence and divided randomly into two thermal acclimation groups with the same light and humidity conditions: 21.5 C (warm-acclimated) and 6 0.5 C (cold-acclimated). Flies were le at their respective acclimation temperatures for ve days to acclimate and to control for age and anaesthesia eects on cold tolerance63,110, before being snap-frozen in liquid nitrogen and held at 80C for later analysis of gene expression and metabolite abundance.

Gene expression. Transcript abundance was quantified using high-throughput mRNA sequencing (RNA-seq). Frozen warm- and cold-acclimated male D. melanogaster (ve biological replicates per acclimation temperature, each containing 25 ies) were homogenized by mortar and pestle over liquid nitrogen and then in 1 mL of TRIzol reagent. Homogenized samples were centrifuged at 12,000 g for 5 min at 4 C and the bottom 200L (lipids and DNA) discarded. Chloroform (200L) was added to the tubes, samples were shaken for 15s, and le at 22 C for 3 min. Samples were centrifuged at 12,000 g for 15 min at 4 C, the upper (aqueous) phase transferred to a new tube with 1 volume of 70% ethanol, and vortexed (30s). For further sample purication, a 700 L aliquot of the sample was centrifuged (8000 g for 15 s at 22 C) through an RNeasy Mini spin column (Qiagen, Hilden, Germany) in a 2mL collection tube. The effluent was discarded and the sample (bound to the column) was washed with 700L of RW1 buer (Qiagen), and 2500L of RPE buer (Qiagen) using the same centrifuge settings. The column was placed into a new collection tube and RNase-free water (50 L) was centrifuged through the column (8000g for 1min at 22C) twice to release the RNA from the column. Total mRNA was puried from RNA, and cDNA libraries were prepared using a TruSeq non-stranded mRNA preparation kit (Illumina Inc., San Diego, CA, USA) following the TruSeq RNA sample prep v2 protocol (Illumina). Briey, mRNA was puried by poly-A selection using magnetic beads bound to poly-T oligo-nucleotides and chemically fragmented. Fragmented mRNA was reverse-transcribed into single stranded cDNA by reverse transcriptase and the RNA on the opposite strand was replaced with DNA to yield double stranded (ds) cDNA. Overhangs at the 3 and 5 ends of ds cDNA were cleaved and repaired by exonuclease and DNA polymerase, respectively. Illumina HiSeq2000 was used to sequence multiplexed (1314 samples per lane, run alongside samples from two other experiments) single-end 50bp cDNA libraries at the High Throughput Genomics Core Facility at the Huntsman Cancer Institute, University of Utah.

Transcriptome data manipulation and statistical analyses were completed in Galaxy111 following the analysis workow of Trapnell et al.112. Briey, Illumina adapter sequences were clipped from reads using FastQ_Clipper (part of the FastX Toolkit), and clipped reads were aligned to the genome of D. melanogaster (Ensembl build 5.25) using TopHat v.2.0.8113, with alignment limited to known splice junctions. Cufflinks112 was used to assemble and count transcript reads for each sample and Cumerge v.2.0.2 was used to merge assembled transcripts for all biological replicates into a single reference transcriptome114. To test for dierential expression of genes that mapped to the transcriptome, we used CuDi 2, which uses a geometric-based algorithm to generate test statistics, P-values, and Benjamini-Hochberg115 false discovery rate (FDR) corrected P-values (Q-values) for each gene mapped116. More than 4,300 genes were signicantly dierentially expressed with thermal acclimation aer FDR correction (Q < 0.05), so to reduce the number of targets we set cuto values for the magnitude and statistical signicant of dierential expression. Genes that had both a 2-fold dierence in expression between the acclimation groups, and a Q-value less than 0.01 were used in all downstream analyses and were deemed dierentially expressed genes. Dierentially expressed genes were mapped to the Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases using goseq in R96. To reduce the number of signicant GO terms for plotting purposes, we used REVIGO, which eliminates semantically redundant terms117.

Metabolite abundance. Metabolites were extracted from whole y samples (n = 40 per acclimation temperature) following Knee et al.118, with minor modications as described below. Samples consisted of groups of 15 male ies that were weighed to the nearest 0.01mg and ash frozen in liquid nitrogen. The volume of extraction

SCIENTIFIC REPORTS

10

www.nature.com/scientificreports/

solvent used for each sample was normalized by the mass of the sample by pipetting 6.35L extraction solvent per mg of y. The extraction solvent consisted of a 3:1:1 ratio of methanol:chloroform:water. Samples were homogenized with a mixer mill (Qiagen) at 30 Hz for 1 min. Homogenates were placed at 20 C for 60 minutes to promote metabolite extraction and increase protein precipitation. The samples were subsequently centrifuged at 16000 RCF to remove protein and debris, and the supernatant was stored at 80C for later analysis.

Metabolomic analysis was carried out as previously described118. Briey, LC-MS analysis was performed on a Dionex UltiMate 3000 Rapid Separation (U)HPLC (ThermoFisher) coupled to a micrOTOF QII mass spectrometer equipped with an Electrospray Ionization source (Bruker Daltonics). Samples (10L) were injected onto a Kinetex 100 mm 2.1 mm, 1.7 m C18 column (Phenomenex). Metabolites were separated using a multi-step, two solvent gradient with diamylammonium acetate, pH 4.95, as solvent A and methanol as solvent B and a ow-rate of 0.4mL/min (Table S2). Data was acquired at 2Hz in negative ionization mode.

Raw metabolite data was subject to automatic intra-run calibration with sodium formate and putative metabolites identied using the Find Molecular Features peak picking soware (Bruker Daltonics). Principal Component Analysis (PCA; with Pareto scaling) was performed in R, and and T-tests were performed on ProleAnalysis soware (Bruker). Only features that were detected in at least 50% of samples of both group types were used, with the exception of features that were completely undetected in one group type and detected in at least 50% of the other (which we refer to as cases of extreme regulation). Integrated peak signals were normalized by the sum of the feature values in the analysis.

Pathway Analysis, a combination of quantitative enrichment analysis and topology analysis, was carried out on MetaboAnalyst online soware using characterized KEGG pathways for D. melanogaster as the back-end knowledge95,119. Enrichment analysis is based on concentration values of compounds loaded into the analysis and uses a Global Test120 to determine if pathways are signicantly altered by an experimental condition. Topology analysis gauges the potential impact metabolite concentrations have on a particular pathway based on where the metabolites fall in its structure. Topology in our analysis was determined with a relative-betweeness centrality and normalized to the sum of the calculated importance of the pathway95,119.

References

1. Overgaard, J., Kearney, M. R. & Homann, A. A. Sensitivity to thermal extremes in Australian Drosophila implies similar impacts of climate change on the distribution of widespread and tropical species. Glob. Chang. Biol. 20, 17381750 (2014).

2. Andersen, J. L. et al. How to assess Drosophila cold tolerance: chill coma temperature and lower lethal temperature are the best predictors of cold distribution limits. Funct. Ecol. 29, 5565 (2015).

3. Warren, R. J. & Chick, L. Upward ant distribution shi corresponds with minimum, not maximum, temperature tolerance. Glob. Chang. Biol. 20822088, doi: 10.1111/gcb.12169 (2013).

4. Parmesan, C. & Yohe, G. A globally coherent ngerprint of climate change impacts across natural systems. Nature 421, 3742 (2003).

5. MacMillan, H. A. et al. Parallel ionoregulatory adjustments underlie phenotypic plasticity and evolution of Drosophila cold tolerance. J. Exp. Biol. 218, 423432 (2015).

6. Overgaard, J., Srensen, J. G., Com, E. & Colinet, H. The rapid cold hardening response of Drosophila melanogaster: Complex regulation across dierent levels of biological organization. J. Insect Physiol. 62, 4653 (2014).

7. Armstrong, G. A. B., Rodrguez, E. C. & Robertson, R. M. Cold hardening modulates K+ homeostasis in the brain of Drosophila melanogaster during chill coma. J. Insect Physiol. 58, 15111516 (2012).

8. Gerken, A. R., Eller, O. C., Hahn, D. a. & Morgan, T. J. Constraints, independence, and evolution of thermal plasticity: Probing genetic architecture of long- and short-term thermal acclimation. Proc. Natl. Acad. Sci. 112, 43994404 (2015).

9. Ransberry, V. E., MacMillan, H. A. & Sinclair, B. J. The relationship between chill-coma onset and recovery at the extremes of the thermal window of Drosophila melanogaster. Physiol. Biochem. Zool. 84, 553559 (2011).

10. Kelty, J. Rapid cold-hardening of Drosophila melanogaster in a eld setting. Physiol. Entomol. 32, 343350 (2007).11. Mackay, T. F. C. et al. The Drosophila melanogaster genetic reference panel. Nature 482, 173178 (2012).12. Morgan, T. J. & Mackay, T. F. C. Quantitative trait loci for thermotolerance phenotypes in Drosophila melanogaster. Heredity (Edinb). 96, 23242 (2006).

13. Markow, T. A. & OGrady, P. M. Drosophila biology in the genomic age. Genetics 177, 126976 (2007).14. Mellanby, K. Low temperature and insect activity. Proc. R. Soc. B 127, 473487 (1939).15. Baust, J. G. & Rojas, R. R. Insect cold hardiness: Facts and fancy. J. Insect Physiol. 31, 755759 (1985).16. Hazell, S. P. & Bale, J. S. Low temperature thresholds: are chill coma and CTmin synonymous? J. Insect Physiol. 57, 10851089 (2011).

17. MacMillan, H. A. & Sinclair, B. J. Mechanisms underlying insect chill-coma. J. Insect Physiol. 57, 1220 (2011).18. MacMillan, H. A., Andersen, J. L., Loeschcke, V. & Overgaard, J. Sodium distribution predicts the chill tolerance of Drosophila melanogaster raised in dierent thermal conditions. Am. J. Physiol. Regul. Integr. Comp. Physiol. 308, 823831 (2015).

19. Yi, S.-X., Moore, C. W. & Lee, R. E. J. Rapid cold-hardening protects Drosophila melanogaster from cold-induced apoptosis. Apoptosis 12, 11831193 (2007).

20. David, R. J. et al. Cold stress tolerance in Drosophila: analysis of chill coma recovery in D. melanogaster. J. Therm. Biol. 23, 291299 (1998).

21. Gilchrist, G. W., Huey, R. B. & Partridge, L. Thermal sensitivity of Drosophila melanogaster: evolutionary responses of adults and eggs to laboratory natural selection at dierent temperatures. Physiol. Zool. 70, 40314 (1997).

22. Anderson, A., Hoffmann, A. A. & McKechnie, S. W. Response to selection for rapid chill-coma recovery in Drosophila melanogaster: physiology and life-history traits. Genet. Res. 85, 1522 (2005).

23. Homann, A. A. Physiological climatic limits in Drosophila: patterns and implications. J. Exp. Biol. 213, 870880 (2010).24. Goto, S. G. & Kimura, M. T. Heat- and cold-shock responses and temperature adaptations in subtropical and temperate species of Drosophila. J. Insect Physiol. 44, 12331239 (1998).

25. Kellermann, V. et al. Phylogenetic constraints in key functional traits behind species climate niches: patterns of desiccation and cold resistance across 95 Drosophila species. Evolution 66, 33773389 (2012).

26. Gibert, P. & Huey, R. B. Chill-coma temperature in Drosophila: Eects of developmental temperature, latitude, and phylogeny. Physiol. Biochem. Zool. 74, 429434 (2001).

27. Overgaard, J., Kristensen, T. N., Mitchell, K. A. & Homann, A. A. Thermal tolerance in widespread and tropical Drosophila species: does phenotypic plasticity increase with latitude? Am. Nat. 178, S80S96 (2011).

28. MacMillan, H. A., Andersen, J. L., Davies, S. A. & Overgaard, J. The capacity to maintain ion and water homeostasis underlies interspecic variation in Drosophila cold tolerance. Sci. Rep. 5, 18607 (2015).

SCIENTIFIC REPORTS

11

www.nature.com/scientificreports/

29. Cooper, B. S., Hammad, L. a. & Montooth, K. L. Thermal adaptation of cellular membranes in natural populations of Drosophila melanogaster. Funct. Ecol. 28, 886894 (2014).

30. Slotsbo, S. et al. Tropical to sub-polar gradient in phospholipid composition suggests adaptive tuning of biological membrane function in drosophilids. Funct. Ecol., doi: 10.1111/1365-2435.12568 (2015).

31. Williams, C. M. et al. Cold adaptation shapes the robustness of metabolic networks in Drosophila melanogaster. Evolution 68, 350523 (2014).

32. Norry, F. M., Gomez, F. H. & Loeschcke, V. Knockdown resistance to heat stress and slow recovery from chill coma are genetically associated in a quantitative trait locus region of chromosome 2 in Drosophila melanogaster. Mol. Ecol. 16, 327484 (2007).

33. Norry, F. M., Scannapieco, A. C., Sambucetti, P., Bertoli, C. I. & Loeschcke, V. QTL for the thermotolerance eect of heat hardening, knockdown resistance to heat and chill-coma recovery in an intercontinental set of recombinant inbred lines of Drosophila melanogaster. Mol. Ecol. 17, 457081 (2008).

34. Qin, W., Neal, S. J., Robertson, R. M., Westwood, J. T. & Walker, V. K. Cold hardening and transcriptional change in Drosophila melanogaster. Insect Mol. Biol. 14, 60713 (2005).

35. Zhang, J., Marshall, K. E., Westwood, J. T., Clark, M. S. & Sinclair, B. J. Divergent transcriptomic responses to repeated and single cold exposures in Drosophila melanogaster. J. Exp. Biol. 214, 40219 (2011).

36. Colinet, H., Overgaard, J., Com, E. & Srensen, J. G. Proteomic proling of thermal acclimation in Drosophila melanogaster. Insect Biochem. Mol. Biol. 43, 35265 (2013).

37. Bing, X., Zhang, J. & Sinclair, B. J. A comparison of Frost expression among species and life stages of Drosophila. Insect Mol. Biol. 21, 319 (2012).

38. Colinet, H., Lee, S. F. & Homann, A. Functional characterization of the Frost gene in Drosophila melanogaster: importance for recovery from chill coma. PLos One 5, e10925 (2010).

39. Udaka, H., Ueda, C. & Goto, S. G. Survival rate and expression of Heat-shock protein 70 and Frost genes aer temperature stress in Drosophila melanogaster lines that are selected for recovery time from temperature coma. J. Insect Physiol. 56, 18891894 (2010).

40. Goto, M., Li, Y. & Honma, T. Changes of diapause and cold hardiness in the Shonai ecotype larvae of the rice stem borer, Chilo suppressalis Walker (Lepidoptera: Pyralidae) during overwintering. Appl. Entomol. Zool. 36, 323328 (2001).

41. Sinclair, B. J., Gibbs, A. G. & Roberts, S. P. Gene transcription during exposure to, and recovery from, cold and desiccation stress in Drosophila melanogaster. Insect Mol. Biol. 16, 435443 (2007).

42. Rako, L., Blacket, M. J., McKechnie, S. W. & Homann, A. A. Candidate genes and thermal phenotypes: identifying ecologically important genetic variation for thermotolerance in the Australian Drosophila melanogaster cline. Mol. Ecol. 16, 29482957 (2007).

43. Udaka, H., Percival-Smith, A. & Sinclair, B. J. Increased abundance of Frost mRNA during recovery from cold stress is not essential for cold tolerance in adult Drosophila melanogaster. Insect Mol. Biol. 22, 54150 (2013).

44. Teets, N. M., Yi, S.-X., Lee, R. E. & Denlinger, D. L. Calcium signaling mediates cold sensing in insect tissues. Proc. Natl. Acad. Sci. USA 110, 91549159 (2013).

45. Colinet, H. & Homann, A. A. Comparing phenotypic eects and molecular correlates of developmental, gradual and rapid cold acclimation responses in Drosophila melanogaster. Funct. Ecol. 26, 8493 (2012).

46. Dennis, A. B., Dunning, L. T., Sinclair, B. J. & Buckley, T. R. Parallel molecular routes to cold adaptation in eight genera of New Zealand stick insects. Sci. Rep. 5, 13965 (2015).

47. Dunning, L. T. et al. Identication of cold-responsive genes in a New Zealand alpine stick insect using RNA-Seq. Comp. Biochem. Physiol. D 8, 2431 (2013).

48. Parker, D. J. et al. How consistent are the transcriptome changes associated with cold acclimation in two species of the Drosophila virilis group? Heredity (Edinb). 115, 1321 (2015).

49. Sulmon, C. et al. Abiotic stressors and stress responses: What commonalities appear between species across biological organization levels? Environ. Pollut. 202, 6677 (2015).

50. Findsen, A., Pedersen, T. H., Petersen, A. G., Nielsen, O. B. & Overgaard, J. Why do insects enter and recover from chill coma? Low temperature and high extracellular potassium compromise muscle function in Locusta migratoria. J. Exp. Biol. 217, 12971306 (2014).

51. Langfeld, K. S., Crockford, T. & Johnston, I. A. Temperature acclimation in the common carp: Force-velocity characteristics and myosin subunit composition of slow muscle bres. J. Exp. Biol. 155, 291304 (1991).

52. Scott, G. R. & Johnston, I. A. Temperature during embryonic development has persistent eects on thermal acclimation capacity in zebrash. Proc. Natl. Acad. Sci. USA 109, 1424752 (2012).

53. Gracey, A. Y. et al. Coping with cold: An integrative, multitissue analysis of the transcriptome of a poikilothermic vertebrate. Proc. Natl. Acad. Sci. USA 101, 169705 (2004).

54. Andersen, J. L., MacMillan, H. A. & Overgaard, J. Muscle membrane potential and insect chill coma. J. Exp. Biol. 218, 24922495 (2015).

55. Hosler, J. S., Burns, J. E. & Esch, H. E. Flight muscle resting potential and species-specic dierences in chill-coma. J. Insect Physiol. 46, 621627 (2000).

56. Frazier, M. R., Harrison, J. F., Kirkton, S. D. & Roberts, S. P. Cold rearing improves cold-ight performance in Drosophila via changes in wing morphology. J. Exp. Biol. 211, 21162122 (2008).

57. Johnston, I. A. Calcium regulatory proteins and temperature acclimation of actomyosin ATPase from a eurythermal teleost (Carassius auratus L.). J. Comp. Physiol. 129, 163167 (1979).

58. Kim, M., Robich, R. M., Rinehart, J. P. & Denlinger, D. L. Upregulation of two actin genes and redistribution of actin during diapause and cold stress in the northern house mosquito, Culex pipiens. J. Insect Physiol. 52, 12261233 (2006).

59. Cottam, D. M. et al. Non-centrosomal microtubule-organising centres in cold-treated cultured Drosophila cells. Cell Motil. Cytoskeleton 63, 88100 (2006).

60. Goto, S. Expression of Drosophila homologue of senescence marker protein-30 during cold acclimation. J. Insect Physiol. 46, 11111120 (2000).

61. Clowers, K. J., Lyman, R. F., Mackay, T. F. C. & Morgan, T. J. Genetic variation in senescence marker protein-30 is associated with natural variation in cold tolerance in Drosophila. Genet. Res. (Camb). 92, 103113 (2010).

62. Collinge, J. E., Anderson, A. R., Weeks, A. R., Johnson, T. K. & McKechnie, S. W. Latitudinal and cold-tolerance variation associate with DNA repeat-number variation in the hsr-omega RNA gene of Drosophila melanogaster. Heredity (Edinb). 101, 26070 (2008).

63. Colinet, H., Siaussat, D., Bozzolan, F. & Bowler, K. Rapid decline of cold tolerance at young age is associated with expression of stress genes in Drosophila melanogaster. J. Exp. Biol. 216, 2539 (2013).

64. Colinet, H. & Homann, A. Gene and protein expression of Drosophila Starvin during cold stress and recovery from chill coma. Insect Biochem. Mol. Biol. 40, 425428 (2010).

65. ttina, T., Kotl, V. & Korbelov, J. The role of inducible Hsp70, and other heat shock proteins, in adaptive complex of cold tolerance of the fruit y (Drosophila melanogaster). PLos One 10, e0128976 (2015).

66. Chintapalli, V. R., Wang, J. & Dow, J. A. T. Using FlyAtlas to identify better Drosophila melanogaster models of human disease. Nat. Genet. 39, 715720 (2007).

67. Terhzaz, S. et al. Insect capa neuropeptides impact desiccation and cold tolerance. Proc. Natl. Acad. Sci. 112, 28822887 (2015).68. MacMillan, H. A., Williams, C. M., Staples, J. F. & Sinclair, B. J. Reestablishment of ion homeostasis during chill-coma recovery in the cricket Gryllus pennsylvanicus. Proc. Natl. Acad. Sci. USA 109, 2075020755 (2012).

SCIENTIFIC REPORTS

12

www.nature.com/scientificreports/

69. Kelkenberg, M., Odman-Naresh, J., Muthukrishnan, S. & Merzendorfer, H. Chitin is a necessary component to maintain the barrier function of the peritrophic matrix in the insect midgut. Insect Biochem. Mol. Biol. 56, 2128 (2015).

70. MacMillan, H. A. & Sinclair, B. J. The role of the gut in insect chilling injury: cold-induced disruption of osmoregulation in the fall eld cricket, Gryllus pennsylvanicus. J. Exp. Biol. 214, 726734 (2011).

71. Vermeulen, C. J., Srensen, P., Kirilova Gagalova, K. & Loeschcke, V. Transcriptomic analysis of inbreeding depression in cold-sensitive Drosophila melanogaster shows upregulation of the immune response. J. Evol. Biol. 26, 1890902 (2013).

72. Xu, J. & James, R. R. Temperature stress aects the expression of immune response genes in the alfalfa leafcutting bee, Megachile rotundata. Insect Mol. Biol. 21, 269280 (2012).

73. Sinclair, B. J., Ferguson, L. V., Salehipour-shirazi, G. & MacMillan, H. A. Cross-tolerance and cross-talk in the cold: relating low temperatures to desiccation and immune stress in insects. Integr. Comp. Biol. 53, 545556 (2013).

74. Ferguson, L. V., Heinrichs, D. E. & Sinclair, B. J. Paradoxical acclimation responses in the thermal performance of insect immunity. Oecologia 181, 7785 (2016).

75. Roberts, D. B. & Brock, H. W. The major serum proteins of Dipteran larvae. Experientia 37, 103110 (1981).76. Poupardin, R. et al. Early transcriptional events linked to induction of diapause revealed by RNAseq in larvae of drosophilid y, Chymomyza costata. BMC Genomics 16, 720 (2015).

77. Overgaard, J. et al. Metabolomic proling of rapid cold hardening and cold shock in Drosophila melanogaster. J. Insect Physiol. 53, 12181232 (2007).

78. Kotl, V. et al. Long-term cold acclimation extends survival time at 0C and modies the metabolomic proles of the larvae of the fruit y Drosophila melanogaster. PLos One 6, e25025 (2011).

79. Colinet, H., Larvor, V., Laparie, M. & Renault, D. Exploring the plastic response to cold acclimation through metabolomics. Funct. Ecol. 26, 711722 (2012).

80. Teets, N. M. & Denlinger, D. L. Physiological mechanisms of seasonal and rapid cold-hardening in insects. Physiol. Entomol. 38, 105116 (2013).

81. Kotl, V., Vambera, J. & Bastl, J. On the nature of pre-freeze mortality in insects: water balance, ion homeostasis and energy charge in the adults of Pyrrhocoris apterus. J. Exp. Biol. 207, 15091521 (2004).

82. MacMillan, H. A., Williams, C. M., Staples, J. F. & Sinclair, B. J. Metabolism and energy supply below the critical thermal minimum of a chill-susceptible insect. J. Exp. Biol. 215, 13661372 (2012).

83. Colinet, H. Disruption of ATP homeostasis during chronic cold stress and recovery in the chill susceptible beetle (Alphitobius diaperinus). Comp. Biochem. Physiol. A 160, 637 (2011).

84. Dollo, V. H., Yi, S.-X. & Lee, R. E. High temperature pulses decrease indirect chilling injury and elevate ATP levels in the esh y, Sarcophaga crassipalpis. Cryobiology 60, 351353 (2010).

85. Saunders, D. S., Richard, D. S., Applebaum, S. W., Ma, M. & Gilbert, L. I. Photoperiodic diapause in Drosophila melanogaster involves a block to the juvenile hormone regulation of ovarian maturation. Gen. Comp. Endocrinol. 79, 174184 (1990).

86. Wijesekera, T. P., Saurabh, S. & Dauwalder, B. Juvenile Hormone Is Required in Adult Males for Drosophila Courtship. PLos One 11, e0151912 (2016).

87. Horwath, K. L. & Duman, J. G. Induction of antifreeze protein production by juvenile hormone in larvae of the beetle, Dendroides canadensis. J. Comp. Physiol. B 151, 233240 (1983).

88. Toxopeus, J., Jakobs, R., Ferguson, L. V., Gariepy, T. D. & Sinclair, B. J. Reproductive arrest and stress resistance in winter-acclimated Drosophila suzukii. J. Insect Physiol. 89, 3751 (2016).

89. Chintapalli, V. R. et al. Transport proteins NHA1 and NHA2 are essential for survival, but have distinct transport modalities. Proc. Natl. Acad. Sci. 112, 1172011725 (2015).

90. Bang, S. et al. Dopamine signalling in mushroom bodies regulates temperature-preference behaviour in Drosophila. PLos Genet. 7, (2011).

91. Ueno, T., Tomita, J., Kume, S. & Kume, K. Dopamine modulates metabolic rate and temperature sensitivity in Drosophila melanogaster. PLos One 7, e31513 (2012).

92. Yamamoto, S. & Seto, E. S. Dopamine dynamics and signaling in Drosophila: An overview of genes, drugs and behavioral paradigms. Exp. Anim. 63, 107119 (2014).

93. Crill, W. D., Huey, R. B. & Gilchrist, G. W. Within- and between-generation eects of temperature on the morphology and physiology of Drosophila melanogaster. Evolution 50, 12051218 (1996).

94. Kutch, I. C., Sevgili, H., Wittman, T. & Fedorka, K. M. Thermoregulatory strategy may shape immune investment in Drosophila melanogaster. J. Exp. Biol. 217, 36649 (2014).

95. Xia, J., Mandal, R., Sinelnikov, I. V., Broadhurst, D. & Wishart, D. S. MetaboAnalyst 2.0-a comprehensive server for metabolomic data analysis. Nucleic Acids Res. 40, W127W133 (2012).

96. Young, M. D., Wakeeld, M. J., Smyth, G. K. & Oshlack, A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14 (2010).

97. Misener, S. R., Chen, C.-P. & Walker, V. K. Cold tolerance and proline metabolic gene expression in Drosophila melanogaster.J. Insect Physiol. 47, 393400 (2001).98. Kotl, V., imek, P., Zahradnkov, H., Cimlov, J. & ttina, T. Conversion of the chill susceptible fruit y larva (Drosophila melanogaster) to a freeze tolerant organism. Proc. Natl. Acad. Sci. USA 109, 32703274 (2012).

99. Yancey, P. H. Organic osmolytes as compatible, metabolic and counteracting cytoprotectants in high osmolarity and other stresses. J. Exp. Biol. 208, 281930 (2005).100. Krishnan, N., Dickman, M. B. & Becker, D. F. Proline modulates the intracellular redox environment and protects mammalian cells against oxidative stress. Free Radic. Biol. Med. 44, 67181 (2008).101. Lalouette, L., Williams, C. M., Hervant, F., Sinclair, B. J. & Renault, D. Metabolic rate and oxidative stress in insects exposed to low temperature thermal uctuations. Comp. Biochem. Physiol. A 158, 22934 (2011).102. Rojas, R. R. & Leopold, R. A. Chilling injury in the housey: evidence for the role of oxidative stress between pupariation and emergence. Cryobiology 33, 447458 (1996).103. Joanisse, D. & Storey, K. Oxidative stress and antioxidants in overwintering larvae of cold-hardy goldenrod gall insects. J. Exp. Biol.

199, 148391 (1996).104. Grubor-Lajsic, G. et al. Eect of cold acclimation on the antioxidant defense system of two larval Lepidoptera (Noctuidae). Arch.

Insect Biochem. Physiol. 36, 110 (1997).105. Shang, Q. et al. Proteomics analysis of overexpressed plasma proteins in response to cold acclimation in Ostrinia furnacalis. Arch.

Insect Biochem. Physiol. 90, 195208 (2015).106. Kocsy, G., Galiba, G. & Brunold, C. Role of glutathione in adaptation and signalling during chilling and cold acclimation in plants.

Physiol. Plant. 113, 158164 (2001).107. Hayward, S. A. L., Manso, B. & Cossins, A. R. Molecular basis of chill resistance adaptations in poikilothermic animals. J. Exp. Biol.

217, 615 (2014).108. Overgaard, J. et al. Eects of acclimation temperature on thermal tolerance and membrane phospholipid composition in the fruit y Drosophila melanogaster. J. Insect Physiol. 54, 619629 (2008).109. Marshall, K. E. & Sinclair, B. J. Repeated stress exposure results in a survival-reproduction trade-o in Drosophila melanogaster.

Proc. R. Soc. B 277, 963969 (2010).

SCIENTIFIC REPORTS

13

www.nature.com/scientificreports/

110. Nilson, T. L., Sinclair, B. J. & Roberts, S. P. The eects of carbon dioxide anesthesia and anoxia on rapid cold-hardening and chill coma recovery in Drosophila melanogaster. J. Insect Physiol. 52, 10271033 (2006).

111. Goecks, J., Nekrutenko, A., Taylor, J. & Team, T. G. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol. 11, R86 (2010).

112. Trapnell, C. et al. Dierential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat.

Protoc. 7, 562578 (2012).113. Trapnell, C., Pachter, L. & Salzberg, S. L. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 11051111 (2009). 114. Trapnell, C. et al. Transcript assembly and quantication by RNA-Seq reveals unannotated transcripts and isoform switching during cell dierentiation. Nat. Biotechnol. 28, 511515 (2010).115. Benjamini, Y. & Hochberg, Y. Controlling the false discuvery rate: a practical and powerful approach to multiple testing. J. R. Stat.

Soc. Ser. B 57, 289300 (1995).116. Trapnell, C. et al. Dierential analysis of gene regulation at transcript resolution with RNA-seq. Nat. Biotechnol. 31, 4653 (2013). 117. Supek, F., Bonjak, M., kunca, N. & muc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLos One 6, e21800 (2011).118. Knee, J. M., Rzezniczak, T. Z., Barsch, A., Guo, K. Z. & Merritt, T. J. S. A novel ion pairing LC/MS metabolomics protocol for study of a variety of biologically relevant polar metabolites. J. Chromatogr. B. Analyt. Technol. Biomed. Life Sci. 936, 6373 (2013).119. Xia, J., Psychogios, N., Young, N. & Wishart, D. S. MetaboAnalyst: a web server for metabolomic data analysis and interpretation.

Nucleic Acids Res. 37, W652W660 (2009).120. Goeman, J. J. & Bhlmann, P. Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23,

980987 (2007).

Acknowledgements

We thank Caroline Williams and Johannes Overgaard for useful discussion on this work, Brian Dalley at the University of Utah for Illumina sequencing and support, and two anonymous referees for comments that improved the ms. This project was supported by Discovery Grants to BJS and TJSM and postgraduate scholarships to HAM and KEM from the Natural Sciences and Engineering Research Council of Canada, as well as a Canada Research Chair to TJSM and a Banting Postdoctoral Fellowship to HAM.

Author Contributions

H.A.M. and B.J.S. designed the experiment. H.A.M., J.M.K., A.B.D., K.E.M. and H.U. conducted the experiments and collected the data. H.A.M. and J.M.K. analyzed the data, H.A.M., J.M.K., B.J.S. and T.J.S.M. draed the manuscript, and all authors revised the manuscript.

Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing nancial interests: The authors declare no competing nancial interests.

How to cite this article: MacMillan, H. A. et al. Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome. Sci. Rep. 6, 28999; doi: 10.1038/srep28999 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

SCIENTIFIC REPORTS

14