Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome

1Low temperature tolerance is a key determinant of insect distribution, because the physiological effects of temperature ultimately determine physical performance and reproductive success. The extensive inter- and intra-specific variation in Drosophila thermal tolerance makes the genus a useful model for linking the effects of temperature on biochemistry and physiology to ecological patterns and processes. Because of its cosmopolitan distribution in the wild, large degree of thermal tolerance plasticity, and status as a model organism, 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). As chill-susceptible insects, adult flies lose ion and water balance during cold exposure, which eventually leads to cell death, tissue damage, physical performance deficits, and mortality. However, if substantial injury has not yet occurred, flies recover the ability to move when returned to favorable temperatures. The response to cold in D. melanogaster is phenotypically plastic, labile to artificial selection, and segregates genetically. Drosophila populations from higher latitudes are generally more cold-tolerant (defined as improvement in any cold tolerance trait), 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 conditions. These evolved differences in cold tolerance appear, like phenotypic plasticity, to be partially driven by modifications in hemolymph ion and metabolite composition, adaptive shifts in membrane composition that maintain membrane fluidity at low temperatures, and an ability to maintain metabolic homeostasis in the face of chilling stress.

2Surprisingly, 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 CCRT, while microarray and RNA-seq studies have identified few genes consistently associated with cold tolerance or upregulated following cold exposure in D. melanogaster. Although adult acclimation induces such large shifts in thermal tolerance phenotypes, one proteomic study identified a handful of proteins (largely muscle- and reproduction-related) that differed from controls (25 °C) after five days acclimation at 11 °C. Making direct biochemical links between these existing candidates and cold tolerance phenotypes has also proven difficult. For example, Frost appears to encode a disordered protein, is consistently upregulated after cold exposure, is present in a QTL associated with chill coma, and shows geographic variation in amino acid sequence consistent with environmental variation. However, RNAi knockdown of Frost expression does not consistently impair cold tolerance. Thus, although there is substantial cold tolerance plasticity in D. melanogaster, this plasticity is associated with strikingly few candidate molecules.

3It 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.) or compared the basal transcriptome among flies with detectable, but often small, phenotypic differences (e.g.). 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 signaling. By contrast, longer-term cold acclimation over a period of days or weeks dramatically improves cold tolerance in almost all insects, including Drosophila, and this robust response occurs over a timescale that is consistent with modifications 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.

4Here, we compare D. melanogaster adults acclimated for six days at 21.5 and 6 °C, which leads to a 2.5 °C shift in CTmin in this population (from 3.4 ± 0.2 to 0.9 ± 0.1 °C5). We compare the transcriptomic and metabolomic effects of acclimation temperature, and show that almost one third of the transcriptome and nearly half of the metabolome are differentially regulated with thermal acclimation. Importantly, we found considerable overlap in the metabolic pathways identified 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 differentially expressed, and thus our results both support and extend the hypothesized mechanisms of insect chilling tolerance.

5We 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 50 bp reads (109 M reads in total) were mapped to the D. melanogaster genome from each of five biological replicates of cold- (6 °C) and warm-acclimated flies (i.e. control flies 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 influenced the gene expression profile (Figs 1a and 2); a total of 4,362 genes (c. 29% of coding genes) were differentially expressed between the two acclimation temperatures based on a geometric algorithm and a false discovery rate (FDR)-corrected α cutoff of 0.05 (see methods for details). To focus our efforts on a smaller subset of (highly differentially-expressed) genes we further subset this to genes with greater than a two-fold difference in expression and Q < 0.01 (P-value following FDR), resulting in 1,567 genes that differed in expression between the two acclimation treatments (c. 10% of the genome; hereafter referred to as “differentially 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 identified following cold exposure in D. melanogaster, and in most other insect species.

6The majority (985; 63%) of differentially-expressed genes mapped to at least one gene ontology (GO) term, and 108 GO terms were significantly enriched in the dataset following FDR correction (Dataset 3). To simplify interpretation of this list, we used REVIGO to cluster GO terms across all differentially-expressed genes, which resulted in 49 semantically-distinct GO terms that were significantly enriched in response to acclimation (Dataset 3). We then separated genes by their direction of differential 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 flies and flies maintained at 21.5 °C, respectively (Fig. 3). Complete lists of significant GO terms with differentially expressed genes associated with each term can also be found in the supplementary material (Dataset 3). Several terms for cellular components and processes identified here were also enriched in the previously-sequenced transcriptomes of either D. montana or D. virilis (both in different subgenera to D. melanogaster) following cold acclimation (for example, cytoskeletal components, ion transporters, calcium regulation and glutathione-s-transferase activity). Despite using similar acclimation conditions to our study, Parker et al. associated far fewer candidate genes and GO terms with cold acclimation in these drosophilid species. This difference may arise from analytical differences 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 differed 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.

7Several of the GO terms that were significantly enriched in cold-acclimated flies related to the muscle contractile apparatus (Fig. 3, Dataset 3), consistent with a previous proteomic survey. Low temperatures strongly suppress maximal muscle force production in both fish and insects, an effect that can be partially compensated for through acclimation. Acclimation to low temperatures by poikilothermic fish involves reorganization of muscle fiber type and abundance, a change that is associated with differential expression of several proteins involved in muscle contraction, including the myosin light and heavy chains. Our results suggest a similar reorganization of the muscle fibers 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 flies 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 cell resting potentials and decreased muscle force generation at low temperatures. Warm-acclimated D. melanogaster lose the ability to fly at relatively mild low temperatures (c. 14–18 °C), and our results suggest that flies could modulate muscle force generation in response to low temperature exposure to maintain flight performance even if development at high temperatures has already fixed wing morphology. Calcium regulatory proteins largely determine the thermal sensitivity of the actin-myosin complex, and are also differentially regulated during cold acclimation in fish. Thus, these modifications 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 fiber intersections in the midgut, and chilling disrupts cytoskeletal organization in primary embryonic cultures of Drosophila cells. Thus, regulation of cytoskeletal function is likely to be an important component of cold acclimation in a variety of tissues beyond muscle.

8Transcripts 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. melanogaster. Many genes encoding heat shock proteins were upregulated in cold-acclimated flies, 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. 0 °C for 2h) and/or cold acclimation in D. melanogaster. Expression of all of these candidate cold tolerance genes is highly enriched in the Malpighian tubules, which supports the notion that the renal epithelia, and the Malpighian tubules in particular, are central to insect cold tolerance, 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 conflicting evidence about the role of these candidates in determining cold tolerance phenotypes, as whole-animal knockdowns often have no effect on cold tolerance phenotypes (e.g.). 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 reflects the initiation of thermal acclimation, rather than an adaptive response to acute injury.

9Cold acclimation also involved significant 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 insects. Further, knockdown of chitin synthesis can disrupt the barrier functions of the midgut peritrophic matrix of flour beetles, and the barrier function of the gut epithelia appear to be of great importance to insect cold tolerance. Several genes that were upregulated with cold acclimation in the present study have important roles in tracheal system development and/or may influence cuticular permeability (see Dataset 3). A relationship between cold exposure and immunity has also been suggested by numerous transcriptomic studies, although the nature of this relationship remains unclear, particularly since cold acclimation does not appear to improve realized immunity in crickets. 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 costata.

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

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

12Adenosine mono- di- and triphosphate levels were also all higher in cold-acclimated flies (Figs 1d 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 effect on energy availability in the form of ATP in the fat body of firebugs, and slightly increases ATP availability in the cricket muscles and beetle larvae. By contrast, chronic chilling decreases whole-body ATP levels in flesh fly larvae, a detrimental state that can be rescued by interrupting the cold exposure with pulses of high temperature. These different effects of chilling on ATP availability may be related to whether chilling more strongly affects 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. melanogaster, 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 specific functions of increased ATP availability and the size of the total adenylate pool in cold tolerance remains to be thoroughly examined.

13Although 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 flies 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 flies) 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 females (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). Furthermore, JH stimulates production of antifreeze proteins in larvae of the freeze avoidant beetle (Dendroides canadensis). 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. suzukii. 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.

14Several 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 reflect 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 flies, a shift that is suggested to protect against ion balance disruption at low temperatures. A transporter (Nha1) recently suggested to function as an H+-Cl− cotransporter, rather than a Na+/H+ exchanger as originally suspected, was also downregulated in cold-acclimated flies. This contrasts with the recent finding that knockdown of Nha1 considerably slows chill coma recovery time of D. melanogaster adults. 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.

15A total of ten of the 34 positively-identified metabolites significantly decreased in abundance following cold acclimation (i.e. were higher in flies maintained at 21.5 °C) 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 preference. Another potential explanation for the decrease we observed in dopamine concentration is that dopamine secreted by epidermal cells is oxidized into melanin, and given that Drosophila invest in melanin production during development at low temperatures, it may be that dopamine is being invested into melanin synthesis during cold acclimation. Trehalose was also less abundant in cold-acclimated than warm-acclimated flies. This result is curious, given that trehalose is a well-described cryo- and osmoprotectant, and increases in abundance in response to cold shock or rapid cold-hardening treatments in Drosophila adults, as well as during long-term cold acclimation of Drosophila larvae. 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 species.

16We used MetaboAnalyst and goseq to find metabolic pathways in the KEGG database that were significantly 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 significantly modified during thermal acclimation (Fig. 5a). A total of 283 of the 1,567 differentially expressed genes (18%) from the transcriptomic analysis were represented in the KEGG database. Of these, 103 genes were more highly expressed in cold-acclimated flies, and 180 were more highly expressed in flies maintained at 21.5 °C. We mapped all differentially-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 identification of overlapping pathways (those enriched in both the transcriptomic and metabolomic data), we pooled significant pathways from all three transcriptomic analyses (prior to FDR correction).

17In total, 25 metabolic pathways from the KEGG database were significantly 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 fits 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 temperatures. Indeed, a proline-supplemented diet under winter-acclimating conditions confers partial freeze tolerance to D. melanogaster larvae. In the present study, cold-acclimated flies had higher levels of both proline (1.3-fold) and arginine (1.8-fold) than warm-acclimated flies, 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 differentially expressed during thermal acclimation (Fig. 5b). Together, these findings 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 flies. This discrepancy may arise because proline has different roles in chilling survival versus chill-coma recovery, or because selection for cold tolerance elicits trade-offs 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 response.

18Glutamate 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 differential regulation in both the transcriptome and the metabolome. Interestingly, we found no difference in the abundance of reduced glutathione (GSH), but rather a higher abundance of oxidized glutathione (GSSG) in cold-acclimated flies, relative to flies maintained at their rearing temperature (Fig. 5a, Table S1). This apparent shift in redox state could imply that cold-acclimated flies are under greater oxidative stress at 6 °C, as cold stress and chilling injury have been associated with oxidative damage and rightward shift in the GSH:GSSG ratio, particularly in chill susceptible and freeze avoiding insect species. In the transcriptome of D. melanogaster, we found 17/39 (29%) of genes mapping to the Glutathione Metabolism term were differentially 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 Lepidoptera, and GSH:GSSG balance plays important roles in cold stress and cold acclimation signaling in plants, suggesting that this observed shift might serve a signaling, as well as protective, role.

19Cold acclimation is thought to improve insect cold tolerance through an improved ability to maintain homeostasis and cellular integrity at low temperatures. The primary mechanisms proposed to underlie this phenotypic plasticity are 1) protection against cold-induced ion and water balance disruption, 2) defense against oxidative stress, 3) maintenance of metabolic balance, 4) refolding or stabilizing proteins, 5) modifications to membrane fluidity to maintain cellular integrity, and 6) inhibition of apoptosis.

20Our results support most, but not all, hypothesized mechanisms of chill tolerance. While the whole-animal metabolome cannot offer substantial insight into adjustments to ion and water balance, a large number of genes involved in ion homeostasis were differentially-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 modifications 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 specifically suggest changes to membrane lipid metabolism specifically. 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 shifts in Ca2+ handling) to suggest changes to apoptotic signaling pathways during long-term adult thermal acclimation.

21Our 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 shifts 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 modifications to thermal tolerance occurring over these longer time scales. As roles of specific 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 phenotypes are now becoming apparent (e.g.), we suggest ‘omics’ approaches can, and should, now be used in a more hypothesis-driven and tissue-specific 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-specific and whole-organism thermal tolerance phenotypes.

22Temperature affects 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 fitness. Here, we generated a dataset that finally reflects 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).

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

24Newly-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 left at their respective acclimation temperatures for five days to acclimate and to control for age and anaesthesia effects on cold tolerance, before being snap-frozen in liquid nitrogen and held at −80 °C for later analysis of gene expression and metabolite abundance.

25Transcript abundance was quantified using high-throughput mRNA sequencing (RNA-seq). Frozen warm- and cold-acclimated male D. melanogaster (five biological replicates per acclimation temperature, each containing 25 flies) 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 200 μL (lipids and DNA) discarded. Chloroform (200 μL) was added to the tubes, samples were shaken for 15 s, and left 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 (30 s). For further sample purification, 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 2 mL collection tube. The effluent was discarded and the sample (bound to the column) was washed with 700 μL of RW1 buffer (Qiagen), and 2 × 500 μL of RPE buffer (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 (8000× g for 1 min at 22 °C) twice to release the RNA from the column. Total mRNA was purified 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). Briefly, mRNA was purified 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 (13–14 samples per lane, run alongside samples from two other experiments) single-end 50 bp cDNA libraries at the High Throughput Genomics Core Facility at the Huntsman Cancer Institute, University of Utah.

26Transcriptome data manipulation and statistical analyses were completed in Galaxy following the analysis workflow of Trapnell et al.. Briefly, 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.8, with alignment limited to known splice junctions. Cufflinks was used to assemble and count transcript reads for each sample and Cuffmerge v.2.0.2 was used to merge assembled transcripts for all biological replicates into a single reference transcriptome. To test for differential expression of genes that mapped to the transcriptome, we used CuffDiff 2, which uses a geometric-based algorithm to generate test statistics, P-values, and Benjamini-Hochberg false discovery rate (FDR) corrected P-values (Q-values) for each gene mapped. More than 4,300 genes were significantly differentially expressed with thermal acclimation after FDR correction (Q < 0.05), so to reduce the number of targets we set cutoff values for the magnitude and statistical significant of differential expression. Genes that had both a 2-fold difference in expression between the acclimation groups, and a Q-value less than 0.01 were used in all downstream analyses and were deemed “differentially expressed genes”. Differentially expressed genes were mapped to the Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases using goseq in R. To reduce the number of significant GO terms for plotting purposes, we used REVIGO, which eliminates semantically redundant terms.

27Metabolites were extracted from whole fly samples (n = 40 per acclimation temperature) following Knee et al., with minor modifications as described below. Samples consisted of groups of 15 male flies that were weighed to the nearest 0.01 mg and flash frozen in liquid nitrogen. The volume of extraction solvent used for each sample was normalized by the mass of the sample by pipetting 6.35 μL extraction solvent per mg of fly. 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 −80 °C for later analysis.

28Metabolomic analysis was carried out as previously described. Briefly, 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 (10 μL) 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 flow-rate of 0.4 mL/min (Table S2). Data was acquired at 2 Hz in negative ionization mode.

29Raw metabolite data was subject to automatic intra-run calibration with sodium formate and putative metabolites identified using the Find Molecular Features peak picking software (Bruker Daltonics). Principal Component Analysis (PCA; with Pareto scaling) was performed in R, and and T-tests were performed on ProfileAnalysis software (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.

30Pathway Analysis, a combination of quantitative enrichment analysis and topology analysis, was carried out on MetaboAnalyst online software using characterized KEGG pathways for D. melanogaster as the back-end knowledge. Enrichment analysis is based on concentration values of compounds loaded into the analysis and uses a Global Test to determine if pathways are significantly 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 pathway.

31How 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).