Skip to main content

Effects of adult temperature on gene expression in a butterfly: identifying pathways associated with thermal acclimation



Phenotypic plasticity is a pervasive property of all organisms and considered to be of key importance for dealing with environmental variation. Plastic responses to temperature, which is one of the most important ecological factors, have received much attention over recent decades. A recurrent pattern of temperature-induced adaptive plasticity includes increased heat tolerance after exposure to warmer temperatures and increased cold tolerance after exposure to cooler temperatures. However, the mechanisms underlying these plastic responses are hitherto not well understood. Therefore, we here investigate effects of adult acclimation on gene expression in the tropical butterfly Bicyclus anynana, using an RNAseq approach.


We show that several antioxidant markers (e.g. peroxidase, cytochrome P450) were up-regulated at a higher temperature compared with a lower adult temperature, which might play an important role in the acclamatory responses subsequently providing increased heat tolerance. Furthermore, several metabolic pathways were up-regulated at the higher temperature, likely reflecting increased metabolic rates. In contrast, we found no evidence for a decisive role of the heat shock response.


Although the important role of antioxidant defence mechanisms in alleviating detrimental effects of oxidative stress is firmly established, we speculate that its potentially important role in mediating heat tolerance and survival under stress has been underestimated thus far and thus deserves more attention.


Most organisms are faced with heterogeneous and fluctuating environmental conditions, warranting the capability to adjust the expression of phenotypes to prevailing environmental conditions [1, 2]. Such adjustments may result from genetic adaptation or phenotypic plasticity, with the latter comprising a fast and effective means to alter fitness-related traits [3, 4]. Consequently, adaptive plastic responses are considered to be an important aspect of population adaptation [4, 5]. Temperature is an important environmental factor warranting phenotypic adjustment [4, 5]. In butterflies, for instance, temperature-induced adaptive plasticity has been found in a variety of traits including temperature stress resistance [6,7,8,9,10,11]. Regarding the latter, there is evidence across several insect taxa that cold tolerance increases after exposure to cool and heat tolerance after exposure to warm temperatures [9, 10, 12, 13]. In recent years, temperature-mediated plasticity in stress tolerance has gained renewed interest, owing to concerns about the potential impact of anthropogenic climate change on extant biodiversity [4, 5, 14]. While phenotypic plasticity may be induced during development or in the adult stage, we here exclusively focus on acclimation to different adult temperatures [10, 15,16,17,18,19,20]. This was motivated by the facts that adult temperature induces clear variation in thermal stress tolerance, and in order to exclude potentially confounding effects of the developmental environment.

To gain a more complete understanding of the prospects and limitations of beneficial acclimation in buffering effects of temperature variation, a better understanding of the underlying mechanisms is important. At the mechanistic level, phenotypic plasticity in thermal tolerance may be caused by environmental effects on gene and protein expression or biosynthetic pathways [21, 22], but currently our understanding of the underlying mechanisms is very limited [15, 23,24,25,26,27]. Recent developments in sequencing techniques though have opened up new opportunities to address this issue, enabling unbiased approaches by targeting entire genomes and transcriptomes [15, 28]. Against this background, we here apply an RNAseq approach to identify pathways affected by adult temperature, which in turn may indicate the mechanisms causing phenotypic changes during acclimation, in the tropical butterfly Bicyclus anynana.

Since we are interested in adult acclimation here, all individuals were reared in a common environment. To subsequently induce acclimation responses, adult butterflies were exposed to different temperatures and also feeding treatments. Including the latter factor was motivated by the facts that (1) in nature animals are frequently faced with transient periods of food shortage negatively affecting fitness components (e.g. [29, 30]), and as (2) plasticity is assumed to involve costs such that inadequate nutrition may interfere with plastic responses [5]. According trade-offs though may be masked by having food access ad libitum. We deliberately decided to use non-stressful temperatures (19 and 27 °C) to induce acclamatory responses, as we were interested in the mechanisms conferring increased heat tolerance after acclimation rather than the stress response per se which is reasonably well understood (e.g. [26]). The temperatures chosen are well within the range of temperature fluctuations typically experienced by adult butterflies in their natural environment within their life spans. Additionally, we explore sexual differences in gene expression, as males and females typically differ in their responses to environmental variation (e.g. [31, 32]). Our study organism seems to be eminently suitable for our purpose as it inhabits a highly seasonal environment, thus relying heavily on phenotypic plasticity to master associated challenges [33]. Moreover, plastic responses to adult acclimation temperatures and feeding regimes have been described in detail. For instance, exposing adult butterflies for a few days to the temperatures used here results in an increase in heat tolerance by ca. 60%, a decrease in cold tolerance by ca. 36%, an increase in fecundity (ca. 22%) but a decrease in egg size (ca. 15%) at the higher compared with the lower temperature [4, 6, 8, 30, 34,35,36].

Our principle aim here is to identify pathways that are affected by adult temperature and which may therefore be mechanistically linked to the increased heat tolerance after acclimation to higher temperatures [10]. Specifically, we address the following hypotheses: (1) Temperature differences, although exclusively experienced in the adult stage, will affect expression profiles. While much of this variation will be attributable to changes in metabolism, we will specifically focus on pathways that may affect heat tolerance. Namely, an upregulation of the heat shock response and of oxidative defense mechanisms at the higher temperature is expected [37, 38]. (2) Owing to costs of plasticity, food stress will interfere with plastic responses as evidenced by the occurrence of temperature by feeding treatment interactions. Additionally, we explore sexual differences in the response to temperature, evidenced by the occurrence of temperature by sex interactions.


Variation in life-history and physiological traits

Our treatments were effective in manipulating phenotypic values. Adult temperature significantly affected thorax mass, abdomen mass, and thorax-abdomen ratio, but not total adult mass (see Table 1 for all statistical results). The lower compared with the higher adult temperature caused higher thorax (20.1 ± 0.3 mg > 19.6 ± 0.3 mg) but lower abdomen masses (37.4 ± 1.7 mg < 42.8 ± 1.9 mg), and higher thorax-abdomen ratios (0.76 ± 0.03 > 0.68 ± 0.03) owing to much larger temperature-induced variation in abdomen than in thorax mass (Fig. 1). All four traits measured differed significantly among sexes, indicating that males compared with females had lower adult body (43.8 ± 0.9 mg < 89.8 ± 1.1 mg), thorax (17.3 ± 0.2 mg < 21.9 ± 0.2 mg) and abdomen masses (17.8 mg ± 1.0 mg < 58.0 mg ± 1.0 mg) but a higher thorax-abdomen ratio (1.12 ± 0.02 > 0.40 ± 0.01). Sexual differences in body and abdomen mass were more pronounced at the higher than at the lower temperature (significant temperature by sex interactions). Adult feeding treatment only affected the thorax-abdomen ratio significantly, being higher in food-restricted than in control animals (0.72 ± 0.03 > 0.69 ± 0.03). Family effects were significant for thorax and abdomen mass, though variance component analyses revealed that they explained less than 0.01% of the variation.

Table 1 ANOVA results for variation in life history and physiological traits in Bicyclus anynana
Fig. 1

Trait variation in Bicyclus anynana butterflies across adult feeding treatments and two temperatures. Given are means + 1 SE for (a) total body mass, (b) abdomen mass, (c) thorax-abdomen ratio, (d) abdomen fat content, and (e) lysozyme activity for 8-day old adults. Food restriction = yes: only 30 min access to food per 48 h; food restriction = no: food access ad libitum throughout. Black bars: males; open bars: females. Group sample sizes range between 33 and 48

For physiological traits, adult temperature significantly affected fat content and lysozyme activity only, with fat content being higher at the higher temperature (20.1 ± 1.0% > 17.7 ± 0.8%) while lysozyme activity was higher at the lower temperature (0.37 ± 0.01 mOD/mL > 0.33 ± 0.01 mOD / mL; Table 1, Fig. 1). The only significant interaction was the one between temperature and sex for fat content, indicating that temperature effects were largely restricted to males (Fig. 1d). Differences between sexes were significant for all physiological traits except protein content. Females compared with males had higher lysozyme (0.39 ± 0.01 mOD/mL > 0.30 ± 0.01 mOD/mL) and ADH (0.023 ± 0.001 mOD/mL > 0.021 ± 0.001 mOD/mL) activities, but a lower relative fat content (11.6 ± 0.4% < 27.6 ± 0.9%). Feeding treatment significantly affected lysozyme activity only, showing that control individuals had a higher lysozyme activity than food restricted ones (0.37 ± 0.01 mOD/mL > 0.33 ± 0.01 mOD / mL). Additionally, there was a non-significant tendency towards a higher fat content in control individuals (19.5 ± 0.8% vs. 18.2 ± 0.9%). Family effects were significant for fat content and ADH activity only, but again explained only a minor fraction of the total variation (4.1 and 0.01%, respectively).

Gene expression

The 10,696 clusters of assembled transcripts mapped primarily to genes known from other insects (99.9% in total), especially to those from other Lepidoptera (93.4%; Additional file 1: Table S1). Searching the BUSCO database (version 9) for insect universal single-copy orthologs in our assembled transcriptome, identified 92.3% of families as at last once completely present, 51.6% as duplicated, 5.4% as fragmented, and 2.3% as missing. These results indicate that a large fraction of universally present insect genes are also represented in our assembled transcriptome. The large number of duplicated genes is expected as any gene is on average represented by multiple alternatively spliced or alternatively assembled transcripts. A transcriptome-wide overview of the influence temperature, sex, and feeding treatment shows that sex is the major factor influencing the expression profile of individuals (Fig. 2). In general, females appear to be less diverse than males with respect to expression profiles. Additionally, temperature explains some diversity in expression profiles among samples, albeit its effect is much smaller than that of sex. To further analyze the effect of temperature – possibly dependent on sex and feeding treatment – on the expression of individual genes, statistical tests were performed for each assembled transcript and results were associated with the transcript’s cluster as a proxy for the gene. We first identified transcripts that are differentially expressed between the two temperatures. This was achieved by testing the null hypothesis that temperature has neither an overall influence on expression, nor is interacting with sex or feeding treatment, nor jointly with sex and feeding treatment, i.e. using the script null hypothesis H0 : α3 = α5 = α7 =  α8 = 0. Accordingly, the expression of 3489 transcripts in total was significantly affected by temperature. When testing the hypotheses α3 = 0, α5 = 0, α7 = 0 and α8 = 0 individually, we obtained 1280 transcripts that were exclusively affected by temperature (α3 ≠ 0), 134 transcripts with a sex-dependent influence of temperature (α5 ≠ 0), 18 transcripts with a feeding-dependent influence of temperature (α7 ≠ 0), and 5 transcripts where the effect of temperature depended on both food and sex (α8 ≠ 0) (Additional file 2: Data S1). Thus, the vast majority of transcripts was influenced by temperature independently of sex and feeding treatment. We here consider only those transcripts, which could be assigned to a NCBI gi number, which explains the slightly reduced numbers in the Additional file 2: Data S1.

Fig. 2

Multidimensional Scaling plot of the expression profiles of all 75 individuals. The sexes are clearly divided and the expression profiles of females appear less diverse than those of males. Moreover, temperature (different colours) seems to have an impact, while an overall effect of feeding treatment (shape) is not apparent

The highest numbers of differentially expressed transcripts and clusters, respectively, were found to be affected by the factor sex, followed by temperature and the interaction between temperature and sex (Table 2; see also Additional file 2: Data S1). In contrast, the factor feeding regime and the remaining interactions caused differential gene expression in 5–29 transcripts/clusters only. In total, 13,570 transcripts and 6497 transcript clusters were significantly up- or down-regulated, showing annotation rates of 84.6% for the factor sex, 87.7% for temperature, and 88.5% for feeding regime.

Table 2 Number of differentially expressed transcripts and transcript clusters for different sources of variation

Out of the 1089 annotated (out of 1242) transcripts that were affected by the factor temperature alone, 659 were down- and 430 up-regulated at the higher temperature. The genes most strongly down-regulated at the higher temperature were a chaperonin complex component and a CUB domain-containing protein, and the genes most strongly up-regulated at the higher temperature were vacuolar H+ ATPase V1 sector, aldehyde dehydrogenase, and trypsin (Fig. 3). Genes associated with oxidative defense (peroxidase/oxygenase) were also among the most strongly upregulated ones. The transcripts being down-regulated at the higher temperature were mainly associated with cellular processes and signaling (31.7%) and information storage and processing (28.7%; Table 3a). 20.0% of down-regulated transcripts were only poorly characterized. Regarding sub-functions, the highest numbers of down-regulated transcripts were associated with posttranslational modification, protein turnover, chaperones (12.3%), translation, ribosomal structure and biogenesis (7.9%), transcription (7.7%), and RNA processing and modification (7.4%). Transcripts up-regulated at the higher temperature were mainly associated with cellular processes and signaling (42.1%) followed by information storage and processing (25.1%) and metabolism (15.8%; Table 3b). 17.0% of up-regulated transcripts were poorly characterized. The most frequently up-regulated sub-functions included signal transduction mechanisms (22.1%), transcription (12.3%), posttranslational modification, protein turnover, chaperones (6.1%), and RNA processing and modification (6.1%). An enrichment analysis indicated 24 enriched functional categories at the higher compared to the lower temperature. Enriched functions mainly included metabolic and biosynthetic pathways, but most importantly also responses to oxidative stress (nitric oxide biosynthetic process; Table 4).

Fig. 3

Overview of the genes being exclusively affected by the factor temperature. In these genes expression was independent of effects of sex and feeding treatment. Only genes with |logFC values| > 5 are shown. Red: up-regulation at the higher temperature; blue: down-regulation at the higher temperature. The darker the colour, the stronger the up- or down-regulation. Genes marked by * are private genes

Table 3 Functional annotation of transcripts (a) down- or (b) up-regulated at the higher temperature according to molecular main and sub-role
Table 4 Significantly enriched functional categories found to be a) down- or b) up-regulated at the higher temperature

Regarding interactive effects, 110 annotated transcripts were significantly affected by temperature and sex. Of these, 49 were significantly more strongly down-regulated and 61 significantly more strongly up-regulated in females than in males at the higher temperature (Additional file 2: Data S1). The specific genes showing the largest sexual differences in their response to temperature were the RAS effector and lipid phosphate phosphatase (both down-regulated), and the allantoicase and CUB domain-containing protein (up-regulated). For combined effects of temperature and feeding treatment we found 14 annotated transcripts of which 9 were significantly more strongly down-regulated and 5 significantly more strongly up-regulated at the higher temperature under control conditions versus food restriction (Additional file 2: Data S1). The RNA-binding protein fusilli, trypsin, and molecular chaperones were more strongly down-regulated, whereas an extracellular matrix glycoprotein and a calcium transporting ATPase were more strongly up-regulated under control conditions versus food restriction. The interaction between sex and feeding regime revealed 12 annotated transcripts. Five of these were significantly more strongly down-regulated in females compared with males when being fed ad libitum, including F-box proteins and prohibitins and stomatins, while 7 were significantly more strongly up-regulated. Finally, 5 annotated transcripts were affected by the 3-way interaction. Three were down-regulated (thioredoxin binding protein, trypsin, myosin) and 2 up-regulated (RAS effector RIN1, transcription factor zerknullt) at the higher temperature in females compared to males when being fed ad libitum.

Regarding sex differences in gene expression, out of 10,265 annotated differentially expressed transcripts 4853 were down- and 5412 up-regulated in females relative to males. The specific genes most strongly down-regulated in females were trypsin and a nuclear pore complex component, and the genes most strongly up-regulated in females were a lipoprotein and DNA polymerase alpha (Additional file 3: Figure S1). The transcripts being down-regulated in females were mainly related to metabolism (36.7%) and cellular processes and signaling (32.1%), while those being up-regulated were mainly associated with information storage and processing (36.3%) and cellular processes and signaling (32.2%; Additional file 3: Figure S2). We found 85 enriched functional categories in females compared with males, the majority of which was associated with metabolism (Additional file 4: Data S2). Furthermore, up-regulated enriched functional categories in females included stress responses such as DNA repair, heme oxidation, and protein refolding. See below for potential effects of dosage compensation on sexual differences in expression patterns.

Only 26 annotated transcripts were differentially expressed among feeding treatments, of which 14 were down-regulated and 12 up-regulated under control conditions relative to food restriction. Genes most strongly up-regulated under control conditions were phosphatidylinositol transfer protein SEC14 and prohibitins and stomatins, and genes most strongly down-regulated were pyruvate kinase and the thiamine pyrophosphate-requiring enzyme (Additional file 3: Figure S1). Seven of the down-regulated transcripts were related to cellular processes and signaling, three to metabolism, and two to information storage and processing. We found five up-regulated transcripts related to metabolism, three to cellular processes and signaling, and two to information storage and processing (Additional file 3: Figure S2). We found four enriched functional categories affected by feeding treatment (Additional file 4: Data S2).


Variation in life history and physiological traits

Our data show that the adult treatments employed were effective in changing phenotypic traits, as evidenced by significant effects on thorax mass, abdomen mass, thorax-abdomen ratio, fat content, and lysozyme activity. Under cooler conditions, butterflies had heavier thoraces but lighter abdomen and consequently a higher thorax-abdomen ratio, and food-restricted animals also had a higher thorax-abdomen ratio. These differences probably reflect variation in energy expenditure and acquisition. Higher temperatures and food access enable higher levels of food acquisition, shifting thorax-abdomen ratios towards lower values. This interpretation is in line with the results on fat content, being or tending to be higher at the warmer temperature and unrestricted access to food [6, 39]. Lysozyme activity was reduced at the higher temperature and under food restriction, corroborating earlier findings indicating detrimental effects of such conditions on immune competence [40, 41]. The temperatures used here also induce substantial variation in other traits in B. anynana, including cold and heat stress resistance, egg size, egg composition, fecundity, and reproductive output [6, 8, 36, 42,43,44].

As expected, B. anynana showed also pronounced sexual differences. Higher female masses are frequently found in insects, which is likely caused by fecundity selection [39, 41, 45]. The males’ higher relative fat content probably reflects their higher energy demand to sustain high flight activity e.g. during mate location [6, 46]. Concomitantly, males also show a higher thorax-abdomen ratio [47, 48]. Regarding physiological traits, the higher lysozyme activity in females may reflect an increased investment into immune function and thereby longevity to maximize reproductive output [49], while their higher ADH activity (cf. [34] may reflect higher metabolic rates which may in turn increase egg production [50]. Three sex by temperature interactions were significant, indicating that females were able to benefit more strongly from the higher temperature and thus exemplifying sex-specific responses to acclimation temperatures. Other examples of sex-specific temperature effects in butterflies include heat shock protein expression and longevity [39, 41].

Effects of acclimation temperature on gene expression

The vast majority of assembled transcript clusters mapped to genes of other insects and especially other Lepidoptera, indicating that we managed to de novo assemble a high-quality transcriptome of the butterfly B. anynana. Although we have deliberately chosen to focus on short-term adult acclimation, acclimation temperature had a significant impact on gene expression (hypothesis 1). The non-stressful adult temperatures used here do cause striking variation in phenotype, and based on the expression patterns found we conclude that those shifts are induced by changes in a large number of pathways.

Overall, fewer genes were up- than down-regulated at the higher temperature, the former including pathways associated with signal transduction, transcription, posttranslational modification, protein turnover, chaperones, RNA processing and modification, cytoskeleton, and lipid metabolism. For all these processes also down-regulation at the higher temperature was found. Thus in contrast to hypothesis 1, we could not find strong evidence for an overall upregulation of genes related to metabolism at the higher temperature. Nevertheless, several genes with particular high logFC values were related to metabolism, such as ATPase, aldehyde dehydrogenase, trypsin, S-adenosylmethionine synthetase, and beta-glucosidase. These findings may reflect the fact that, in ectotherms, metabolism and enzymatic activity generally increase with increasing temperature [51, 52]. Metabolic activity may also increase at (stressfully) low temperatures [53] or may remain unaffected by temperature variation [54, 55], which may explain the rather high level of down-regulation at the higher temperature found here.

Our principal aim though was to identify candidate pathways that may confer higher heat tolerance after exposure to non-stressful warm temperatures. Interestingly, we indeed found a strong upregulation of genes related to the oxidative stress response at the higher temperature as predicted, including for instance peroxidase/oxygenase (Fig. 3), the Cytochrome P450 superfamily (Additional file 2: Data S1; cf. [56]), heme binding, and oxidoreductase activity (Additional file 2: Data S1). Also, the functional category ‘nitric oxide biosynthetic process’ was overrepresented at the higher temperature. Higher temperatures elevate the production of reactive oxygen species as a by-product of the aerobic metabolism, which may cause oxidative damage [57, 58]. A concomitant up-regulation of antioxidant enzymes may in turn alleviate the resulting oxidative stress, and has therefore been implied to play a crucial role for heat stress resistance in insects [59,60,61,62]. Antioxidant enzymes, in particular superoxide dismutase, have also been implied to play an important role in longevity and starvation resistance in B. anynana and other animals [63,64,65]. Furthermore, responses to oxidative stress may mediate the trade-off between reproduction and longevity, at least under stressful thermal conditions [66]. Thus, upregulation of antioxidant enzymes comprises a straightforward candidate mechanism that may underlie increased heat tolerance after acclimation to warmer temperatures, helping to safeguard homeostasis if subsequently exposed to stressful conditions.

In contrast, we found no support for a meaningful involvement of the heat shock response in acclimating to warmer temperatures. The heat shock response is basically an emergency reaction under acute stress [37]. Hence, we assume that the non-stressful acclimation temperatures used here were not high enough to induce the heat stress response (cf. [67]). Interestingly, similar temperatures did up-regulate heat-shock proteins in a temperate-zone Copper butterfly [68]. This suggests that the tropical butterfly investigated here, in contrast to the aforementioned temperature-zone species, did not experience 27 °C as being stressful and that consequently higher temperatures are needed for up-regulation [53, 67, 69].

Effects of acclimation temperature in interaction with other factors on gene expression

At least some transcripts were involved in the interaction between feeding treatment and temperature, in line with our second hypothesis. Thus, we indeed found evidence that inadequate nutrition may interfere with plastic responses, as evidenced for instance by a stronger up-regulation of molecular chaperones at the lower temperature in the feeding as compared with the food-restricted group. Also the fusilli protein, which is important in embryogenesis and other developmental processes [70], was significantly more strongly down-regulated at 27 °C in the feeding as compared with the food-restricted group. This may in turn impact food intake, ensuring growth and eventually maturation [10].

As expected, sexes responded in part differently to temperature, evidenced by significant sex by temperature interactions for expression patterns. This interaction affected the expression of 111 annotated genes. Examples include genes associated with energy production and lipid transport (lipid phosphate phosphatase, Di-hydrolipoamide acetyltransferase). These were more strongly down-regulated in females than in males at the higher temperature. While sex specific responses to temperature clearly do occur, many traits including temperature stress resistance show similar patterns in both sexes (e.g. [10, 40, 41, 67]), which may explain the relatively low incidence of interactive effects.

Only five transcripts were affected by the 3-way interaction. We found a strong down-regulation of the thioredoxin binding protein, which plays an important role in redox homeostasis by increasing reactive oxygen species and oxidative stress [71], at the higher temperature in females compared with males when being fed ad libitum. This suggests that females, under normal feeding conditions at 27 °C, were less exposed to oxidative stress than the other groups. Furthermore, the RAS effector RIN1 was up-regulated at the higher temperature in females compared with males when being fed ad libitum, which plays a pivotal role in activating various cellular processes [72]. Thus, females activate cellular processes under favourable conditions.

Effects of feeding treatment and sex on gene expression

The fact that only a small number of transcripts was differentially expressed among feeding treatments suggests that butterflies could largely compensate for the temporary lack of food (see also [73]). However, at least some functions seem to have been compromised by limited food access, evidenced by changes in thorax-abdomen ratio, lysozyme activity, and fat content (see above). Accordingly, at least a few metabolic functions such as carbohydrate and amino acid transport and metabolism were down-regulated under food restriction.

We found pronounced sexual differences in expression patterns as expected, being at least partly related to reproduction. For example, the genes most strongly up-regulated in females as opposed to males were lipoproteins, being important for egg production [74, 75]. Furthermore, metabolism (e.g. trypsin, enolase, fatty acid desaturase) and innate immunity (e.g. serpin) appeared to be down-regulated in females compared with males. Whether the overall down-regulation of genes related to innate immunity and metabolic pathways indicates a trade-off between investment into maintenance versus reproduction (cf. [66, 76]), requires further investigation. To what extant sex-specific differences in gene expression are affected by (a lack of) dosage compensation in butterflies, in which females are the heterogametic sex, is currently unknown [77,78,79,80]. Only very few genes were differentially expressed between the sexes in interaction with feeding regime, suggesting that food stress had similar effects in both sexes.


Based on a high-quality transcriptome, our study revealed that, as expected, sex had overall the largest impact on expression patterns. At the higher temperature, several antioxidant markers were up-regulated, which may play an important role in the acclamatory processes subsequently providing increased heat tolerance. In contrast, the heat shock response did not seem to be involved in the acclamatory response. We thus believe that the role of anti-oxidative responses in increasing heat tolerance may have been hitherto underestimated. Their upregulation may provide a crucial mechanism to prepare organisms for later more stressful conditions. Such conditioning may arise from higher basal levels of anti-oxidant enzymes or a quicker response to deteriorating conditions as the respective ‘machinery’ has been already turned on. While the important role of antioxidant enzymes in alleviating detrimental effects of oxidative stress is firmly established [58, 81], we believe that its potentially crucial role for mediating increased stress tolerance, even under non-stressful conditions, deserves more attention (e.g. [59, 62, 65, 82, 83]). Additionally, we provide molecular evidence that plastic responses may be compromised by inadequate nutrition. Our study though shows that even a few days at different temperatures in the adult stage affect expression patterns and phenotypes including thermal stress resistance [10]. Under natural conditions, the phenotypic changes induced (namely increased heat tolerance) are likely to sustain longer lifespans, higher levels of activity, and a higher reproductive output under stress (cf. [12, 13]).


Study organism

Bicyclus anynana, a tropical fruit-feeding butterfly, is distributed from southern Africa to Ethiopia [84]. As an adaptation to alternate wet-dry seasonal environments and the associated changes in resting background and predation, this species exhibits striking phenotypic plasticity (two seasonal morphs; [7]). Reproduction is confined to the favorable wet season during which oviposition plants are abundantly available [35, 85]. In 1988, a laboratory stock population was established at Leiden University, the Netherlands, from over 80 gravid females collected at a single locality in Malawi. Several hundred adults are reared in each generation, maintaining high levels of heterozygosity at neutral loci [86]. From the Leiden stock population, a laboratory population was established at Greifswald University, Germany, in 2007 from which all animals for this experiment originated.

Experimental design

To start this experiment, freshly eclosed virgin males and females (50 each) from the stock population were mated. Single-mated females were afterwards kept individually in 1 L translucent plastic containers for oviposition to create full-sib families. They were provided with cuttings of maize as oviposition substrate and moist banana for adult feeding on a daily basis. All females and their offspring were kept in a common environment inducing wet season phenotypes, i.e. at 23 °C, 70% relative humidity, and a photoperiod of L12:D12. The eggs laid by individual females were removed from the containers and transferred to petri-dishes. After hatching, larvae were transferred individually to small plastic pots (125 ml) lined with moist paper tissue and were fed on cuttings of young maize plants being replaced every other day. Following adult eclosion, butterflies were randomly allocated to one out of four treatment groups, differing in temperature (19 °C or 27 °C) and feeding treatment, using a split-family design (n = 12 families, n = 494 individuals). Butterflies were, based on pilot experiments, either fed with moist banana ad libitum (control) or had access to banana for only 30 min/48 h, resulting in four treatment groups: (1) 19 °C control (n = 122), (2) 19 °C food-restricted (n = 127), (3) 27 °C control (n = 123), (4) 27 °C food-restricted (n = 122). The above temperatures were chosen due to their ecological relevance and because they were known to induce large differences in thermal tolerance and other traits within a few days [8, 10]. On day 8 of adult life, all butterflies were shock-frozen in liquid nitrogen and stored at − 80 °C for later analyses. Food-restricted butterflies were not fed for ca. 44 h before being frozen to maximize effects on phenotypes and expression patterns.

Life-history traits and physiological measurements

As proxies of condition, we scored adult body mass, thorax mass, abdomen mass, thorax-abdomen ratio, fat content, lysozyme activity, alcohol dehydrogenase (ADH) activity, and protein content for all individuals. Lysozyme activity was included as an indicator of investment into immune function [49], while ADH activity may reflect metabolic rates [34, 50]. For subsequent measurements, head, legs, and wings were removed, and abdomen and thorax separated. All masses were weighed as fresh mass to the nearest 0.01 mg (Sartorius microscale LE225D, Sartorius AG, Goettingen, Germany). Butterfly abdomens were used for measuring relative fat content, except for individuals used for RNA extractions (see below). Fat content was determined as difference in abdomen dry mass before and after two fat extractions in percent. Therefore, abdomens were dried to constant mass at 60 °C for 48 h. For fat extractions, 0.4 ml of dichloromethane (CH2Cl2) per abdomen was applied for 24 h (cf. [87]).

For the determination of lysozyme activity, ADH activity, and protein content, the thoraces were used. Each thorax was homogenized in 100 μl phosphate-buffered saline (PBS; 11.9 mM Na2HPO4 * 2 H2O, 137 mM NaCl, 2.7 mM KCl, pH 8.0) and centrifuged at 4 °C for 20 min (14,000 rpm). Lysozyme activity was determined following [88]. The wells of a 96-well plate were loaded with 20 μl supernatant plus 80 μl Micrococcus luteus solution (3 mg/ml in PBS) or blanks containing 20 μl PBS plus 80 μl M. luteus solution. The optical density was measured at 490 nm and 30 °C for 5 h (microplate reader BioTekELx 808, Bad Friedrichshall, Germany). For the determination of ADH activity, wells were loaded with 10 μl supernatant plus 190 μl reaction solution (30 mM isopropanol [23 μl/10 ml TRIS HCl], 3 mM NAD+ [0.0199 g/10 ml TRIS HCl], 0.15 M TRIS HCl buffer pH 8.5) or blank samples (with 10 μl PBS plus 190 μl reaction solution). ADH activity was measured as change in the optical density at 340 nm and 30 °C for 10 min (15 s steps; BioTekELx 808). Both activity measures were calculated by subtracting the final from the initial value, and by additionally correcting for the mean of the blank values. Total protein content was quantified using the BioRad protein assay based on the Bradford method [89]. Therefore, 1 μl of the supernatant was diluted in 160 μl aqua dest. After adding 40 μl BioRad solution and 10 min of incubation, the absorbance was read at 595 nm and 30 °C (BioTekELx 808). Four replicates were measured per individual. A standard curve was constructed with albumine bovine serum in cacodylate buffer, using a concentration series (0–2 mg/ml).

Data on life-history and physiological traits were analyzed using analyses of variance (ANOVAs) with adult temperature, adult feeding treatment, and sex as fixed effects and family as random factor. To test for normality and variance homogeneity, we used the Kolmogorov-Smirnov and the Levene test, respectively. If necessary, data were transformed as appropriate. Pair-wise comparisons were performed employing Tukey’s HSD. Minimum adequate models were constructed by sequentially removing non-significant interaction terms. The above statistical tests were performed using Statistica (8.0). Throughout the text, estimates are stated as mean ± standard error.

RNA isolation and purification

Whole abdomens were used for RNA extraction (n = 80), as the thoraces were used for the above analyses. We selected 10 butterfly families (i.e. the offspring of 10 single-mated females). RNA was extracted from one individual per family and treatment group (i.e. per 2 temperatures * 2 feeding treatments * 2 sexes, resulting in 8 individuals per family). Total sample size was 77 instead of 80 individuals (8 * 10 families), as RNA quality was insufficient in 3 individuals from 27 °C. RNA isolation was carried out using TRIZOL (Invitrogen) according to the manufacturer’s instructions, using one ml of TRIZOL per abdomen. For RNA purification, the RNeasy Mini kit (Qiagen 74,106) and the RNAase free DNAase kit (Qiagen 792,454) were used. The RNA was eluted in 100 μl and stored at − 80 °C for later analysis (for a detailed description see Additional file 5).

Sequencing, transcriptome assembly, and expression analysis

After RNA extraction, mRNA of eukaryotes was enriched by using oligo(dT) magnetic beads. By adding fragmentation buffer, the mRNA was interrupted to short fragments (about 200 bp). Then, the first strand of cDNA was synthesized by a random hexamer-primer using the mRNA fragments as templates. Buffer, dNTPs, RNase, and DNA polymerase were added to synthesize the second strand. The double strand cDNA was purified with the QiaQuick PCR extraction kit and washed with EB buffer for end repair and single nucleotide A (adenine) addition. Finally, sequencing adaptors were ligated to the fragments. The required fragments were purified by agarose gel electrophoresis and enriched by PCR amplification.

Each sample was sequenced as a separate paired end library on a HiSeq2000 Illumina platform by BGI Hong Kong Co. Ltd., with 9 to 10 samples per lane. In total we obtained 2,066,491,822 unstranded paired-end reads (100 bp read length), with 8–14 million reads per individual. As its usefulness is under debate, no filtering and trimming was applied (except for adaptors) [90]. Additionally, using quality trimming in the Trinity assembler resulted in only 0.5% fewer contigs, and was therefore deemed unnecessary. Given that a reference genome of B. anynana is not available, we performed a de novo assembly. We used and compared two popular software packages, SOAP denovo-trans (version 1.02; [22]) and Trinity (version r2013-02-25; [91]). We used two quality measures to choose between both assemblies and to reduce the number of contigs to meaningful transcripts. First, we aligned the transcripts against the proteins of the model organism Drosophila melanogaster as quality criterion (Additional file 6: Table S2). The analyses indicated that using SOAP 1574 proteins aligned with at least 80% coverage, whereas in the case of Trinity 3694 proteins did so. Second, the most likely coding regions were extracted from both assemblies. Accordingly, a total of 25,016 SOAP and 66,712 Trinity transcripts that correspond to likely coding regions were obtained. Based on both quality measures Trinity produced a better assembly, such that the SOAP assembly was discarded.

For the Trinity assembly, transcript abundances were computed with eXpress (version 1.3.1; Reads were mapped to the transcript assembly with Bowtie 2 [92]. Default values of eXpress and Bowtie 2 were used except for parameter B (additional number of EM iterations) which was set to 20. Abundance levels are reported as reads per kilobase of transcript per million mapped reads (RPKM; [93]). To identify differential expression of assembled transcripts, eXpress [94] was used to estimate the relative expression count of each assembled transcript for each individual. ‘Effective counts’, comprising the value recommended for input to count-based differential expression tools (, were rounded to the nearest nonnegative integer. For 30 out of the 66,712 assembled transcripts this count could not be computed by eXpress such that they were subsequently ignored. The resulting 66,682 × 77 matrix was statistically analyzed with edgeR 3.16.5 [95].

First, samples were visualized in two dimensions using multidimensional scaling (function plotMDS of edgeR), which revealed that two samples (one male, one female) were likely wrongly labeled (Additional file 3: Figure S3). These two samples were subsequently removed, leaving 75 samples. Very lowly expressed transcripts were not considered in the differential expression analysis: we ignored transcripts that did not have more than 50 reads in at least 5 samples. Libraries were normalized for statistical comparison using the standard trimmed mean of M-values of edgeR (calcNormFactors). The count data was assumed to follow a negative binomial distribution. A full-factorial generalized linear model was used to model the expression counts, with the binary factors temperature t (19 °C or 27 °C), sex s (male or female), and feeding treatment f (ad libitum or food stress) including all possible interactions.

$$ \ln {\mu}_{g,j}={\alpha}_1+{\alpha}_2{s}_j+{\alpha}_3{t}_j+{\alpha}_4{f}_j+{\alpha}_5{s}_j{t}_j+{\alpha}_6{s}_j{f}_j+{\alpha}_7{t}_j{f}_j+{\alpha}_8{s}_j{t}_j{f}_j $$

Here, μg, j is the mean expression of transcript g in sample j and the variables sj, tj and fj assume either the value 0 or 1. For an alternative approach with separate analyses for each sex see Additional file 7. P-values were adjusted for multiple testing with the Benjamini-Hochberg method [96]. Afterwards, we explored the effects of temperature on gene expression, which was the principle aim of this study, in more detail. Therefore, we compiled for each null hypothesis outlined below (see Results) a list of significant transcripts. Significance was determined using a false discovery rate (FDR) threshold of 5%.

The 66,712 assembled transcripts were aligned against the hexapoda section of the non-redundant protein database (February 2015) using blastx version 2.2.26 [97]. 96% of the assembled transcripts (64,047) had a significant hit at an e-value threshold of 10− 5, indicating that a large fraction of assembled transcripts contained protein-coding sequences. We further applied BUSCO version 3 [98] to assess the completeness of the set of assembled transcripts that were likely to contain coding sequences (see Results). As several transcripts may map to the same gene (e.g. alternative splicing, incomplete assembly), they were assigned to clusters when having a hit to the same target protein (single-linkage clustering). For clustering, only hits with (1) an e-value below 10− 5 that (2) belong to the top ten hits of a transcript and for which (3) the e-value is not larger than a thousand times the smallest e-value for that transcript were considered. These filters were introduced to prevent too broad clusters resulting from domains that are shared between many proteins. Each of the resulting 10,696 clusters thus likely represents assembled transcripts from the same gene or group of homologs.

The assembled transcript clusters were annotated against cluster of orthologous groups (COG) and eukaryotic orthologous groups (KOG; RPS-BLAST algorithms) using the program Prophane [99] and NCBI gi numbers (Additional file 2: Data S1). The COG/KOG system distinguishes among three protein main roles and 23 sub-roles. Additionally, Blast2GO ( [100]) was used to assign transcript clusters to gene ontology (GO) terms (; Additional file 8: Data S3). As both annotations gave very similar results, we only present the results based on the COG/KOG system here (for results based on the GO system see Additional file 9). We finally performed enrichment analyses based on GO terms. GO IDs were retrieved using Interproscan v5.27–66.0. For the detection of significantly enriched GO terms, we used a Fishers exact test as implemented in the R package TopGo [101]. The complete contig set was taken as reference, with the sets of treatment-specific expressed genes comprising the test sets.


  1. 1.

    Walther GR, Post E, Convey P, Menzel A, Parmesan C, Beebee TJC, et al. Ecological responses to recent climate change. Nature. 2002;416:389–95.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Dahlhoff EP, Rank NE. The role of stress proteins in responses of a montane willow leaf beetle to environmental temperature variation. J Biosci. 2007;32:477–88.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  3. 3.

    Pigliucci M. Phenotypic plasticity: beyond nature and nurture. Baltimore: John Hopkins University Press; 2001.

    Google Scholar 

  4. 4.

    Fischer K, Karl I. Exploring plastic and genetic responses to temperature variation using copper butterflies. Clim Res. 2010;43:17–30.

    Article  Google Scholar 

  5. 5.

    Whitman DW, Agrawal AA. What is phenotypic plasticity and why is it important? In: Whitman DW, Ananthakrishnan TN, editors. Phenotypic plasticity of insects. Enfield, NH: Science Publishers; 2008. p. 1–63.

    Google Scholar 

  6. 6.

    Fischer K, Eenhoorn E, Bot ANM, Brakefield PM, Zwaan BJ. Cooler butterflies lay larger eggs: developmental plasticity versus acclimation. Proc R Soc Lond B. 2003;270:2051–6.

    Article  Google Scholar 

  7. 7.

    Lyytinen A, Brakefield PM, Lindström L, Mappes J. Does predation maintain eyespot plasticity in Bicyclus anynana? Proc R Soc Lond B. 2004;271:279–83.

    Article  Google Scholar 

  8. 8.

    Geister TLK, Fischer K. Testing the beneficial acclimation hypothesis: temperature effects on mating success in a butterfly. Behav Ecol. 2007;18:658–64.

    Article  Google Scholar 

  9. 9.

    Karl I, Janowitz SA, Fischer K. Altitudinal life-history variation and thermal adaptation in the copper butterfly Lycaena tityrus. Oikos. 2008;117:778–88.

    Article  Google Scholar 

  10. 10.

    Fischer K, Dierks A, Franke K, Geister TL, Liszka M, Winter S, et al. Environmental effects on temperature stress resistance in the tropical butterfly Bicyclus anynana. PLoS One. 2010;5:e15284.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. 11.

    Luo S, Wong SC, Xu C, Hanski I, Wang R, Lehtonen R. Phenotypic plasticity in thermal tolerance in the Glanville fritillary butterfly. J Therm Biol. 2014;42:33–9.

    CAS  PubMed  Article  Google Scholar 

  12. 12.

    Kristensen TN, Loeschcke V, Hoffmann AA. Can artificially selected phenotypes influence a component of field fitness? Thermal selection and fly performance under thermal extremes. Proc R Soc Lond B. 2007;274:771–8.

    Article  Google Scholar 

  13. 13.

    Overgaard J, Sørensen JG, Jensen LT, Loeschcke V, Kristensen TN. Field tests reveal genetic variation for performance at low temperatures in Drosophila melanogaster. Funct Ecol. 2010;24:186–95.

    Article  Google Scholar 

  14. 14.

    Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM. Beyond predictions: biodiversity conservation in a changing climate. Science. 2011;332:53–8.

    CAS  PubMed  Article  Google Scholar 

  15. 15.

    Zhou S, Campell TG, Stone EA, Mackay TFC, Anholt RRH. Phenotypic plasticity of the Drosophila transcriptome. PLoS Genet. 2012;8:e1002593.

    PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Piersma T, Drent J. Phenotypic flexibility and the evolution of organismal design. Trends Ecol Evol. 2003;18:228–33.

    Article  Google Scholar 

  17. 17.

    Terblanche JS, Chown SL. The relative contributions of developmental plasticity and adult acclimation to physiological variation in the tsetse fly, Glossina pallidipes (Diptera, Glossinidae). J Exp Biol. 2006;209:1064–73.

    PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Bowler K. Acclimation, heat shock and hardening. J Therm Biol. 2005;30:125–30.

    Article  Google Scholar 

  19. 19.

    Gerken AR, Eller OC, Hahn DA, Morgan TJ. Constraints, independence, and evolution of thermal plasticity: probing genetic architecture of long- and short-term thermal acclimation. Proc Natl Acad Sci U S A. 2015;112:4399–404.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Huey RB, Berrigan D, Gilchrist GW, Herron JC. Testing the adaptive significance of acclimation: a strong inference approach. Am Zool. 1999;39:323–37.

    Article  Google Scholar 

  21. 21.

    Foray V, Desouhant E, Voituron Y, Larvor V, Renault D, Colinet H, et al. Does cold tolerance plasticity correlate with the thermal environment and metabolic profiles of a parasitoid wasp? Comp Biochem Physiol A. 2013;164:77–83.

    CAS  Article  Google Scholar 

  22. 22.

    Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. GigaScience. 2012;1:18.

    PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Franssen SU, Gu J, Bergmann N, Winters G, Klostermeier UC, Rosenstiel P, et al. Transcriptomic resilience to global warming in the seagrass Zostera marina, a marine foundation species. Proc Natl Acad Sci U S A. 2011;108:19276–81.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  24. 24.

    MacMillan HA, Sinclair BJ. Mechanisms underlying insect chill-coma. J Insect Physiol. 2011;57:12–20.

    CAS  PubMed  Article  Google Scholar 

  25. 25.

    Geerts AN, Vanoverbeke J, Vanschoenwinkel B, Van Doorslaer W, Feuchtmayr H, Atkinson D, Moss B, Davidson TA, Sayer CD, De Meester L. Rapid evolution of thermal tolerance in the water flea Daphnia. Nat Climate Change. 2015;5:665–8.

    Article  Google Scholar 

  26. 26.

    Liu Y, Su H, Li R, Li X, Xu Y, Dai X, Zhou Y, Wang H. Comparative transcriptome analysis of Glyphodes pyloalis Walker (Lepidoptera: Pyralidae) reveals novel insights into heat stress tolerance in insects. BMC Genomics. 2017;18:974.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  27. 27.

    de Jong MA, Saastamoinen M. Environmental and genetic control of cold tolerance in the Glanville fritillary butterfly. J Evol Biol. 2018;31:636–45.

    PubMed  PubMed Central  Article  Google Scholar 

  28. 28.

    Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2010;12:87–98.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  29. 29.

    Metcalfe NB, Monaghan P. Compensation for a bad start: grow now, pay later? Trends Ecol Evol. 2001;16:254–60.

    PubMed  Article  Google Scholar 

  30. 30.

    Bauerfeind SS, Fischer K. Effects of food stress and density in different life stages on reproduction in a butterfly. Oikos. 2005;111:514–24.

    Article  Google Scholar 

  31. 31.

    Fischer K, Fiedler K. Sex-related differences in reaction norms in the butterfly Lycaena tityrus (Lepidoptera: Lycaenidae). Oikos. 2000;90:372–80.

    Article  Google Scholar 

  32. 32.

    Fischer K, Fiedler K. Dimorphic growth patterns and sex-specific reaction norms in the butterfly Lycaena hippothoe sumadiensis. J Evol Biol. 2001;14:210–8.

    Article  Google Scholar 

  33. 33.

    Brakefield PM, Pijpe J, Zwaan BJ. Developmental plasticity and acclimation both contribute to adaptive responses to alternating seasons of plenty and of stress in Bicyclus butterflies. J Biosci. 2007;3:465–75.

    Article  Google Scholar 

  34. 34.

    Fischer K, Klockmann M, Reim E. Strong negative effects of simulated heat waves in a tropical butterfly. J Exp Biol. 2014;217:2892–8.

    PubMed  Article  Google Scholar 

  35. 35.

    Brakefield PM. Phenotypic plasticity and fluctuating asymmetry as response to environmental stress in the butterfly Bicyclus anynana. In: Bijlsma RR, Loeschcke V, editors. Environmental stress: adaptation and evolution. Basel: Birkhäuser; 1997. p. 65–78.

    Chapter  Google Scholar 

  36. 36.

    Fischer K, Perlick J, Galetz T. Residual reproductive value and male mating success: older males do better. Proc R Soc Lond B. 2008;275:1517–24.

    Article  Google Scholar 

  37. 37.

    Sørensen JG, Kristensen TN, Loeschcke V. The evolutionary and ecological role of heat shock proteins. Ecol Lett. 2003;6:1025–37.

    Article  Google Scholar 

  38. 38.

    Lalouette L, Williams CM, Hervant F, Sinclair BJ, Renault D. Metabolic rate and oxidative stress in insects exposed to low temperature thermal fluctuations. Comp Biochem Physiol A. 2011;158:229–34.

    CAS  Article  Google Scholar 

  39. 39.

    Karl I, Fischer K. Why get big in the cold? Towards a solution of a life-history puzzle. Oecologia. 2008;155:215–25.

    PubMed  Article  Google Scholar 

  40. 40.

    Karl I, Stoks R, De Block M, Janowitz SA, Fischer K. Temperature extremes and butterfly fitness: conflicting evidence from life history and immune function. Glob Change Biol. 2011;17:676–87.

    Article  Google Scholar 

  41. 41.

    Franke K, Fischer K. Effects of inbreeding and temperature stress on life history and immune function in a butterfly. J Evol Biol. 2013;26:517–28.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Steigenga MJ, Fischer K. Ovarian dynamics, egg size and egg number in relation to temperature and mating status in a butterfly. Entomol Exp Appl. 2007;125:195–203.

    Article  Google Scholar 

  43. 43.

    Geister TL, Lorenz MW, Hoffmann KH, Fischer K. Adult nutrition and butterfly fitness: effects of diet quality on reproductive output, egg composition, and egg hatchling success. Front Zool. 2008;5:10.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  44. 44.

    Geister TL, Lorenz MW, Hoffmann KH, Fischer K. Energetic of embryonic development: effect of temperature on egg and hatchling composition in a butterfly. J Comp Physiol B. 2009;179:87–98.

    PubMed  Article  Google Scholar 

  45. 45.

    Honěk A. Intraspecific variation in body size and fecundity in insects: a general relationship. Oikos. 1993;66:483–92.

    Article  Google Scholar 

  46. 46.

    Roff DA. Life history evolution. Sunderland: Sinauer Associates Inc; 2010.

  47. 47.

    Berwaerts K, Van Dyck H, Aerts P. Does flight morphology relate to flight performance? An experimental test with the butterfly Pararge aegeria. Funct Ecol. 2002;16:484–91.

    Article  Google Scholar 

  48. 48.

    Van Dyck H, Wiklund C. Seasonal butterfly design: morphological plasticity among three developmental pathways relative to sex, flight and thermoregulation. J Evol Biol. 2002;15:216–25.

    Article  Google Scholar 

  49. 49.

    Schwartz A, Koella JC. The cost of immunity in the yellow fever mosquito, Aedes aegypti depends on immune activation. J Evol Biol. 2004;17:834–40.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Cosmides N, Loukas M, Zouros E. Differences in fitness components among alcohol dehydrogenase genotypes of the olive fruit fly (Diptera: Tephritidae) under artificial rearing. Ann Entomol Soc Am. 1997;90:363–71.

    CAS  Article  Google Scholar 

  51. 51.

    Gillooly JF, Brown JH, West GB, Savage VM, Charnow EL. Effects of size and temperature on metabolic rate. Science. 2001;293:2248–51.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Clarke A, Fraser KPP. Why does metabolism scale with temperature? Funct Ecol. 2004;18:243–51.

    Article  Google Scholar 

  53. 53.

    Xiao R, Wang L, Cao Y, Zhang G. Transcriptome response to temperature stress in the wolf spider Pardosa pseudoannulata (Araneae: Lycosidae). Ecol Evol. 2016;

  54. 54.

    Li S, Huang J, Liu C, Liu Y, Zheng G, Xi L, et al. Interactive effect of seawater acidification and elevated temperature on the transcriptome and biomineralization on the pearl oyster Pinctada fucata. Environ Sci Technol. 2016;50:1157–65.

    CAS  PubMed  Article  Google Scholar 

  55. 55.

    Levine MT, Eckert ML, Begun DJ. Whole-genome expression plasticity across tropical and temperate Drosophila melanogaster populations from eastern Australia. Mol Biol Evol. 2011;28:249–56.

    CAS  PubMed  Article  Google Scholar 

  56. 56.

    Torson AS, Yocum GD, Rinehart JP, Kemp WP, Bowsher JH. Transcriptional responses to fluctuating thermal regimes underpinning differences in survival in the solitary bee Megachile rotundata. J Exp Biol. 2015;218:1060–8.

    PubMed  Article  Google Scholar 

  57. 57.

    Tumminello RA, Fuller-Espie SL. Heat stress induces ROS production and histone phosphorylation in celomocytes of Eisenia hortensis. Invert Surviv J. 2013;10:50–7.

    Google Scholar 

  58. 58.

    Ju RT, Wei HP, Wang F, Zhou XH, Li B. Anaerobic respiration and antioxidant responses of Corythucha ciliata (say) adults to heat-induced oxidative stress under laboratory and field conditions. Cell Stress Chaper. 2014;19:255–62.

    CAS  Article  Google Scholar 

  59. 59.

    Yao P, Lu W, Meng F, Wang X, Xu B, Guo X. Molecular cloning, expression and oxidative stress response of a mitochondrial thioredoxin peroxidase gene (AccTpx-3) from Apis cerana cerana. J Insect Physiol. 2013;59:273–82.

    CAS  PubMed  Article  Google Scholar 

  60. 60.

    Hu T, Sun X, Zhang X, Nevo E, Fu J. An RNA sequencing transcriptome analysis of the high-temperature stressed tall fescue reveals novel insights into plant thermotolerance. BMC Genomics. 2014;15:1147.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  61. 61.

    Beaulieu M, Gillen E, Hahn S, Pape JM, Fischer K. Behavioural antioxidant strategies to cope with high temperatures: a study in a tropical butterfly. Anim Behav. 2015;109:89–99.

    Article  Google Scholar 

  62. 62.

    Zhang S, Fu W, Li N, Zhang F, Liu TX. Antioxidant responses of Propylaea japonica (Coleoptera: Coccinellidae) exposed to high temperature stress. J Insect Physiol. 2015;73:47–52.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Arking R, Burde V, Graves K, Hari R, Feldman E, Zeevi A, et al. Forward and reverse selection for longevity in Drosophila is characterized by alteration of antioxidant gene expression and oxidative damage patterns. Exp Gerontol. 2000;35:167–85.

    CAS  PubMed  Article  Google Scholar 

  64. 64.

    Landis GN, Tower J. Superoxide dismutase evolution and life span regulation. Mech Ageing Dev. 2005;126:365–79.

    CAS  PubMed  Article  Google Scholar 

  65. 65.

    Pijpe J, Pul N, van Duijn S, Brakefield PM, Zwaan BJ. Changed gene expression for candidate ageing genes in long-lived Bicyclus anynana butterflies. Exp Gerontol. 2011;46:426–34.

    CAS  PubMed  Article  Google Scholar 

  66. 66.

    Beaulieu M, Geiger RE, Reim E, Zielke L, Fischer K. Reproduction alters oxidative status when it is traded-off against longevity. Evolution. 2015;69:1786–96.

    PubMed  Article  Google Scholar 

  67. 67.

    Franke K, Fischer K. Inbreeding interferes with heat-shock response. Heredity. 2015;114:80–4.

    CAS  PubMed  Article  Google Scholar 

  68. 68.

    Karl I, Sørensen JG, Loeschcke V, Fischer K. HSP70 expression in the copper butterfly Lycaena tityrus across altitudes and temperatures. J Evol Biol. 2009;22:172–8.

    CAS  PubMed  Article  Google Scholar 

  69. 69.

    Sørensen JG, Schou MF, Kristensen TN, Loeschcke V. Thermal fluctuations affect the transcriptome through mechanisms independent of average temperature. Sci Rep. 2016.

  70. 70.

    Wakabayashi-Ito N, Belvin MP, Bluestein DA, Anderson KV. Fusilli, an essential gene with a maternal role in Drosophila embryonic dorsal-ventral patterning. Dev Biol. 2001;229:44–54.

    CAS  PubMed  Article  Google Scholar 

  71. 71.

    Zhou J, Chng WJ. Roles of thioredoxin binding protein (TXNIP) in oxidative stress, apoptosis and cancer. Mitochondrion. 2013;13:163–9.

    CAS  PubMed  Article  Google Scholar 

  72. 72.

    Wang Y, Waldron RT, Dhaka A, Patel A, Riley MM, Rozengurt E, Colicelli J. The RAS effector RIN1 directly competes with RAF and is regulated by 14-3-3 proteins. Mol Cell Biol. 2002;22:916–26.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  73. 73.

    Karl I, Lorenz MW, Fischer K. Energetics of reproduction: consequences of divergent selection on egg size, food limitation, and female age for egg composition and reproduction effort in a butterfly. Biol J Linn Soc. 2007;91:403–18.

    Article  Google Scholar 

  74. 74.

    Kawooya JK, Law JH. Role of lipophorin in lipid transport to the insect egg. J Biol Chem. 1988;263:8748–53.

    CAS  PubMed  Google Scholar 

  75. 75.

    Ziegler R, Van Antwerpen R. Lipid uptake by insect oocytes. Insect Biochem Mol. 2006;36:264–72.

    CAS  Article  Google Scholar 

  76. 76.

    Monaghan JR, Epp LG, Putta S, Page RB, Walker JA, Beachy CK, et al. Microarray and cDNA sequence analysis of transcription during nerve-dependent limb regeneration. BMC Biology. 2009;7:1.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  77. 77.

    Suzuki MG, Shimada T, Kobayashi M. Bm kettin, homologue of the Drosophila kettin gene, is located on the Z chromosome in Bombyx mori and is not dosage compensated. Heredity. 1999;82:170–9.

    CAS  PubMed  Article  Google Scholar 

  78. 78.

    Harrison PW, Mank JE, Wedell N. Incomplete sex chromosome dosage compensation in the indian meal moth, Plodia interpunctella, based on de novo transcriptome assembly. Genome Biol Evol. 2012;4:1118–26.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  79. 79.

    Mank JE. Sex chromosome dosage compensation: definitely not for everyone. Trends Genet. 2013;29:677–83.

    CAS  PubMed  Article  Google Scholar 

  80. 80.

    Walters JR, Hardcastle TJ, Jiggins CD. Sex chromosome dosage compensation in Heliconius butterflies: global yet still incomplete? Genome Biol Evol. 2015;7:2545–59.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  81. 81.

    Constantini D, Monaghan P, Metcalfe NB. Loss of integration is associated with reduced resistance to oxidative stress. J Exp Biol. 2013;216:2213–20.

    Article  CAS  Google Scholar 

  82. 82.

    Cziesielski MJ, Liew YJ, Cui G, Schmidt-Roach S, Campana S, Marondedze C, Aranda M. Multi-omics analysis of thermal stress response in a zooxanthellate cnidarian reveals the importance of associating with thermotolerant symbionts. Proc R Soc B. 2018;285:20172654.

    PubMed  Article  CAS  Google Scholar 

  83. 83.

    Wang W, Teng F, Lin Y, Ji D, Xu Y, Chen C, Xie C. Transcriptomic study to understand thermal adaptation in a high temperature-tolerant strain of Pyropia haitanensis. PLoS One. 2018;13:e0195842.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  84. 84.

    Larsen TB. The butterflies of Kenya and their natural history. Oxford: University Press; 1991.

    Google Scholar 

  85. 85.

    Brakefield PM, Reitsma N. Phenotypic plasticity, seasonal climate and the population biology of Bicyclus butterflies (Satyridae) in Malawi. Ecol Entomol. 1991;16:291–303.

    Article  Google Scholar 

  86. 86.

    Van’t Hof AE, Zwaan BJ, Saccheri IJ, Daly D, Bot ANM, Brakefield PM. Characterization of 28 microsatellite loci for the butterfly Bicyclus anynana. Mol Ecol Notes. 2005;5:169–72.

    CAS  Article  Google Scholar 

  87. 87.

    De Block M, Slos S, Johansson F, Stoks R. Integrating life history and physiology to understand latitudinal size variation in a damselfly. Ecography. 2008;31:115–23.

    Article  Google Scholar 

  88. 88.

    Drayton JM, Jennions MD. Inbreeding and measures of immune function in the cricket Teleogryllus commodus. Behav Ecol. 2011;22:486–92.

    Article  Google Scholar 

  89. 89.

    Bradford MM. A rapid and sensitive method for the quantification of microgram quantities of protein utilizing the principle of protein-dye binding. Anal Biochem. 1976;71:248–54.

    Article  Google Scholar 

  90. 90.

    Mbandi SK, Hesse U, Rees DJG, Christoffels A. A glance at quality score: implication for de novo transcriptome reconstruction of Illumina reads. Front Genet. 2014;5:17.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  91. 91.

    Haas BJ, Papanicolaou A, Yassour MJ, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.

    CAS  PubMed  Article  Google Scholar 

  92. 92.

    Langmead B, Salzberg S. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  93. 93.

    Mortazavi A, Brian AW, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  94. 94.

    Pachter L. Models for transcript quantification from RNA-seq. arXiv:1104.3889 [q-bio.GN], 2013.

  95. 95.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    CAS  PubMed  Article  Google Scholar 

  96. 96.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc B. 1995;57:289–300.

    Google Scholar 

  97. 97.

    Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  98. 98.

    Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, Kriventseva EV, Zdobnov EM. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2017;35:543–8.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  99. 99.

    Schneider T, Schmid E, de Castro JV Jr, Cardinale M, Eberl L, Grube M, et al. Structure and function of the symbiosis partners of the lung lichen (Lobaria pulmonaria L. Hoffm.) analyzed by metaproteomics. Proteomics. 2011;11:2752–6.

    CAS  PubMed  Article  Google Scholar 

  100. 100.

    Conesa A, Götz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

    CAS  PubMed  Article  Google Scholar 

  101. 101.

    Alexa A, Rahnenfuhrer J. topGO: Enrichment analysis for Gene Ontology. R package version 2.30.0; 2016.

Download references


We thank Dr. Stephanie Bauerfeind for help with animal rearing, Vivian Herkules for assistance with RNA extractions, and Prof. Karlhans Endlich and Dr. Martin Haase for granting access to their laboratory facilities. We thank two anonymous reviewers for providing constructive criticism that helped to improve the quality of the manuscript.


We acknowledge financial support from the German Research Council (DFG Fi 846/6–1) to KF. The funding body played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

All data generated or analysed during this study are included in this published article [and its supplementary information files], except for sequence data. All sequence data will be deposited in the NCBI Sequence Read Archive (SRA) database.

Author information




CW, KFi, and VO designed and conceived the experiment. IK performed all rearing experiments, collected life history data, performed physiological measurements, and isolated and purified RNA. IK and KFi analyzed life history and physiological data. MS and TBC assembled the transcriptome, performed the transcriptome quality checks, and computed transcript abundances. MS statistically identified differential expression. KFr and MS performed alignments and assignments to transcript clusters. KFr performed annotations and interpreted expression data. BF performed the enrichment analyses and assisted in the interpretation of expression data. CL and KR constructed the tree maps. IK, KFi, and KFr interpreted the data and wrote the manuscript. All other authors contributed to manuscript writing and editing. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Klaus Fischer.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file1:

Table S1. Distribution of transcript clusters across taxa. (DOCX 13 kb)

Additional file 2:

Data S1. This data file contains all original data from the differential expression analysis based on cluster of orthologous groups and eukaryotic orthologous groups. (XLSX 1045 kb)

Additional file 3:

Figures S1-S3. Figure S1. Gives an overview of the specific genes most strongly affected by the factors sex and feeding regime. Figure S2 depicts the functional annotation of transcripts being down- or up-regulated in females relative to males, and being down- or up-regulated under food ad libitum relative to food restriction. Figure S3. shows a multidimensional scaling plot for detecting outliers. Two outlier samples were assumed to have a mislabelled sex and were therefore excluded from further analyses. (DOCX 875 kb)

Additional file 4:

Data S2. This data file contains all original data from the enrichment analyses. (XLSX 17 kb)

Additional file 5:

RNA isolation and purification are described here in detail. (DOCX 13 kb)

Additional file 6:

Table S2. The table shows the distribution of protein coverage for the Trinity and SOAP assemblies of the Bicyclus anynana transcriptome. Trinity produced a better assembly such that only Trinity was used in further analyses, whereas the SOAP assembly was discarded. (DOCX 13 kb)

Additional file 7:

Here a justification is given for analysing males and females in a joint analysis. (DOCX 14 kb)

Additional file 8:

Data S3. This data file contains all original data from the differential expression analysis based on gene ontology terms. (XLSX 2402 kb)

Additional file 9:

This file gives an overview of the results based on gene ontology terms. The assembled transcript clusters were annotated against both, cluster of orthologous groups/eukaryotic orthologous groups (COG/KOG) as well as gene ontology terms. As both annotations gave very similar results, we only present the results based on the COG/KOG system in the main text. (DOCX 461 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Franke, K., Karl, I., Centeno, T.P. et al. Effects of adult temperature on gene expression in a butterfly: identifying pathways associated with thermal acclimation. BMC Evol Biol 19, 32 (2019).

Download citation


  • Bicyclus anynana
  • Heat tolerance
  • Oxidative stress
  • Phenotypic plasticity
  • RNAseq
  • Transcriptome