Developmental plasticity in metabolism but not in energy reserve accumulation in a seasonally polyphenic butterfly

Temperate and high latitudes are characterized by the seasonal cycle between favourable and unfavourable conditions. Organisms living in these environments year-round must be able to tolerate adverse conditions. In short-lived organisms, survival over the adverse season is usually possible only in a specific life stage that is able to enter a state of diapause, dormancy or hibernation that can tolerate adverse conditions and survive long periods on stored energy reserves (Tauber et al., 1986; Danks, 1987; Hahn and Denlinger, 2007; Varpe, 2017). In organisms that can complete several generations during the favourable season, not all generations experience the adverse winter conditions. Hence, a developmental switch (Nijhout, 2003) resulting in alternative developmental pathways – one leading to overwintering and delayed maturation and one leading to direct maturation – evolves. Because individuals completing their life cycles during the favourable season face different selection pressures compared with overwintering individuals, the evolution of alternative developmental pathways should result in phenotypic divergence (polyphenism) between directly developing and overwintering individuals (Moran, 1992; Kivelä et al., 2013, 2016a).

Seasonal polyphenism – the expression of different phenotypes in different generations through the season – in morphology and life history is widespread in arthropods, particularly in insects (e.g. West-Eberhard, 2003; Arbačiauskas, 2004; Simpson et al., 2011). Insect seasonal life history polyphenisms are associated with alternative developmental pathways of diapause (overwintering) and direct development (e.g. Wiklund et al., 1991; Blanckenhorn and Fairbairn, 1995; Aalberg Haugen et al., 2012; Friberg et al., 2012; Välimäki et al., 2013; reviewed by Kivelä et al., 2013). Alternative life histories in the diapause and direct development pathways have a physiological basis, like life history phenotypes in general (e.g. see contributions in Flatt and Heyland, 2011). These physiological differences culminate in the life stage that is capable of entering diapause, because surviving a long adverse period in the dormant state of diapause necessitates a very different physiology compared with direct development into a reproductive adult (e.g. Han and Bauce, 1998; Hahn and Denlinger, 2007, 2011; Liu et al., 2007, 2016; Rozsypal et al., 2013; Lehmann et al., 2016, 2018).

Life history divergence between the diapause and direct development pathways may, however, start accumulating before the attainment of the developmental stage that is capable of diapause. For instance, individuals entering direct development often have shorter development times and higher growth rates compared with individuals that will enter diapause later in development (Wiklund et al., 1991; Blanckenhorn and Fairbairn, 1995; Aalberg Haugen et al., 2012; Friberg et al., 2012; Välimäki et al., 2013; reviewed by Kivelä et al., 2013). Such pre-diapause trait differentiation is expected because the developmental switch between diapause and direct development is typically already triggered by environmental cues of seasonal change before the attainment of the diapause-capable life stage (Tauber et al., 1986; Friberg et al., 2011). Hence, we would expect divergence in life history and physiology to start accumulating immediately after the induction of the developmental pathway, but the developmental timing of trait divergence between diapause and direct development remains insufficiently understood. Depending on species, life history traits start diverging between the alternative developmental pathways either early (Nylin et al., 1989; see also Lindestad et al., 2019) or late in juvenile development (Friberg et al., 2012; Esperk et al., 2013). Development and growth rates are expected to be associated with other biological rates, such as metabolic rate (e.g. Brown et al., 2004; but see Glazier, 2015 for a potential lack of association between metabolic rate and other biological rates). Nevertheless, it remains unknown whether the relatively short development time and high growth rate under direct development are associated with a high metabolic rate that also diverges after the induction of alternative developmental pathways. Alternatively, directly developing individuals may just follow the same metabolic and growth trajectories as diapausing individuals, but metamorphose earlier.

Besides metabolic rate, the diapause and direct development pathways are expected to also show other pre-diapause physiological differences. Overwinter survival necessitates energy reserves that can fuel metabolism through the entire dormant period during which feeding is either not possible at all or insufficient to fully fuel metabolism (Tauber et al., 1986; Danks, 1987; Hahn and Denlinger, 2007, 2011; Lehmann et al., 2016, 2018; see Liu et al., 2016 for energy reserves for summer diapause). These energy reserves must be accumulated before the beginning of diapause, which may result in differences between pathways in the allocation of resources to building reserves (Tauber et al., 1986; Danks, 1987; Hahn and Denlinger, 2007, 2011; Liu et al., 2007).

Here, we address potential pre-diapause physiological divergence between the diapause and direct development pathways using the green-veined white butterfly, Pieris napi (Linnaeus 1758), as study system. Pieris napi produces at least two generations per year in a large part of its Palearctic distribution (Wiklund et al., 1991; Tolman and Lewington, 2002; Kivelä et al., 2015). It overwinters exclusively as a pupa, but the induction of the developmental pathway takes place during the penultimate larval instar in response to photoperiod and temperature (Friberg et al., 2011). Individuals following the diapause and direct development pathways have different life histories, including a shorter larval development time under direct development than under diapause (Wiklund et al., 1991; Friberg et al., 2011; Kivelä et al., 2015). However, we do not know whether the relatively high rates of development and growth under direct development are associated with a high metabolic rate, and whether metabolic rate diverges between the diapause and direct development pathways at the same time as the divergence in development and growth rate appears. This information is vital to fully understand the physiology of alternative developmental pathways. In P. napi, diapause pupae contain more storage lipids than pupae that develop directly into adult when measured 6 days after pupation and, in this same phase of pupal development, lipidomes and metabolomes also differ between the developmental pathways (Lehmann et al., 2016, 2018). It remains unknown whether these differences appear already in the larval stage, during which lipid reserves are accumulated in P. napi (Kivelä et al., 2016b) and many other insects (Tauber et al., 1986; Hahn and Denlinger, 2007, 2011). Interestingly, however, metabolic rate and the metabolome are similar in diapausing and directly developing P. napi pupae directly after pupation (Lehmann et al., 2016, 2018). This is surprising given that the induction of the developmental pathway takes place already in the larval stage (Friberg et al., 2011) and some physiological differences appear in the final larval instar as the physiology of pupation induction differs between diapause and direct development in this species (Kivelä et al., 2017). Hence, P. napi clearly has potential for inter-pathway physiological differentiation already in the larval stage, but whether there are differences in metabolic rate and reserve accumulation in the larval stage needs to be elucidated.

We measured metabolic rates of growing P. napi larvae around the induction of either pupal diapause or direct development in the penultimate larval instar and continued the measurements through the final instar. This design enabled us to assess whether metabolic rate is higher under direct development than diapause and whether potential metabolic divergence between diapause and direct development begins directly after the induction of alternative developmental pathways has taken place in the penultimate instar or later. We also measured body lipid content through the final larval instar to compare energy reserve accumulation in the alternative developmental pathways. The fat body is the major energy reserve in insects (Hahn and Denlinger, 2007, 2011; Arrese and Soulages, 2010), and because lipid is the main fuel during diapause in P. napi (Lehmann et al., 2016), body lipid content should be a good indicator of energy reserve size. Together, these two traits provide insight into the metabolic strategy of diapause in P. napi, and whether the energetic challenge is met through strong metabolic suppression, extensive lipid accumulation or a combination of both.

Six diapause-generation (i.e. overwintered) adult P. napi females were captured in Stockholm, Sweden (59°22′N, 18°4′E), and individually placed in 1 litre containers provided with the natural larval host plant Alliaria petiolata as an oviposition substrate in a laboratory. Eggs were frequently monitored and neonate larvae were individually moved to 0.5 litre containers with ad libitum A. petiolata. Each family (offspring of a single female) was split between two climate cabinets (Termaks Series KB8400L; Termaks, Bergen, Norway): one with diapause-inducing short-day conditions (12 h day length at 20°C) and the other with long-day conditions that induce direct development (23 h day length at 20°C). The larvae were monitored daily, and fresh A. petiolata was provided ad libitum. Larvae were measured as soon as they finished the moult to the penultimate instar. At this stage, siblings from both diapause- and direct development-inducing conditions were further split between control and respirometry groups: the first larva to finish the moult to instar IV within each family was randomized to one group; thereafter, every second sibling to finish the moult was assigned to alternative groups.

Larvae in the control group were weighed (Precisa XB 120A; Precisa, Dietikon, Switzerland) twice per day with a ca. 8 h interval during photophase in the short-day climate cabinet. This procedure continued through instars IV and V until the larvae became prepupae. In the respirometry group, the procedure of measuring the larvae was otherwise similar as in the control group, but weighing of the larvae was always preceded by a measurement in a respirometer (larvae were weighed immediately after the respirometric measurement; see below for details of respirometric measurements). Larvae were kept in the climate cabinets under the assigned diapause- or direct development-inducing conditions all the time, except during the measurements.

Prepupae were monitored twice a day for pupation and, approximately 48 h after pupation, the pupae were weighed and sex was determined. Pupae were housed under the same conditions in the same climate cabinets as during the larval period. Pupae were monitored daily for adult eclosion. At 20°C, directly developing individuals typically eclose within 2 weeks after pupation. Therefore, we considered all individuals that did not eclose and showed no signs of development within 3 weeks after pupation to be in diapause. All such individuals were still in diapause after 6 months.

The experiment was repeated with two temporally separate cohorts so that the start of the experiment in the second cohort was 14 days later than that in the first cohort. Different families were included in the two cohorts. The first and the second cohorts included two and four families, respectively. In the first cohort, 10 control group larvae and 10 respirometry group larvae were exposed to both diapause- and direct development-inducing conditions (5 larvae per family in each combination of developmental pathway and group). In the second cohort, the group-specific number of larvae was 12 (2–5 larvae per family in each combination of developmental pathway and group). Owing to mortality in the first cohort (two direct-development control larvae), we obtained data from 42 control and 44 respirometry larvae through the instars IV and V in total.

We followed the ethical guidelines of using animals in research (ASAB, 2006) and Swedish legislation. No permissions were needed in this study, as we used a non-endangered invertebrate.

Differential mode flow-through respirometry was used to measure carbon dioxide (CO₂) production and oxygen (O₂) consumption of individual larvae. Outdoor air was sucked into the system by an SS-4 pump (Sable Systems, Las Vegas, NV, USA) set to 500 ml min⁻¹, and scrubbed of H₂O using Drierite (WA Hammond Drierite, Xenia, OH, USA). Then, the air inflow was split into two lines. Both lines were pushed at a steady rate of 150 ml min⁻¹ using two separate mass flow controllers (840 Series; Sierra Instruments, Monterey, CA, USA). The first line went into a MUX multiplexer (Sable Systems) and through the measurement chamber (volume 10 ml), located in a climatic cabinet (Panasonic MIR-154-PE, Panasonic, Osaka, Japan) set to 20°C. Thereafter, it was again scrubbed of H₂O with magnesium perchlorate (Sigma-Aldrich) before entering the sample line of a Li-7000 CO₂ analyser (LiCor, Lincoln, NE, USA). The second line proceeded in the same way, mimicking the exact length of the sample line (including an empty measurement chamber), before entering the reference line of the Li-7000 CO₂ analyser. Both lines then proceeded through ascarite CO₂ scrubbers (Acros Organics, NJ, USA), and entered an Oxzilla FC-2 O₂ analyser (Sable Systems), after which air was ejected to the room. Both CO₂ and O₂ analysers were calibrated before the start of the experiment. Differential CO₂ and O₂ were calculated by subtracting the output of the reference line from the output of the sample line. The sampling rate was 1 Hz for all measurements. The measurement protocol entailed a 2 min baseline measurement of an empty chamber (preliminary experiments showed that this was long enough), followed by an 8 min measurement of the chamber containing the animal, and ending with a 3 min measurement of the empty chamber (see Fig. S1). In the program ExpeData (version 1.7.30), the raw data were first lag corrected to synchronize CO₂ and O₂ traces and then drift corrected (there was minimal drift during the 15 min measurements). The trace was then baseline corrected against the empty chamber value, fractioned and multiplied with flow rate to convert the raw output to ml CO₂ h⁻¹ or ml O₂ h⁻¹ (Lighton, 2008).

We aimed to minimize the time larvae had to be without access to food, and the consequent negative effects on growth. Therefore, larvae were not starved before respirometric measurements, yet they were without food during the measurement (lasting approximately 15 min including handling). The number of faecal pellets produced during the measurement was recorded to take into account potential variation in the digestive cycle and CO₂ released from the faeces. We did not measure movement of larvae during measurements as a previous study using a similar setup showed that larvae move little in the tubes and, in that study, the movement effect was either statistically insignificant or very small (Kivelä et al., 2016b).

We calculated the respiratory quotient (RQ) only from measurements where the signal-to-noise ratio was >5, as these measurements were considered to be reliable. The signal-to-noise ratio was calculated as the difference in measurement mean and baseline mean, divided by the baseline standard deviation. All CO₂ measurements exceeded this threshold but, for O₂, the signal was generally too weak until the end of the penultimate instar (the smallest larva with a successful O₂ measurement was 18.4 mg).

In addition to larvae used in control and respirometric measurements, we reared extra larvae under the above-mentioned diapause and direct development conditions. Once these larvae entered the final (V) instar, a time point (day during the instar V or a 2-day-old pupa) for freezing was randomized for each individual. The larvae were weighed twice a day until they were frozen in liquid nitrogen at the pre-determined time points. The numbers of frozen final-instar larvae and pupae were 99 and 14 in the diapause pathway, respectively, the corresponding direct development pathway numbers being 94 and 17.

Samples were kept at −80°C until analysis. Water and total lipid content of larvae was measured after chloroform:methanol (chl:meth) extraction (Folch et al., 1957) according to the following protocol. First, larvae were weighed (fresh mass; Precisa XB 120A balance), and then dried for 72 h at 55°C. Larvae were then reweighed (dry mass), placed in small glass vials (20 ml) containing 4 ml chl:meth [2:1 (vol:vol)], crushed thoroughly with a glass rod and left for 72 h in a fume hood at 20°C. Chl:meth effectively dissolves lipids in the sample into the liquid phase. This includes neutral lipids and phospholipids, but also glycerol, carbohydrates and amino acids (Newman et al., 1972), so this method overestimates total lipids (Williams et al., 2011). In P. napi, there appears to be no significant accumulation of these compounds right after pupation regardless of pathway (Lehmann et al., 2018), and, therefore, the method should be suitable for between-pathway comparisons even though the absolute lipid values need to be treated with caution. A previous study, using mass-spectrometric methods, verified that neutral lipids constitute the absolute major component of total lipids in young pupae of P. napi (Lehmann et al., 2016). After the chl:meth submersion, 250 µl of the suspension was retrieved with a syringe (Hamilton, model 1825 RN SYR, 22 s ga) and applied in pre-weighed (Cahn 28 Automatic Electrobalance) aluminium cups (350 µl VWR aluminium round micro weigh). Between each sample, the syringe was rinsed with chl:meth twice. The aluminium cups were left to evaporate in a fume hood for 4 h and reweighed.

Total lipid mass in the sample was calculated by multiplying the mass of lipids in the cup with the dilution factor (4000/250=16). Relative lipid content was calculated as the ratio of total lipid mass to dry mass. Relative water content was calculated as 1–(dry mass/fresh mass).

We used linear mixed-effects models fitted with the maximum likelihood method using the function ‘lme’ (package ‘nlme’; http://cran.r-project.org/package=nlme) in R version 3.3.2 (https://www.r-project.org/) when analysing the data (Kivelä et al., 2019). The multimodel inference approach (e.g. Burnham and Anderson, 2002) was used throughout.

We compared life history traits associated with growth and development (development time through instars IV and V, growth rate in instars IV and V, time from the beginning of instar IV until the attainment of peak larval mass in instar V, peak larval mass, and pupal mass) between the diapause and direct development pathways, and between the control and respirometry groups. Individual-specific growth rate was calculated as the regression slope of ln(mass) on time (days) in the focal instar by including all observations in the regression in instar IV, and the first half of them in instar V to obtain a growth rate estimate for the free growth period before growth deceleration preceding pupation (see Fig. S2). In analyses of life history traits, the ‘global model’ (i.e. the most extensive model considered) included group (control/respirometry), developmental pathway (diapause/direct), sex (female/male) and all interactions among them, as well as cohort (first/second; not included in any interactions) as fixed effects. Family was set as a random effect to control for non-independence of siblings. The set of all meaningful models (i.e. following the principle of hierarchy) simpler than the global model was derived with the function ‘dredge’ (package ‘MuMIn’; https://cran.r-project.org/package=MuMIn; see Table S1, for the sets of models with Akaike weights >0.05). The models, including the global model, were averaged with the function ‘model.avg’ (in package ‘MuMIn’). When needed, we took heteroscedasticity into account by adding weights determined by a variance function (either ‘varIdent’ or ‘varExp’ in package ‘nlme’) into the models. The sets of models with Akaike weights >0.05 are presented in Table S1.

For CO₂ production rate, we analysed instars IV and V separately and included only observations during the growth period in the analysis. Data were ln transformed to facilitate allometric inferences (see Kilmer and Rodríguez, 2016 for arguments favouring this approach over reduced major axis regression). Fixed effects in the global model included centered ln-transformed mass, developmental pathway, sex and all interactions among them, as well as the square of centered ln-transformed mass included in interactions with developmental pathway and sex. Cohort (not included in any interactions) and the number of frass produced during the respirometric measurement (a covariate) were also included in the fixed effects. To control for repeated within-individual measures, random effects included random individual-specific intercepts and slopes in relation to centered ln-transformed mass. The inclusion of random slopes significantly improved model goodness-of-fit over the random-intercepts-only model both in instar IV and V data when measured by Akaike’s information criterion (ΔAIC>2). The inclusion of random individual-specific slopes is also important to avoid overconfident estimates of fixed effects (Schielzeth and Forstmeier, 2009). To reduce correlations among random intercepts and slopes, ln-transformed mass was centered (subtracting the mean from each observation). Centering a variable changes only the mean but not the variance. Therefore, centering of ln(mass) does not affect the regression slopes of other traits in relation to ln(mass), and allometric inferences remain, consequently, unaffected by centering. Model averaging and modelling of heteroscedasticity was performed as explained above, and the sets of models with Akaike weights >0.05 are presented in Table S2.

Variation in RQ was analysed as explained above for CO₂, except that the models included only random individual-specific intercepts because the inclusion of random individual-specific slopes did not improve model goodness-of-fit (AIC values were lower for the random intercept models). Because there was no need to include random individual-specific slopes in the models for RQ, mass was not centered, to make interpretation of the results easier. In instar IV, we added weights calculated using the ‘varIdent’ variance function (in package ‘nlme’) to the models, because modelling the higher residual variance in females than in males improved the estimation accuracy of random effects, yet AIC did not suggest a difference between the weighted and unweighted models (ΔAIC=0.919). The sets of models with Akaike weights >0.05 are presented in Table S3.

In the analyses of body composition data (body water and lipid content) from instar V larvae, fixed effects of the global model included developmental pathway, ln-transformed fresh mass, the square of ln-transformed fresh mass, and interactions between developmental pathway and each of ln-transformed fresh mass and its square. Family was set as a random effect. Weights calculated using the ‘varExp’ variance function (in package ‘nlme’) significantly improved model goodness-of-fit (ΔAIC=43.0 for absolute water content, ΔAIC=11.5 for absolute lipid content, and ΔAIC=22.1 for relative lipid content). We also compared relative water and lipid contents between similarly sized instar V larvae and pupae. Here, only larvae with mass between the minimum and maximum of pupal masses in these data were included in the analysis. The global model included developmental pathway, developmental stage (instar V larva/pupa) and their interaction as fixed effects, and family as a random effect. Model averaging and modelling of heteroscedasticity was performed as explained above, and the sets of models with Akaike weights >0.05 are presented in Tables S4 and S5.

In both the control and respirometry groups, all individuals exposed to the short day length manipulation entered diapause, whereas the long day manipulation invariably induced direct development. Individuals developing directly into adults had higher larval growth rates than individuals entering pupal diapause in both instars IV and V (Table 1, Fig. 1). Development time through instars IV and V was shorter under direct development than under diapause, and directly developing individuals reached peak larval mass in instar V sooner after entering instar IV than those entering diapause (Table 1, Fig. 2). Instar duration was shorter under direct development than under diapause in both instar IV (direct development mean=2.11 days, 95% CI=2.03, 2.19 days; diapause mean=2.62 days, 95% CI=2.52, 2.72 days) and instar V (direct development mean=3.38 days, 95% CI=3.19, 3.57 days; diapause mean=4.64 days, 95% CI=4.49, 4.79 days). The relative difference in instar duration was more pronounced in instar V compared with instar IV (instar duration under direct development is 73% and 81% of that under diapause development in instars V and IV, respectively). Despite having higher growth rate, the fast development of directly maturing individuals led to a low peak larval mass compared with that attained by individuals entering diapause (Table 1, Fig. 2). Males had higher pupal masses than females, and there was also a weak tendency for the pupal mass to be higher under diapause than direct development (Table 1, Fig. 2).

We used whole-organism CO₂ production rate as a proxy for metabolic rate because there is a strong and positive correlation between these two traits (estimated correlation coefficient based on resampling to avoid pseudoreplication=0.961, 95% CI=0.934, 0.980; Fig. S3). Furthermore, O₂ consumption rate, which is required in the calculation of metabolic rate, could not be reliably measured early in instar IV with the set-up used here. Therefore, to facilitate inter-pathway comparisons through the two investigated instars, we base our inferences on CO₂ production rate.

There was no difference between diapause and direct development in CO₂ production in instar IV, but in instar V, CO₂ production was, on average, 1.14 ml CO₂ h⁻¹ higher under direct development than under diapause (Table 2, Fig. 3). CO₂ production rate increased non-linearly through instar IV; a hyperallometric increase early in the instar changed to a hypoallometric increase later in the instar, and the ontogenetic increase in CO₂ production ceased altogether at the end of the instar (Table 2, Fig. 3A,C; see also Fig. S4). In instar V, CO₂ production had a qualitatively similar relationship to body mass as in instar IV, yet the non-linearity of the relationship was weaker (Table 2, Fig. 3B,D).

The RQ in instar IV was independent of developmental pathway (diapause mean=0.873, 95% CI=0.805, 0.940; direct mean=0.884, 95% CI=0.807, 0.961) and body mass, the latter being apparently due to the concentration of data in the later portion of the instar (Table 2, Fig. 4A). The larvae became large enough for reliable measurement of O₂ consumption (and RQ) only in the latter part of the instar IV, which is why we do not have RQ data from the beginning of the instar IV. In instar V, RQ showed a peaked relationship to body mass, RQ values being higher in the second than in the first cohort (Table 2, Fig. 4B). There was no difference in RQ values between diapause and direct development (diapause mean=0.957, 95% CI=0.933, 0.981; direct mean=0.966, 95% CI=0.941, 0.991).

The absolute amount of water in the body did not differ between diapause and direct development and increased nearly isometrically, yet non-linearly, over instar V, the scaling becoming increasingly hypoallometric at the end of the instar (Table 3, Fig. 5A). Absolute mass of lipids in the body did not differ between diapause and direct development either and increased linearly and hyperallometrically through instar V (allometric exponent=1.37; Table 3, Fig. 5B). Relative lipid content (proportion of lipids in dry body matter) peaked in the middle of instar V and did not differ between diapause and direct development (Table 3, Fig. 6).

Relative water content was lower in 2-day-old diapause pupae than in similarly sized instar V larvae or 2-day-old directly developing pupae (Fig. 5C), as indicated by the significant interaction between developmental pathway (diapause/direct) and stage (instar V larva/pupa; Table 3). A similar analysis in the proportion of lipids in dry body matter revealed that the proportion of lipids was higher in instar V larvae than in pupae (Table 3, Fig. 5D), but there was no statistical support for a difference between the developmental pathways.

We found the expected life history divergence between the diapause and direct development pathways through the last two larval instars in P. napi; larval growth and development were faster and body size smaller in individuals destined for direct development rather than pupal diapause. This life history differentiation between the alternative developmental pathways was associated with physiological differentiation, metabolic rate (defined as CO₂ production) being higher under direct development than under diapause. As expected, the metabolic rate of larvae entering direct development and pupal diapause differed only in the final larval instar, and not in the penultimate instar. This is because the penultimate instar is the typical time window during which the developmental pathway is induced in P. napi (Friberg et al., 2011), resulting in a limited post-induction period for phenotypic divergence to become apparent in the penultimate instar. Putting these results together with the observation that the inter-pathway difference in development time became more pronounced in the final larval instar and the fact that the induction of a developmental pathway is largely independent of rates of growth and development in P. napi (Friberg et al., 2012; Kivelä et al., 2015), a high metabolic rate seems to be associated with fast growth and rapid development in individuals set for direct maturation after the induction of the developmental pathway. A high metabolic rate may be the cause for a high growth rate, or vice versa (Glazier, 2015). Either way, it seems that a relatively high metabolic rate under direct development is an integral part of the direct development phenotype. This difference in metabolic rate between the alternative developmental pathways is relatively small in magnitude in terms of CO₂ production rate, but the mass-specific metabolic rate (μW mg⁻¹) is approximately 10% lower in final-instar larvae entering diapause than in direct development (based on a model-averaged estimate; data not shown). Therefore, the metabolic divergence is not trivial, yet there was no association between metabolic rate and accumulation of energy storage. The lack of exaggerated lipid accumulation in larvae entering pupal diapause suggest that energy saving during pupal diapause is primarily attained by metabolic suppression in the pupal stage in P. napi. Although we cannot rule out that the lower metabolic rate found in final-instar larvae preparing for diapause could have some energetic benefit, it seems more likely that actual diapause-related metabolic suppression occurs in the diapausing pupae (Lehmann et al., 2016).

Ontogenetic allometric scaling of metabolism to body size was complex. The scaling of metabolic rate to body size was non-linear, early-instar hyperallometry turning into hypoallometry at the end of an instar (see also Glazier, 2005). This non-linear pattern was consistent both in the penultimate and final instars, yet the non-linearity was exaggerated in the penultimate instar. In an earlier study concentrating only on the diapause pathway, we found a similar non-linear scaling of metabolic rate to body mass in the penultimate instar in the diapause pathway and a linear scaling in the final instar (Kivelä et al., 2016b). A larger sample size in the present study is a likely explanation of why we could detect a clearer levelling off of metabolic rate in the penultimate instar (Fig. 3A,C, Fig. S4) and a clear but weak non-linearity also in the final instar (Fig. 3B,D). Non-linearity of scaling may arise because of an increasing relative amount of metabolically inactive matter in the body towards the end of instars (see Glazier, 2005). Accordingly, the proportion of lipids in body matter largely follow this pattern. Moreover, the present data from the penultimate instar (see also Fig. S4) are consistent with the oxygen-dependent induction of moulting hypothesis, predicting that the tracheal respiratory system, which primarily grows during moults, reaches its maximum capacity at the end of the instar and so sets an upper limit for respiration rate (Greenlee and Harrison, 2004, 2005; Callier and Nijhout, 2011; Kivelä et al., 2016c).

Even though metabolic rate differed between the developmental pathways, there was no difference in the accumulation of lipids between diapause and direct development, although lipids constitute the main energy reserves used during diapause in many insects (Hahn and Denlinger, 2007, 2011; Arrese and Soulages, 2010), including P. napi (Lehmann et al., 2016). This is not exceptional (Hahn and Denlinger, 2007), yet it is unexpected given the ecology of P. napi, where diapausing pupae must survive 5 to 10 months on the energy reserves acquired during the larval stage. In addition, many diapause-generation individuals pupate early in the summer and therefore also experience a long pre-winter warm period in pupal diapause. This may be particularly stressful as higher temperatures should lead to increased fat reserve depletion in ectotherms (Han and Bauce, 1998; Irwin and Lee, 2003; Bosch et al., 2010; Williams et al., 2014; see also Liu et al., 2016). Although there was no inter-pathway difference in the accumulation of energy reserves, final-instar P. napi larvae accumulate a large quantity of lipids in their bodies during this instar, as indicated by the hyperallometric increase in total lipid mass in the body (Fig. 5B). RQ values exceeding 1 suggest that either metabolism is partially anaerobic (e.g. Nielsen and Christian, 2007) or larvae convert carbohydrates to lipids (Zebe, 1953; Wigglesworth, 1972). Conversion of carbohydrates to lipids seems more likely as the proportion of lipids in dry body matter peaks similarly to RQ halfway through instar V (compare Figs 4B and 6; see also Kivelä et al., 2016b). In the pupal stage, diapause pupae had lower water (and possibly higher lipid; see Fig. 5D) contents than directly developing pupae soon after pupation, and this difference emerged during pupation as larval water (and lipid) content did not differ between the developmental pathways. Reduced water content of diapause pupae may be a strategy to improve cold tolerance (see Danks, 2000; Bennett et al., 2005; Rozsypal et al., 2013; Lehmann et al., 2018).

Seasonal polyphenisms in life history and morphology associated with the diapause and direct development pathways have been demonstrated repeatedly in insects (Wiklund et al., 1991; Blanckenhorn and Fairbairn, 1995; West-Eberhard, 2003; Simpson et al., 2011; Aalberg Haugen et al., 2012; Friberg et al., 2012; Välimäki et al., 2013; Kivelä et al., 2013). Here, we show that the phenotypic divergence between the developmental pathways extends to metabolism, and the metabolic divergence appears after the induction of a developmental pathway has taken place. Pre-diapause metabolic divergence between larvae entering diapause and direct development that we demonstrate here persists after pupation (Lehmann et al., 2016). This is necessary because diapause pupae must survive up to 10 months with reserves similar to those of directly developing pupae that instead eclose as adults within a few weeks.

The generally extensive phenotypic divergence between alternative developmental pathways suggests that many components of the physiological and developmental machinery are independent enough for pathway-specific evolution. It is possible that pathway-specific selection changes metabolic rate, which would imply that metabolic rate is a key component of physiology affecting the evolution of life history polyphenisms (see also Arbačiauskas, 2004; Furness et al., 2015; but see Saastamoinen et al., 2013), and possibly life history evolution across species (Brown et al., 2004; Maino et al., 2014). Yet, it also remains possible that metabolic rate responds to changes in growth rate (Glazier, 2015), and that growth rate is the key to phenotypic differentiation between the alternative developmental pathways. Further studies are necessary for a rigorous assessment of these alternatives.

We thank Loke von Schmalensee for help in running the experiment. We are also grateful to two anonymous reviewers and an anonymous editorial board member for commenting on an earlier version of the manuscript.