Genotype‐environment interactions rule the response of a widespread butterfly to temperature variation

Understanding how organisms adapt to complex environments is a central goal of evolutionary biology and ecology. This issue is of special interest in the current era of rapidly changing climatic conditions. Here, we investigate clinal variation and plastic responses in life history, morphology and physiology in the butterfly Pieris napi along a pan‐European gradient by exposing butterflies raised in captivity to different temperatures. We found clinal variation in body size, growth rates and concomitant development time, wing aspect ratio, wing melanization and heat tolerance. Individuals from warmer environments were more heat‐tolerant and had less melanised wings and a shorter development, but still they were larger than individuals from cooler environments. These findings suggest selection for rapid growth in the warmth and for wing melanization in the cold, and thus fine‐tuned genetic adaptation to local climates. Irrespective of the origin of butterflies, the effects of higher developmental temperature were largely as expected, speeding up development; reducing body size, potential metabolic activity and wing melanization; while increasing heat tolerance. At least in part, these patterns likely reflect adaptive phenotypic plasticity. In summary, our study revealed pronounced plastic and genetic responses, which may indicate high adaptive capacities in our study organism. Whether this may help such species, though, to deal with current climate change needs further investigation, as clinal patterns have typically evolved over long periods.

As temperature is an important ecological factor, thermal selection is believed to cause or at least contribute to local adaptation (Reeve, Fowler, & Partridge, 2000). Potential target traits of thermal selection include body size, locomotory performance and coloration. Animals are generally predicted to be larger in cooler environments based on both the Bergmann's rule and the temperature-size rule (Arnett & Gotelli, 1999b;Chown & Gaston, 1999).
Although living in cool environments may constrain life histories (see above), high temperatures may also involve physiological costs.
They have direct effects on cellular processes and metabolism, and may induce oxidative stress and concomitantly costly antioxidant responses (Gillooly, Brown, West, Savage, & Charnov, 2001). Oxidative stress is an imbalance between the production of reactive oxygen species (ROS), typically increasing with metabolism and thus temperature in ectotherms, and antioxidant defences, resulting in higher oxidative damage (Ju, Wei, Wang, Zhou, & Li, 2014;Tumminello & Fuller-Espie, 2013). Important parameters in this context are consequently metabolic rate potentially affecting the production of ROS, the degree of oxidative damage and antioxidant defence mechanisms (Ju et al., 2014;Žagar, Carretero, Marguč, Simčič, & Vrezec, 2018).
In an earlier study using exclusively field-collected individuals of the butterfly Pieris napi, we found clinal variation in body size and wing coloration (Günter et al., 2019). Specifically, we showed that body size and the wings' yellow reflectance decreased whereas wing melanization increased with increasing latitude and altitude.
However, the sources of variation, namely the relative contribution of genetic versus plastic responses, remained unknown in this earlier study and are a scientific field that still has open questions and needs more investigation (for a review, see Crispo, 2008). Against this background, we here set out to investigate clinal variation in P. napi using a common garden design with replicated populations from Italy to Sweden. We investigate traits related to body size, flight performance, wing coloration, physiological status and stress resistance, as local adaptation in such traits is at least partly still not well understood. We hypothesize that genetic variation in such traits is widespread as is plasticity therein and that responses are modulated by interactions between both sources of variation. Specifically, we predict (a) reduced growth rates (countergradient variation), metabolic potential, flight performance and wing melanization but increased heat tolerance and defence mechanisms to fight oxidative stress at lower latitudes (genetic variation), and (b) higher growth rates, metabolic potential, heat tolerance and defence mechanisms to fight oxidative stress but reduced body size, flight performance and wing melanization at a higher compared with a lower developmental temperature (plastic responses).

| Study organism and experimental populations
The green-veined white butterfly Pieris napi L. is a widespread temperate-zone butterfly, occurring throughout Europe and the temperate zone of Asia (Ebert & Rennwald, 1993). In Europe, it is one of the most common butterflies. Nevertheless, it is predicted to suffer from anthropogenic climate change because of its association with moist habitats (Fox et al., 2015). Larvae feed on a variety of Brassicaceae, with Alliaria petiolata Cavara & Grande being the most important one. Brassicaceae are also the preferred nectar plants. Pieris napi is a polyandrous species in which males transfer large nuptial gifts to their female partners (Bergström & Wiklund, 2002). Accordingly, males are larger than females (Wiklund & Kaitala, 1995).

| Sampling locations
For this study, we collected freshly eclosed, spring generation females from a total of nine populations along a latitudinal gradient from northern Italy to Sweden (Figure 1) Sweden: 14ºC) followed a temperature gradient across latitudes, whereas precipitation is higher in Italy than in Germany and Sweden (Table S1). The minimal straight distance between two populations was 73 km, and the total latitudinal gradient spanned ca. 1,660 km.
We collected a total of 74 females from Italy, 94 from Germany and 76 from Sweden between 19 April and 14 June 2016. All females were afterwards transferred to Greifswald University for egg laying.
Of these, 36, 37 and 37 females, respectively, produced sufficient offspring and thus contributed to the experiments below, resulting in a total number of 1,829 offspring.

| Experimental design
Field-caught females were kept individually in translucent 1-litre plastic pots covered with gauze, which were placed into a climate chamber set at a constant temperature of 25°C, 65% relative humidity and a photoperiod of L18:D6. Females were fed ad libitum with water, a 10% sugar solution and additionally flowers (e.g. Sambucus nigra, Taraxacum spec. and Senecio spec.). For oviposition, they were provided daily with a fresh cutting of A. petiolata. The resulting eggs were collected daily, counted and kept separated by oviposition day and female under the above conditions until hatching. Eggs were collected until the death of the females. After hatching, all larvae were transferred individually to translucent plastic pots (250 ml) lined with moist tissue and cuttings of A. petiolata ad libitum. Host plants were replaced on a daily basis. On day 3 after hatching, larvae were randomly divided among two rearing temperatures (18 and 25°C), using a split-brood design. These temperatures were chosen to mimic summer conditions (July) in Germany/Sweden and Italy (Table S1).
Larval growth rate was calculated as the natural logarithm of pupal mass/larval time. One day after adult eclosion, butterflies were subjected to a heat knock-down assay. Therefore, they were individually placed into translucent plastic pots (100 ml) in a randomized block design and exposed to a constant temperature of 43°C (climate cabinet Sanyo MIR-553). The time until physical knock-down, characterized by an inability to move in a coordinated manner, was recorded.
Butterflies were afterwards frozen at −80°C for later analyses.
Frozen butterflies were later thawed, and their adult fresh mass was measured with a digital scale (KERN ABJ-120-4M). Then, the head, wings and legs were removed, and the thorax and abdomen were separated before being weighed. The thorax-abdomen ratio was calculated as the thorax mass divided by the abdomen mass.
Butterfly wings were used to measure wing morphology and melanization. Therefore, we took a photograph of one dorsal forewing and one ventral hindwing per individual under standardized illumination with a PC microscope camera (Veho MS-004 Discovery Deluxe USB Microscope). The ventral hindwing was used as its melanization is known to influence heat gain during lateral basking (Heinrich, 1996), whereas dorsal forewing melanization influences the heat gain during dorsal basking (Kingsolver, 1987). To score wing area, we used the 'lasso' function in Adobe Photoshop CS6.
Wing melanization was defined as the percentage of black wing area ( Figure S2). We used a threshold approach with a fixed value of 128 (on a scale of 0 to 255), turning each pixel on the butterfly wing into either black or white (Adobe Photoshop CS6). We additionally calculated forewing-to-hindwing ratio (forewing area divided by hindwing area), wing loading (total body mass divided by forewing area) and wing aspect ratio to examine wing shape (4 × forewing length 2 divided by forewing area; Berwaerts et al., 2002).

F I G U R E 1
Map of sampling locations of Pieris napi individuals used in the present study. Italy (light grey circles), Germany (dark grey circles) and Sweden (black circles)

| Physiological parameters
We measured the following parameters related to oxidative stress:

| St atis tic al analyses
All data on developmental and morphological traits were analysed using nested ANOVAs with country (Italy, Germany and Sweden), sex, temperature and all resulting interactions as fixed factors, and replicate population as well as family (i.e. the offspring of an individual female) as random effects. Population was nested within country and family within country and population. The physiological parameters PMA, OXY and ROM were analysed with similar ANOVAs, but without the factors sex (as only males were used here) and family (owing to small sample sizes per family). Likewise, MDA, GSH and DNA damage were analysed using similar ANOVAs, including sex but once again excluding family effects because of small sample sizes per mother. Pupal time was ln-transformed, PMA, ROM and OXY were log-transformed, and DNA damage as well as GSH was 1/x-transformed prior to analyses to meet ANOVA requirements. Models were constructed by stepwise backward removal of nonsignificant interactions (Table S4). Owing to the high number of statistical tests, we corrected all p-values for table-wide false discovery rates using the Bonferroni-Holm method (Table 1).
To reduce the complexity of the data set, we additionally performed PCAs. Resulting PC scores were analysed with nested ANOVAs as indicated above (Table S5). Throughout the text, means are given as ± 1 SE. All statistical tests were performed with Statistica 12.0 (Tulsa, StatSoft, OK).

| Developmental traits, morphology and heat knock-down time
After correcting for multiple testing, country significantly affected wing melanization only (Table 1; for details see Table S4).

Country
Popul Family Sex

TA B L E 1
Results of nested ANOVAs for the effects of country (fixed factor), population (random, nested within country), family (random, nested within country and population), sex and temperature (both fixed) on various traits in Pieris napi from a latitudinal gradient Forewing and hindwing melanization were higher in Swedish followed by German and finally Italian butterflies (S > G> I; Tukey HSD after ANOVA; Figure 2a). Additionally, the PCA revealed that Italian individuals were significantly heavier/larger than German and finally Swedish ones (PC1: I > G> S; F 2,6 = 13.5, p = .0060; Figure 2b). PC1 was strongly correlated with all measures of body size (all r-values > .73; Table S5). Despite being larger, Italian larvae had the shortest development time, realized by higher growth rates compared with German or Swedish larvae (PC2: I > G = S; F 2,6 = 43.3, p = .0003; Figure 2c). PC2 was positively related to larval time (r = .713) but negatively to larval growth rate (r = −.749).
Significant temperature and sex differences were found for all traits investigated except for pupal time (no sex difference)

F I G U R E 2
Hindwing melanization (a), abdomen mass (b), larval growth rate (c), wing loading (d) and heat knock-down time (e) of Pieris napi in relation to latitude (Italy, I; Germany, G; and Sweden, S) and rearing temperature (18°C and 25°C). Different letters above boxes indicate significant differences between groups. All data are pooled across sexes. Group samples sizes range between 90 and 333 individuals for all traits and forewing-to-hindwing and wing aspect ratio (no temperature difference) (Table 1). Rearing individuals at the higher compared with the lower temperature resulted in shorter larval and pupal development times; higher larval growth rate, thorax-abdomen ratio and heat stress resistance; lower pupal, adult, thorax and abdomen masses, forewing and hindwing melanization and wing loading; and smaller hindwings and forewings (Figure 2; Figure S6).
Females had significantly longer larval development times; a lower larval growth rate, pupal and thorax mass, thorax-abdomen ratio, hindwing melanization and wing aspect ratio; smaller forewings and hindwings; and higher adult and abdomen masses, forewing melanization, forewing-hindwing ratio, wing loading and heat stress resistance than males.
The above-mentioned general patterns were partly modulated by interactive effects (Table 1). Country-by-temperature interactions were significant for larval time, pupal mass, larval growth rate, abdomen mass, thorax-abdomen ratio and forewing melanization.
In addition, these interactions were often caused by variation in absolute values, while relative changes were similar. For instance, the higher developmental temperature reduced larval time by 4.5 (Italy), 5.1 (Germany) and 5.5 days (Sweden), causing the significant interaction. In relative terms, this is equivalent to a reduction in 30.0, 29.3 and 32.3%. We therefore mention below only those interactions, which reflect differences in relative changes > 5% among countries for a specific trait (see also Figure S6). The higher temperature reduced abdomen mass more strongly in German (−21.8%) than in Italian (−11.6%) or Swedish (−8.8%) individuals ( Figure 2b). Moreover, thorax-abdomen ratio increased by 11.0% in German individuals at the higher temperature, whereas it remained virtually unaffected in Italian (−0.4%) or Swedish (−2.7%) individuals. Regarding heat knock-down time, plasticity was highest in Italian animals, showing an increase by 32.1% at the higher temperature compared with 17.9% in Swedish and 3.0% in German individuals (Figure 2f). Note that the latter interaction was not significant after Bonferroni-Holm correction, but well in the PCA (PC6: F 2,717 = 11.2, p < .0001). Sexby-temperature, country-by-sex and three-way interactions were nonsignificant throughout.

| Latitudinal variation
Wing melanization, development time, wing loading and wing aspect ratio increased, whereas body size, larval growth rate and heat tolerance decreased with increasing latitude. The patterns found for wing melanization, body size and wing loading are in agreement with data based on field-caught individuals from the same cline (Günter et al., 2019). Thus, variation in these traits seems to be buffered across environments. For wing aspect ratio, though, a cline opposite to the one reported earlier was found. This may indicate a large environmental impact on wing aspect ratio.
However, overall there was once again no evidence for strong selection on flight-related morphology across latitudes (Günter et al., 2019). For instance, there was no variation in thorax-abdomen ratio and contradictory evidence for wing aspect ratio versus wing loading. Although the increase in wing aspect ratio with latitude is associated with better manoeuvrability in insects, a higher wing loading is clearly detrimental for flight (Berwaerts et al., 2002;Hassall, 2015).
The finding that individuals from warmer rather than cooler environments were larger challenges both the Bergmann's and the temperature-size rule (Angilletta & Dunham, 2003;Chown & Gaston, 1999). As converse Bergmann size clines have been repeatedly reported (Blanckenhorn & Demont, 2004;Blanckenhorn et al., 2006), we suggest that body size is not directly under thermal selection but perhaps more closely associated with season length. Warmer climates are more beneficial for the development of Pieris species than cooler ones (Angilletta & Dunham, 2003), probably resulting in the pattern observed. Additionally, German individuals appeared to be more plastic than Italian or Swedish ones with regard to abdomen mass. However, this was not found in other body size-related traits, such that no general conclusions seem possible. Despite being larger, Italian larvae exhibited the shortest larval development time, realized by a higher growth rate compared with German or Swedish ones. Again, this finding challenges our a priori prediction (countergradient variation). Thus, populations from warm climates seem to be selected for rapid growth, enabling simultaneously short development time and large size. Perhaps, rapid growth enables squeezing in more generations per year, resulting in compound interest benefits (Fischer & Fiedler, 2002).
The increase in wing melanization towards cooler environments was expected as a darker coloration increases heat absorption, in turn facilitating higher levels of activity (Ellers & Boggs, 2004;Stoehr & Goux, 2008). Therefore, darker individuals will likely more readily reach mating partners, nectar and host plants or escape predators under thermally challenging conditions (Berger, Walters, & Gotthard, 2008). In contrast, strong melanization may cause overheating in warm environments (Kingsolver & Watt, 1983;Van Dyck, Matthysen, & Dhondt, 1997). Thus, clinal variation in melanization results from both, selection to increase it in cooler but to decrease it in warmer environments.
Italian individuals tended to be more heat stress-resistant than German or Swedish ones and also showed higher levels of plasticity in thermal tolerance. Both were expected based on a higher risk of exposure to detrimentally high temperatures in warmer climates in southern Europe. Thus, cold-adapted populations may suffer from future anthropogenic climate change due to both a lower baseline heat tolerance and low levels of plasticity therein (cf. Karl, Janowitz, & Fischer, 2008). The higher heat resistance of Italian butterflies might be related to higher antioxidant defences, as suggested by a trend towards a higher GSH concentration than in German butterflies. This can prevent oxidative damage caused by reactive oxygen species (Hermes-Lima, 2004;Pamplona & Costantini, 2011).
However, the parallel increase in GSH in Swedish animals seems surprising. We speculate that stress imposed by cold temperatures may yield similar effects and therefore adaptations as those imposed by high temperatures (Grim, Simonik, Semones, Kuhn, & Crockett, 2013;Lalouette, Williams, Hervant, Sinclair, & Renault, 2011).

| Effects of developmental temperature
Higher temperatures generally speed up biochemical processes and metabolic rates, resulting in higher growth rates and thus shorter development times as found here (Angilletta, 2009;Blanckenhorn, 1997;Kingsolver & Woods, 1997). Likewise, a reduction in body size respectively mass at the higher temperature was expected according to the temperature-size rule, a near-universal 'law' in ectotherms (Angilletta & Dunham, 2003;Atkinson, 1994).
The temperature-size rule in turn results from thermal effects on behavioural and physiological mechanisms, namely a reduced food intake and a lower efficiency in converting ingested food into body matter at higher temperatures .
The reduced wing loading at the higher temperature likely increases flight performance (Berwaerts et al., 2002;Hassall, 2015).
Why this should be the case is not obvious, as physical limitations on flight are less pronounced at higher temperatures. Thus, the pattern does not seem to imply adaptive plasticity, but is a consequence of an apparently stronger effect of higher temperatures on body mass than body (wing) size. At the higher temperature, thoraxabdomen ratio was overall higher, indicating that especially abdomen mass and therefore storage reserves were reduced. This likely reflects higher metabolic losses at the higher temperature Kingsolver & Woods, 1997). Wing aspect ratio was lower at the higher temperature, which may reduce flying speed and manoeuvrability but increase performance during long-distance flights (Hassall, 2015). The increased wing melanization at lower temperatures likely reflects adaptive plasticity (Karl, Hoffmann, & Fischer, 2010;Peñuelas et al., 2017). Likewise, reduced melanization at higher temperatures may be adaptive for preventing overheating (Kingsolver & Watt, 1983;Van Dyck et al., 1997;Watt, 1968).
Individuals reared at the higher temperature showed an increased heat tolerance, reflecting a well-known pattern of adaptive phenotypic plasticity .
The decrease in potential metabolic activity at the higher temperature is counterintuitive, though similar patterns have been documented before (e.g. Žagar, Holmstrup, Simčič, Debeljak, & Slotsbo, 2018). The mitochondrial respiratory chain is a major source of reactive oxygen species (ROS), which can cause cellular damage (Hermes-Lima, 2004). The reversible inactivation of enzymes, which limits enzymatic processes outside the usual temperature range (Simčič, Pajk, Jaklič, Brancelj, & Vrezec, 2014), may additionally decrease PMA. Such a reduction in enzyme activity is adaptive by reducing the risk of cellular damage due to ROS at higher temperatures and may therefore be an 'upstream' antioxidant response to avoid oxidative stress. These antioxidant responses appear efficient, as none of the considered markers of oxidative damage was significantly affected by thermal conditions.

| Sexual differences
The sexual differences found largely followed expected patterns.
Males showed shorter development times, facilitated by higher growth rates, than females, which is typically explained by selection for earlier male emergence to maximize mating opportunities (Fagerström & Wiklund, 1982;Wiklund & Fagerström, 1977).
Females had a lower pupal and thorax mass as well as smaller forewings and hindwings than males, reflecting the strong selection on large male body size in P. napi, enabling large nuptial gifts (Bergström & Wiklund, 2002). Nevertheless, adult mass was actually higher in females, owing to a much higher mass loss upon metamorphosis in males, which has been described as a cost of the males' rapid development (Fischer, Zeilstra, Hetz, & Fiedler, 2004). Sexual differences in the relative investment into locomotion versus reproduction explain the females' lower thorax-abdomen ratio but higher wing loading (Berwaerts et al., 2002;. With regard to wing coloration, the females' higher degree of forewing melanization results from a second black spot which is missing in males (Tolman & Lewington, 2008). However, hindwing melanization, being relevant to heat absorption (Wasserthal, 1975), was more pronounced in males. This may facilitate flight ability and thus mate location under thermally challenging conditions (Åhman & Karlsson, 2009;Wiklund & Forsberg, 1991).

| CON CLUS IONS
Our common garden experiment revealed clinal patterns in wing melanization, body size and wing loading in P. napi, which are in strong agreement with field-based data (Günter et al., 2019).
Likewise, laboratory and field data concur in that there is no clear evidence for selection on flight-related morphology across latitudes (Günter et al., 2019). Our present data further indicate shorter development times, enabled by higher growth rates, and increased heat tolerance at lower latitudes. These findings suggest that, in P. napi, warmer climates (1) are more beneficial for development, (2) select for rapid growth to enable additional generations per year, (3) select for increased heat tolerance (and potentially antioxidant capacities) and (4) relax selection for thermal melanization. Although (3) and (4) likely reflect thermal adaptation, variation in growth trajectories and body size seems to be more closely related to season length. We thus found clear evidence for genetic adaptation to local conditions. One caveat here is that our study did not control for cross-generational parental effects, but at least all offspring was reared under controlled conditions. Although, in general, parental effects may have important implications (Huestis & Marshall, 2006;Mousseau & Dingle, 1991), their influence is often rather weak as compared with those of the developmental environment (Robinson & Partridge, 2001), in particular in butterflies (Steigenga & Fischer, 2007). Note also that most patterns documented here are in agreement with predictions based on theory and the results of other studies (growth trajectories, Van Doorslaer & Stoks, 2005; heat tolerance, Fischer and wing melanization, Watt, 1968;Van Dyck et al., 1997). In general, plastic responses to temperature were more pronounced than differences among origins, stressing the importance of phenotypic plasticity for dealing with environmental variation. Still, several patterns prevailed under laboratory and field conditions. Some plastic responses were indicative of adaptive phenotypic plasticity, most notably in wing melanization and heat tolerance. Our study highlights that temperature is an important selective agent along latitudinal clines, though the importance of other forces such as season length should not be neglected.

ACK N OWLED G M ENTS
We thank Kristin Franke and Christin Park for help with processing butterflies. This research was funded by the DFG Research Training Group RESPONSE (DFG GRK 2010).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.