Introduction

Insects have a range of tactics to deal with the onset of seasonal cold (Salt, 1961). One widespread tactic, cold acclimation, is an increase in an organism’s cold tolerance that occurs over days or weeks when exposed to sub-lethal low temperatures, allowing insects to remain active for longer as the seasons change (Lee, 2010). The physiological changes leading to increased cold tolerance during cold acclimation have been well studied in insects and are known to primarily involve a shift in the metabolic profile of the insect, particularly sugars, polyols and amino acids (Lee, 1991; Lee, 2010). Accompanying this is a compositional change of the cell membrane primarily involving changes in phospholipids (Hazel, 1995; Koštál et al., 2003; Overgaard et al., 2008). These changes are cryoprotective, allowing cells to maintain their osmotic balance and thus to continue to function at low temperatures.

Cold acclimation is likely to be a homologous trait for many species of closely related insects as it has clear fitness benefits (Salt, 1961; Lee, 2010) and is thus unlikely to be lost and gained independently. Homology can be defined at several levels of biological organisation including genes, gene networks, and phenotypes, which may not necessarily coincide as the different levels may evolve at different rates. As a result it could be that although a trait is homologous and phenotypically similar across species, the genes that contribute to the production of the trait could be different as a result of evolutionary turnover, whereby genes involved in similar biological processes are co-opted in or out of influencing a trait as a result of orthologous selection or drift. This process of evolutionary turnover can thereby allow the genetic basis of a trait to change, whilst the trait itself remains conserved (see Canestro et al. (2007) and McCune and Schimenti (2012) for reviews).

Candidate gene approaches, whereby genes of interest identified in a model species are functionally examined in others, have often been successful and have indicated that gene function tends be conserved across species (Fitzpatrick et al., 2005; Reaume and Sokolowski, 2011; Martin and Orgogozo, 2013). In contrast, quantitative genetics studies have often suggested that traits important to adaptation are polygenic, and quantitative trait loci (QTL) may have low repeatability between species (Arbuthnott, 2009; Huang et al., 2012; Rockman, 2012). One reason for this disparity may be that candidate genes are typically those that have a large influence on a particular trait, and thus the evolutionary forces affecting them are likely to be stronger and more consistent than those for smaller effect loci. In particular, we expect that genes with a large influence on a particular trait are more likely to be conserved in a particular role (by purifying selection) than those with a smaller influence, which can be more easily co-opted in or out of influencing a trait.

As cold acclimation is an induced response, the gene expression changes involved with the trait can be studied using an RNA-seq approach. RNA-seq allows examination of the changes in gene expression independent of the strength of influence they have on cold acclimation, identifying both genes that have a major controlling influence on cold acclimation and genes involved in smaller secondary changes. Furthermore, it allows assessment of the proportion of genes that respond to cold acclimation that are the same between the two species, and what biological processes they are involved in.

Here we investigated gene expression changes in response to cold acclimation in Drosophila (D.) montana and D. virilis using an RNA-seq approach. These species belong to the D. virilis species group, and diverged from each other approximately 9–11 mya (Reis et al., 2008; Morales-Hojas et al., 2011). Both species show a high level of cold tolerance compared to other drosophilids, and D. montana is significantly more cold tolerant than D. virilis (Kellermann et al., 2012; Vesala et al., 2012b). These species have adapted to live in quite different environmental conditions (D. montana is found at high latitudes (30–70°N) and altitudes while D. virilis is a human commensal found at lower latitudes (south from 35°N) (Throckmorton, 1982)) but both exhibit a similar cold acclimation phenotype, whereby flies kept at sub-optimal temperatures show increased cold tolerance as measured by chill-coma recovery time (Vesala et al., 2012b). When insects are placed in cold (<0 °C) temperatures a chill-coma is induced whereby insects lose the ability to move temporarily, due to a loss of nerve and muscle excitability (MacMillan and Sinclair, 2011). The time it takes for an insect to regain its ability to move is referred to as chill coma recovery time, which can be used as an estimate of an insect’s cold tolerance (see David et al. (1998)). The cold acclimation responses of D. montana and D. virilis are likely to represent homologous traits as this trait is present in both species and across Drosophila species in general (see Hoffmann et al. (2003)). Thus our aims were: (i) to identify the genes that are differentially expressed in response to cold acclimation in D. montana and D. virilis and (ii) to identify the molecular pathways and processes underlying the cold acclimation response in each species, allowing us to determine to what extent to which these are consistent between the species.

Materials and Methods

Samples

We used a single isofemale line established from the progenies of wild-caught fertilized females for each species: D. montana line 175OJ8 (originated from Oulanka, Finland, 2008) and D. virilis line TOY3F9 (originated from Toyama, Japan, 2003). The use of isofemale lines will reduce confounding effects of genetic variation within lines and cold acclimation responses (see discussion). Stock cultures have been maintained since their establishment on malt bottles in continuous light at 19±1 °C, 65% humidity. Flies for the experiment were collected within one day of eclosion using light CO2 anesthesia. Female flies were put into malt medium vials and transferred into a climate chamber (Sanyo MLR-351H) for 15 days at 16 °C for D. montana and 19 °C for D. virilis; see rationale for the temperature selection in discussion. Flies of each species were then split into two groups: a control group where flies were left in the same conditions for an additional 6 days, and a cold acclimation group where flies were maintained for 6 days at +5 °C. The light:dark (LD) cycle was 22:2 for all experimental groups and the samples for RNA extraction were flash frozen in liquid nitrogen 5 h after lights were turned on in the chambers (Zeitgeber=5).

RNA extraction and sequencing

Frozen flies were pooled into 6 samples for D. montana (3 acclimation group samples and 3 control group samples, with 10 whole flies in each sample) and 4 samples for D. virilis (2 acclimation group samples and 2 control group samples, with 20 whole flies in each sample). Different numbers of pooled individuals per sample is unlikely to influence our analysis as all flies are from isofemale lines, but this design could have more power to detect changes in D. montana than D. virilis (but, see later). RNA was extracted from each sample using Tri Reagent (Sigma-Aldrich) followed by RNeasy Mini kit purification with DNase treatment (Qiagen). Purity of the RNA was checked using NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies) and integrity with 2100 Bioanalyzer (Agilent Technologies). Total RNA for each sample was approximately equal prior to sequencing.

Extracted RNA was sequenced using the SOLiD platform. For D. montana we used the SOLiD 5500 XL to produce 46 million 75+35 bp paired end reads and for D. virilis we used SOLiD V4 to produce 49 million 50 bp single end reads. Raw sequence reads were then trimmed using SOLiD TRIM (with run options: -p 3 -q 22 -y y -e 2 -d 10) to remove polyclonal errors from the data (Sasson and Michael, 2010). The reads that passed this filter were then error corrected using SOLiD Accuracy Enhancer Tools (SAET) to reduce the amount of color calling errors, or erroneous bases, in the sequence. Remaining low quality bases at the end of the reads were then trimmed using CLC Genomics Workbench 5.0.1 (CLC Bio http://www.clcbio.com/) (quality score: 0.02). Differences in the number of replicates and read type used between samples could influence our power to detect differential expression (DE), however we expect these issues to be small as we used a single model to detect DE and interactions. In addition we also repeated the analysis for D. montana using two replicates and with reads trimmed to be 50 bp and single end only (to make the samples comparable to those obtained for D. virilis). This produced very similar results as those with three replicates and paired reads (numbers of genes found to be DE 162 vs 177; correlation of the gene expression for DE genes from both analyses (rs)=0.94, P<2.2 × 10−16, see Supplementary Materials 1 for more details).

Mapping

Reads for each sample were mapped individually to the D. virilis genome (r1.2, available from http://flybase.org/) using CLC Bio. We also used this method to map reads to a de novo assembly of the D. montana transcriptome (see Supplementary Materials 2 for more details). HTSeq (Anders et al., 2014) was used to quantify the read counts mapping uniquely to the reference using the reference annotation available from Flybase (r1.2).

Expression analysis

Gene expression analysis was performed using the Bioconductor package EdgeR (Robinson et al., 2010) in R (R Core Team, 2013). Normalisation factors for each sample were computed using the TMM method. TMM normalisation allows variation in read depth (due to RNA quality, variation in sequencing reaction, etc.) to be accounted for, to prevent differences in read depth from influencing the detection of DE (Robinson and Oshlack, 2010). We then fitted a generalized linear model (GLM) with negative binomial distribution with the terms species, treatment and species * treatment (full model), and estimated dispersion using the Cox-Reid profile-adjusted likelihood (CR) method. We used a GLM likelihood ratio test to determine significance of a treatment effect for each gene by comparing appropriate model contrasts: the effect of treatment on D. montana, the effect of treatment on D. virilis and an interaction between species and treatment. The interaction term tests the extent to which gene expression changes differ between the species. The P-values from the GLM likelihood ratio tests were corrected for multiple testing using Benjamini and Hochberg's algorithm to control for false discovery rate (FDR) (Benjamini and Hochberg, 1995) with significance taken here to be <5% (FDR<0.05).

Functional classification

In order to functionally classify genes we used Gene Ontology (GO) annotation for orthologous genes in D. melanogaster (available from Flybase, version: FB2013_06). We used this approach instead of using GO terms from the D. virilis annotation due to the superior GO annotation available in D. melanogaster both in terms of the number of annotations and their specificity (Tweedie et al., 2009). Significant enrichment of single GO terms were determined using the ‘Gene Ontology Enrichment’ function in FlyMine (www.flymine.org, v.37). The D. melanogaster orthologs of DE genes were also analysed using DAVID (Database for Annotation, Visualization and Integrated Discovery) v. 6.7 (Huang et al., 2009a; Huang et al., 2009b). DAVID clusters genes into functional groups using a ‘fuzzy’ clustering algorithm, and then uses a Fisher’s exact test to identify significantly enriched functional groups. A functional group was considered to be significantly enriched if its enrichment score (the geometric mean (in -log scale) of the P-values of the GO terms in the group) was >1 (P<0.1).

Results

We obtained approximately 46 million reads for D. montana and 49 million reads for D. virilis, of which approximately 10% were found to map uniquely to a gene in the D. virilis reference genome. The number of reads mapping to each gene was highly correlated between the species (rs=0.90), suggesting the mapping efficiency per gene is approximately equal for the two species. The number of genes that were differentially expressed (DE) during cold acclimation was fewer in D. montana (177) than in D. virilis (458) at 0.05 FDR (representing 593 total unique genes DE across both species, see Supplementary Materials 3). Among these genes only 42 genes (7%) were DE in both species. The genes DE in both species had very similar expression changes in response to cold acclimation in terms of direction (rs=0.85, P<2.2 × 10−16) (Figure 1a) and fold change (Figure 1b) though changes were on average somewhat higher for D. virilis. One gene (Dvir\GJ10437) showed a significant interaction effect, being up-regulated in D. montana but down-regulated in D. virilis. The remaining genes, i.e. the genes that were DE in only one of the species, showed a much smaller correlation in expression levels (rs=0.20, P=1.18 × 10−6) (Figure 2) and 132 (24%) of them showed a significant species by treatment interaction, whereby genes which showed a strong response to acclimation in one species do not show a response in the other (Figure 3) (i.e. these genes showed species-specific responses).

Figure 1
figure 1

(a) Plot of log2 fold changes between control and cold acclimated flies for each of the genes that were DE in both D. montana and D. virilis (b) Average change in gene expression (Δ log2 fold change) between control and cold acclimated flies for genes that were DE in both species. Error bars represent approximate 95% confidence intervals.

Full size image
Figure 2
figure 2

Plot of log2 fold changes between control and cold acclimated flies for each of the genes that are DE in either D. montana or D. virilis.

Full size image
Figure 3
figure 3

Log2 fold change between control and cold acclimated flies for the top 20 genes which showed a significant species by treatment interaction in D. montana and D. virilis. D. virilis gene names are given in the legend with D. melanogaster orthologs in parentheses.

Full size image

Gene function

To examine the biological processes and pathways of the genes DE due to cold acclimation, we examined the genes in four subsets: those DE in D. montana, those DE in D. virilis, those DE in both species, and those that showed a significant treatment by species interaction. Broad level (2) GO term functional classification revealed that the majority of genes that showed significant differential expression due to cold acclimation in each of the subsets were involved in metabolic and cellular processes. We also found that the proportion of genes involved in a particular biological process was similar for each of the subsets (Figure 4). When the metabolic GO term was split into its constitute parts (Figure 5, Supplementary Materials 4–8) we found the proportion of genes annotated with each metabolic GO term was similar for each of the subsets.

Figure 4
figure 4

Proportion of genes annotated for each of the level 2 GO terms identified for each of the gene subsets (D. montana, D. virilis, both species and interaction). Black dashed lines indicate proportion of genes annotated for each of the level 2 GO terms in the D. virilis genome.

Full size image
Figure 5
figure 5

Proportion of genes annotated for each of the primary metabolic processes for each of the gene subsets (D. montana, D. virilis, both species and interaction). Black dashed lines indicate proportion of genes annotated for each of the primary metabolic processes in the D. virilis genome.

Full size image

GO term enrichment analysis

The single GO terms oxidation-reduction process, and single-organism metabolic process were significantly enriched in D. montana (P=0.002 and P=0.016 respectively). No single GO terms were significantly enriched in D. virilis. By using the functional clustering program DAVID we were able to identify significantly enriched functional clusters in three of the four subsets (D. montana, D. virilis, and interaction (Table 1)) but none were identified for the genes DE in both species. Significant clusters involved in metabolism or producing metabolites were found in the D. montana, D. virilis, and interaction subsets (Table 1, yellow) as well as in transmembrane transport/ion transport in D. montana and D. virilis (Table 1, green). 14 (8%) of the genes DE in D. montana and 39 (9%) of the genes DE in D. virilis were annotated with the GO terms for transmembrane transport/ion transport, with 2 genes DE in both species and 9 showing a significant interaction (Supplementary Materials 9).

Table 1 Functional clusters identified by DAVID for each of the gene groups for enrichment scores with P<0.1
Full size table

In addition we identified a number of enriched gene clusters which were different between the species. In D. montana, we found two gene clusters which showed significant enrichment for muscle protein (Table 1, orange), while in D. virilis there were several clusters significantly enriched for protein signalling (Table 1, pink). Although the above mentioned functional clusters were only significantly enriched in one species, these processes are likely to be a component of cold acclimation in both of the species as we also found genes annotated with these functions DE in both species (Supplementary Materials 10).

Discussion

While our understanding of the genetic basis of many traits has advanced greatly in many model systems, whether these results translate across species remains an open question. Numerous candidate gene studies have demonstrated that genes appear to be conserved in their function between taxa (Fitzpatrick et al., 2005; Reaume and Sokolowski, 2011; Martin and Orgogozo, 2013). Candidate gene studies however generally focus upon a subset of genes which have a large influence on a particular trait. By using an RNA-Seq approach we were able to examine genes which altered expression in response to cold acclimation independent of the strength of influence they have on cold acclimation in two closely related Drosophila virilis group species where the trait is likely to be homologous. This approach captures both the genes that have a major influence on cold acclimation but also genes of smaller influence and those involved in secondary changes, which may be more susceptible to evolutionary turnover. Using this approach enabled us to identify genes that show differential expression in response to cold acclimation in both or only one of the species and to identify the molecular pathways and processes involved. Our study found evidence for both conservation and divergence in gene expression in response to cold acclimation between D. montana and D. virilis.

We identified 42 genes that were differentially expressed in both species. These genes (with one exception) showed very similar changes in terms of fold change and direction of expression changes, and thus may represent a ‘core set’ of genes, which appear to be evolutionary conserved in response to cold acclimation. We also identified a large number of genes that were DE in one species but not the other with around a quarter of these showing a significant species by treatment interaction. Gene expression changes in these genes were only weakly correlated with each other, with only one of the 132 genes which showed a significant interaction DE in both species. This suggests that a large proportion of the genes that were DE in response to cold acclimation were different between D. montana and D. virilis. We also observed that many more genes were DE in D. virilis than D. montana. One potential reason for this may be that D. virilis is less cold tolerant than D. montana (Kellermann et al., 2012; Vesala et al., 2012b) and thus the cold acclimation treatment may be more stressful for D. virilis than D. montana, promoting a stronger cold acclimation response. In support of this we also note that the average fold change in genes DE in both species was also slightly larger in D. virilis than D. montana (Figure 1b).

Our focus in this study was to look at temperature-induced changes in gene expression in D. montana and D. virilis. Therefore we chose the conditions for the control treatments to represent non-cold-acclimating ‘summer’ conditions for each species. A complication is that identical long day conditions for both species may have led to D. montana flies entering reproductive diapause which would have complicated the results markedly as D. virilis does not have such photoperiodic diapause (Throckmorton, 1982). According to temperature data obtained for the collection site of D. montana flies used in this experiment (years 2000–2009, Oulanka research station, University of Oulu, Finland), the average temperature of summer months was: June: 12.2 °C, July: 15.5 °C and August 12.6 °C. Based on this data we chose 16 °C to represent non-acclimating (control) temperature for D. montana strain and 19 °C for more southern and less cold tolerant D. virilis. As such, the temperature difference between the control (non-acclimating) and the acclimation temperatures was slightly larger for D. virilis than D. montana (14 and 11 °C respectively) which could partly influence the larger number of genes DE in D. virilis compared to D. montana. However, given the relatively small difference between the two control non-acclimating temperatures compared to the large difference between the control and the acclimation temperatures, we expect the potential influence of this to be small.

As we were primarily interested in species-specific responses to cold acclimation we used isofemale lines for both species to minimise confounding intra-population genetic variation within species. However, as we have only one line per species we cannot distinguish intra- and inter-specific sources of genetic variation. Yet, as the species separated 9–11 mya (Reis et al., 2008; Morales-Hojas et al., 2011) it seems most likely that interspecific differences will dominate the main differences described here. Future population-level work is required to explore the extent of differences within species.

Functional processes

As we found that the genes DE due to cold acclimation were often different between the species, we examined the functional processes involved in four subsets: genes DE in D. montana, in D. virilis, in both species, and genes which showed a significant treatment by species interaction. We discuss the main functional processes found in the DE gene subsets below (for full list of processes for all the genes found to be DE see Supplementary Materials 11).

Metabolic profile

Cold acclimation is known to involve a shift in the metabolic profile as well as the production of cyroprotectants which act to maintain osmotic balance and stabilise the membrane structures of a cell (Lee, 1991; Lee, 2010). Consistent with this previous work, we found that the majority of genes DE in both species were metabolic in function. Splitting the metabolic GO term into its constitute parts showed that there were a similar proportion of genes involved in each of the metabolic processes for the studied subsets (Figure 4). The single GO-terms ‘oxidation-reduction processes’ and ‘single-organism metabolic processes’ were significantly enriched in D. montana. We did not find any significantly enriched single GO terms for D. virilis, though D. virilis showed a similar (though slightly lower) proportion of metabolic genes involved in the cold acclimation response. Examination of the enriched functional clusters identified clusters of genes involved in metabolism/production of metabolites for three of the subsets (D. montana, D. virilis, and the species by treatment interaction, Table 1). We did not find any significantly enriched functional clusters for genes DE in both species, likely due to the small number of genes in this group. Although we identified a number of metabolic gene clusters for three of the subsets, the metabolic pathways implicated for each set were different. This suggests that although both species show significant enrichment for metabolic changes, the exact metabolic pathways involved differ.

Cell membrane composition

Changes to the composition of phospholipids in the cell membrane are thought to be particularly important for cold acclimation, as it allows cells to maintain their physiological function in sub-optimal temperatures (Hazel, 1995; Koštál et al., 2003). Our previous work has shown the major metabolite produced in overwintering D. montana flies to be myo-inositol (Vesala et al., 2012a), which functions in many processes including regulation of cell development and growth (Loewus and Loewus, 1983). It also is a precursor of inositol phospholipids leading to production of inositol phosphates that function as second messengers (Downes and Macphee, 1990). Myo-inositol has not been found to play a role in cold acclimation in other Drosophila species, but it has been shown to accumulate seasonally in few other insect species (Block and Sømme, 1983; Koštál et al., 1996). It has also been correlated with an increase in cold tolerance in some Coleoptera and Lepidoptera species (Watanabe and Tanaka, 1999; Watanabe, 2002; Liu et al., 2009) and we hypothesise that its accumulation in overwintering D. montana would act as a cryoprotectant. In the present study the ortholog of inos (Dvir\GJ20549) was upregulated in both D. montana and D. virilis. This is interesting as inos encodes the enzyme myo-inositol-1-phosphate synthase, which is part of the inositol biosynthetic pathway (GO:0006021), catalysing the conversion of D-glucose-6-phosphate into 1L-myo-inositol-1-phosphate, the first committed step into production of all inositol compounds (Majumder et al., 1997). It is not known what the major metabolite produced in overwintering D. virilis flies is but given the DE of inos ortholog in response to cold acclimation we suggest that the production of 1L-myo-inositol-1-phosphate is likely to be important for both species. Another gene annotated in the inositol metabolic pathway, Dvir\GJ15346 (D. melanogaster ortholog: CG6910), showed a significant species by treatment interaction during acclimation. This gene, which was significantly down-regulated only in D. montana, is annotated with the GO term inositol catabolic process and has been indicated to have the same function as myo-inositol oxygenase, an enzyme that is the first committed step in inositol catabolism in eukaryotes (Jones et al., 2012). Thus the downregulation of this gene may lead to higher concentrations of inositol compounds, enabling D. montana to accumulate proportionally more inositol compounds, than D. virilis. D. virilis, on the other hand, may increase its cold tolerance with the aid of sphingolipid and sterol compounds, which have previously been shown to act together to alter cell membrane fluidity (Guan et al., 2009), as genes annotated for the metabolism of these compounds were found to be DE in D. virilis but not in D. montana (Supplementary Materials 5).

An additional group of genes that may influence cell membrane composition are Niemann-Pick type C (npc) genes. In D. melanogaster npc genes have been shown to be involved in the homeostatic regulation of cholesterol, which influences the permeability and fluidity of cell membranes (Huang et al., 2007; Niwa and Niwa, 2011). We found that npc2e was DE in both species, though it was upregulated in D. montana and downregulated in D. virilis (interaction P-value<0.001). In addition, npc1b and npc2d were DE in D. virilis but not in D. montana (interaction P-values for both genes<0.001). This suggests that the changes in the regulation of cholesterol are important for cold acclimation in both species, though the genetic basis used to produce the change is different in each of the species.

Ion transport/transmembrane transport

Changes in temperature are known to affect the transport mechanisms involved in the maintenance of cellular ion balance (Heitler et al., 1977; Kivivuori et al., 1990). Failure to maintain the ionic balance of cells can lead to metabolic perturbations which can cause a wide range of negative consequences, in particular, the loss of nerve excitation (Hochachka, 1986; Kostal et al., 2004). We found enrichment of several functional gene clusters involved in ion transport/transmembrane transport in both species (Table 1) which suggests that changes in gene expression of these genes are acting to maintain cellular ion concentration at low temperatures. Although both species showed enrichment for ion transport/transmembrane transport the genes responsible were largely different, with only two genes shared between the species and nine showing a significant interaction (Supplementary materials 9).

Rhodopsin

The ortholog of a key gene in the rhodopsin bio-synthesis pathway santa-maria, Dvir\GJ17608, was up-regulated in both species in response to cold acclimation. Rhodopsin has long been known to be the primary pigment for phototransduction (see Katz and Minke (2009) for a review), however recently it has been shown that the rhodopsin signalling pathway also has a several additional light-independent roles including hearing (Senthilan et al., 2012) and thermosensory signalling (Shen et al., 2011). Shen et al. (2011) showed that by knocking out santa-maria in D. melanogaster, flies were unable to discriminate between differences in temperature. Our finding that Dvir\GJ17608 is upregulated in both species is intriguing as it suggests the rhodopsin pathway may act to detect changes in temperature and thus may cue the cold acclimation response. None of the other genes involved in the rhodopsin bio-synthesis pathway were found to be DE in either of the species suggesting that santa-maria may be the key gene involved in the detection of temperature change. We did, however, find that an ortholog of CdsA (Dvir\GJ13151), a gene involved in the regulation of lipid storage, was DE in D. virilis but not D. montana. CdsA also belongs to a rhodopsin mediated signalling pathway and thus could play a role in linking rhodopsin signalling and downstream metabolic changes.

Circadian clock

Consistent with previous candidate gene approaches in D. montana and D. virilis (Vesala et al., 2012b) we found that two of the core circadian clock genes, period and vrille, were DE in both species. These genes have a well-studied role in the regulation and maintenance of circadian rhythms in insects (Bell-Pedersen et al., 2005). In addition, we also found four other peripheral clock genes to be DE in either D. montana or D. virilis including one which showed a significant species by treatment interaction (Table 2). The circadian clock has been previously implicated in cold acclimation in several taxa (Fowler and Thomashow, 2002; Magnone et al., 2005; Vesala et al., 2012b). It is not clear however if the changes in these genes have a direct influence on the increase in cold tolerance or if changes of the expression levels of these genes affect changes to the circadian rhythms of the fly in anticipation of the change in seasonal photoperiod. Further functional genetic studies will allow disentanglement of these hypothesis but given extensive work showing how the clock genes have a large influence on the regulation of metabolic processes in Drosophila (Xu et al., 2008; Xu et al., 2011; Sahar and Sassone-Corsi, 2012) we suggest that clock genes are good candidates for orchestrating the shifts in metabolic profile seen during cold acclimation.

Table 2 Circadian genes differentially expressed in response to cold acclimation
Full size table

Differences between the species

Despite the overall pattern of broadly similar processes occurring in both species we also identified differences between the species. In D. montana we found significant enrichment for actin filaments/muscle protein, while in D. virilis we found significant enrichment for several clusters of genes involved in protein signalling/modification (Table 1). Changes in actin filaments /muscle protein have previously been implicated in insects to prevent cold injury by preventing the depolymerisation of the cell membrane during diapause (Kim et al., 2006; Robich et al., 2007; Koštál, 2010). Thus, these changes may help to maintain cellular physiology at sub-optimal temperatures. The four clusters of genes involved in protein signalling/modification indicate that these processes may also be important in cold acclimation. These differences seen between the species are unlikely to represent the absence of one process in one species and presence in another, but rather a greater or lesser relative importance in the context of the other processes occurring during cold acclimation for that species. For each of these clusters we also found several genes in the other species involved in the same biological process (Supplementary Materials 10) suggesting that these processes are likely to be a component of cold acclimation in both of the studied species.

Cold acclimation and cold shock

Previous work on cold tolerance in insects (particularly in D. melanogaster) has mostly focused upon the response to short-term cold shock rather than longer periods of cold acclimation (Colinet and Hoffmann, 2012). These processes, although related, are distinct and likely to produce differing transcriptional responses. This is important to consider when comparing our findings to previous work. Broadly, work on short-term cold shock gene expression changes has shown that the predominant group of genes showing expression changes are those involved in the stress and immune response (e.g. heat shock proteins (Hsps), Turandot genes, etc) (Zhang et al., 2011; Sinclair et al., 2013; Storey and Storey, 2013). In addition, the neuropeptide CAPA has been shown to play an important role in both cold and desiccation resistance in several species of Drosophila, including D. montana, when assessed using cold shock (Terhzaz et al., submitted). We do not find these genes or others involved in heat-shock response to be DE during cold acclimation. This is consistent with the findings by Colinet et al. (2013) who did not detect changes in heat shock protein levels after a 5-day cold acclimation in D. melanogaster. One explanation for this may be that the Hsps are important in the early part of thermal stress but are not involved in longer-term cold acclimation per se (Colinet and Hoffmann, 2012; Teets and Denlinger, 2013; but see Vesala et al. 2012b).

Functional modules

We were able to identify several functional processes that contribute to the cold acclimation response, the most important of which are likely to include: cold detection, circadian genes, and metabolic and cellular membrane profile shifts. For each of these functional processes we found genes that are conserved between the species and others that have diverged. We suggest that the genes that were found to be DE in both species are likely to represent key genes with a large influence on the cold acclimation phenotype. Further functional genetic studies are needed to confirm this, although it should be noted that several of the genes DE in both species have previously been demonstrated to have a large influence on processes implicated in cold acclimation. For example, we found that orthologs of two of the major genes involved in the circadian clock (period and vrille) were DE in both species, whereas genes more peripheral to the clock were DE in one species but not in the other (Table 2). We also observed the ortholog of the gene inos, which is responsible for producing the main overwintering metabolite in D. montana, and thus presumably has a large influence on the cold acclimation, is DE in both species. Interestingly, we found santa-maria (Dvir\GJ17608) was DE in both species. As stated above, santa-maria has a large role in the rhodopsin bio-synthesis pathway (Shen et al., 2011). Although this pathway has been previously linked to the ability to discriminate between differences in temperature (Shen et al., 2011), it has not been implicated in the cold acclimation response. One possibility is that rhodopsin bio-synthesis pathway may act to detect the onset of cold and to subsequently cue the DE of genes which increase cold tolerance.

Conclusions

We examined the effect of cold acclimation on gene expression in two Drosophila virilis group species to determine what genes and pathways were involved and if these were the same for both species. We found that the transcriptional changes associated with cold acclimation were broadly different with only 42 genes DE in both species and 132 showing a significant species by treatment interaction. This suggests that when comparing homologous phenotypes the underlying genetic basis of the trait may differ. The genes that were DE in both species are likely to have been conserved in their role in cold acclimation since at least the time of the species split. These evolutionarily conserved genes are likely to have a large influence on cold acclimation, as the selective forces maintaining cold acclimation will be stronger on these genes, than on genes that have a smaller influence.

Although we found that many of the genes DE due to cold acclimation in each of the species were different, we also found that the biological processes they were involved in were broadly similar. This is consistent with the idea of evolutionary turnover co-opting functionally related genes of smaller influence in and out of the cold acclimation response. The result of such evolutionary turnover is that although a trait may be homologous and appear phenotypically similar between species, many of the genes underlying the trait can be different.

Data Archiving

All data has been deposited in NCBI’s Gene Expression Omnibus (GEO). Accession ID: GSE58925 Databank URL: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58925