The impact of metabolic plasticity on winter energy use models

Understanding the energetic impact of climate change is an important aspect of global change biology, which can reveal critical patterns of vulnerability that temperature alone cannot detect (Dillon et al., 2010; Fitzpatrick et al., 2019). Ecophysiological energy use models predict energy use by mapping body temperature (often estimated from microclimate temperature) onto energetic expenditure using metabolic rate–temperature relationships (Kearney and Porter, 2020; Sinclair et al., 2016, 2013). However, these models assume that metabolic rate–temperature relationships are static, an assumption that is clearly violated given the large amount of phenotypic plasticity that occurs in response to microclimate variation or developmental change (Sinclair, 2015; Sinclair et al., 2016). The impact of plastic variation in metabolic rate–temperature curves on energy use modelling has not been well explored, which may give rise to misleading energy estimates, or overlook important periods of energetic drain. Here, we examined the importance of incorporating metabolic plasticity into winter energy use models in dormant overwintering insects.

Energy conservation is particularly relevant to insect fitness in winter, when finite energy stores must sustain life and also fuel subsequent reproduction (Hahn and Denlinger, 2011; Sinclair, 2015). To conserve limited energy stores, many insects overwinter in diapause – a programmed dormancy characterized by arrested (or slowed) development and metabolic suppression (Koštál, 2006; Wilsterman et al., 2021). Insect metabolic rate during diapause is a fraction (varying from 10% to 85%) of active metabolic rate of comparable developmental stages (Hahn and Denlinger, 2011; Ragland et al., 2009). Metabolic rate changes dynamically as diapause progresses (Koštál, 2006), reaching a minimum several weeks after onset (Toxopeus et al., 2021), followed by a gradual increase towards the time of diapause termination, usually in spring (Gray et al., 1995; Lester and Irwin, 2012).

Metabolic rate–temperature relationships during winter can be altered in two main ways: (1) by changing the intercept, which represents overall metabolic intensity, or (2) by changing the slope, representing thermal sensitivity (Terblanche et al., 2009). Both metabolic intensity and thermal sensitivity change as a result of phenotypic plasticity, including developmental plasticity (e.g. by proceeding through stages of diapause) and acclimation or acclimatization to winter microclimates (Gray et al., 1995; Lester and Irwin, 2012; Williams et al., 2012b). Variable winter microclimates lead to decreased thermal sensitivity, while warmer winter microclimates lead to decreased metabolic intensity, both leading to energetic savings in energetically demanding environments (Sgolastra et al., 2010; Williams et al., 2015, 2012b). Developing a general method to incorporate these widespread patterns of phenotypic plasticity into energy use models may improve our ability to accurately predict winter energy use.

Chrysomela aeneicollis (Schaeffer 1928) populations in the Sierra Nevada mountains are an important model system for understanding the physiological and genetic basis of responses to climate change (Dahlhoff et al., 2019, 2008; Rank and Dahlhoff, 2002). Freeze-tolerant adults overwinter in diapause in the soil, often beneath snow, for up to 8 months of their 1 year life cycle (Boychuk et al., 2015). Snow cover modifies the thermal environment of soil by buffering from cold air temperatures (Pauli et al., 2013), providing a relatively warm and stable thermal environment for overwintering insects that can impact energetics and fitness (Irwin and Lee, 2003). Ecophysiological energy use models suggest that energetic costs of winter decrease across elevation, and that snowy years are more energetically demanding than dry years as a result of longer winter periods (Roberts et al., 2021), but previous energy use models used a single metabolic rate–temperature curve with no plasticity taken into account. In this study, we addressed two main objectives: first, we empirically tested how metabolic rate–temperature relationships of C. aeneicollis change throughout dormancy and in response to thermal acclimation to ecologically relevant microclimates associated with winter snow. Second, we used ecophysiological models to assess how these modifications will impact energy use estimates, by predicting energy use using different metabolic rate–temperature curves. We then compared these predictions with actual energy use derived from empirical measurements of energy reserves.

Beetles were collected in August 2018 from Mosquito Flat (3067 m) in the Rock Creek Drainage (37°26′25.8″N, 118°44′46.0″W) in the Eastern Sierra Nevada Mountains in CA, USA. Beetles were then housed in incubators (MIR-154-PA incubators, Panasonic Scientific, Wood Dale, IL, USA) in the laboratory at a 20/4°C 12 h:12 h day:night cycle until they entered dormancy, following protocols in Roberts et al. (2021). Once beetles entered dormancy (determined by cessation of feeding in the presence of food; between 21 September and 11 October 2018), they were housed in 50 ml conical tubes filled to around 40% with coconut husk substrate, and moved to an incubator held at 1°C in constant darkness. On 1 November 2018, the beetles (N=120) were haphazardly divided between a constant and variable acclimation treatment in separate incubators under constant darkness. The constant treatment simulated conditions beneath the snow (constant 1°C), while the variable treatment simulated uncovered ground (−2.5/−1/0/−1°C in a 6/6/6/6 h cycle). Temperature regimes were selected to mimic late-winter conditions, based on temperature measurements from plots subject to an experimental snow manipulation (Roberts et al., 2021). Beetles were kept in these acclimation treatments, with moisture being added monthly to keep the substrate from drying (around 2–5 ml per month), until their respiration rates were measured.

To capture changes in the metabolic rate–temperature relationship through the end of dormancy, oxygen consumption rate (⁠⁠) was measured using a Sable Systems FoxBox respirometer (Sable Systems International, North Las Vegas, NV, USA) stop-flow respirometry system once a month from February to May, following protocols from Roberts et al. (2021). We chose these time points to best reflect metabolic rate through the transition out of the coldest portion of dormancy. Briefly, at each time point a subset of 15 beetles from each acclimation treatment (constant and variable) were individually placed into a 10 ml syringe and flushed with CO₂- and H₂O-free air, generated using a Drierite–Ascarite–-Drierite column, and then incubated at −1, 4 or 9°C for 48 h, or 20°C for 24 h before oxygen consumption was measured. The order of temperature exposures for the three lowest temperatures was randomized, but always finished with the 20°C measurement to prevent any downstream impacts of prolonged exposure to warm temperatures on respiration rate. Beetles were given a minimum of 48 h to recover between each measurement, during which time they were returned to their acclimation treatment following Lake et al. (2013). All measurements were taken at room temperature, with each syringe being removed from a given incubator within 20 s of the measurement being taken. was calculated from oxygen measurements by integration of the area under the curve in Expedata (Sable Systems International), then corrected for air left in the syringe and divided by time spent in the syringe [=(5 ml/3 ml)/time] (Lighton, 2018). After each beetle was measured at all temperatures, beetles were weighed then frozen for biochemical analysis.

All analyses were performed in R 4.0.2 (http://www.R-project.org/) unless otherwise specified.  was log-transformed (log_e) prior to analyses to approximate a linear relationship with temperature. Metabolic rate–temperature curves were fitted to the data using a linear mixed effect model with the lmer function from the lme4 package. We used a Kenward–Roger approximation in the lmerTest package to approximate degrees of freedom and P-values (Bates et al., 2015; Kuznetsova et al., 2017). Model selection was done by starting with a full factorial model and eliminating terms using a backwards reduced Kenward–Roger degrees of freedom method with the lmerTest package, which eliminates fixed effects via marginal contrasts of degrees of freedom (Table S1). The initial model included  as the dependent variable and month, measurement temperature and acclimation treatment as fixed effects, and mass as a covariate, along with their interaction terms. As  measurements were repeated measures, beetle identity was included in the model as a random effect. The final model included month, measurement temperature, acclimation treatment, mass and the interaction of acclimation treatment and temperature as fixed effects, and individual as a random effect (Table S2). Once the best fit model was chosen and significant terms were identified, we then made models that identified the slope and intercepts of monthly values and treatment independently. Slopes and intercepts were extracted and used to approximate thermal sensitivity and metabolic intensity in energy use models.

Hourly oxygen consumption throughout the winter was calculated using Eqns 1–3, using temperatures (t) from the constant and variable acclimation treatments. Monthly energy use was calculated by summing hourly oxygen consumption, and converted to lipid consumed assuming 2 l oxygen per 1 g lipid metabolized (Schmidt-Nielsen, 1997), giving rise to monthly lipid use estimates under each model, for each acclimation treatment.

We quantified lipid (triacylglyceride) stores of beetles (15 per month per treatment) using thin-layer chromatography coupled to a flame ionization detector (TLC-FID, Iatroscan MK-6s TLC-FID Analyzer, Shell-USA, Spotsylvania, VA, USA), with cholesterol as an internal standard (Roberts et al., 2021; Williams et al., 2011). We tested whether lipid stores changed through time by ANOVA in which month, acclimation and their interaction were included as independent variables and lipid stores as the dependent variable.

To validate model predictions of lipid use against these empirical measurements of lipid stores, we subtracted the monthly lipid use estimates for each model (Eqns 1–3) from the average February lipid stores for constant and variable acclimation separately, giving us predicted monthly lipid stores. We expressed cumulative model errors (over the entire period of February to May) as the absolute difference between lipid stores of individual beetles in May and predicted lipid stores in May under each model, separately for constant and variable incubators. Model accuracy was compared between models using ANOVA with absolute errors as the response variable and model as the predictor, separately for each acclimation treatment.

Beetle metabolic rate–temperature relationships changed through time, and in response to thermal acclimation. Beetle metabolic rate increased between February and May, resulting in a gradual increase in metabolic intensity with no corresponding change in thermal sensitivity (F_3,108=32.27, P<0.0001; Fig. 1; Table S2). A similar increase of metabolic intensity during diapause development has been observed in pine beetles (Lester and Irwin, 2012), gypsy moths (Gray et al., 1995) and solitary bees (Sgolastra et al., 2010) but not goldenrod gall flies (Irwin et al., 2001), and occurs during the transition from diapause to post-diapause quiescence (Hahn and Denlinger, 2011; Lester and Irwin, 2012; Ragland et al., 2009). We thus consider it likely that the increase in metabolic rate through time results from developmental plasticity due to diapause progression (Koštál, 2006), although we cannot rule out some influence of other factors that change through time over winter (e.g. long-term low temperature acclimation). Paralleling our results of no change in thermal sensitivity, thermal sensitivity remains consistent throughout diapause in gypsy moth eggs (Gray et al., 1995), but increases towards the end of diapause in fall webworm pupae (Williams et al., 2015), suggesting that changes in thermal sensitivity throughout diapause may be less important and consistent than changes in overall metabolic intensity.

Beetles acclimated to variable temperatures had higher metabolic intensity than beetles acclimated to constant conditions (F_1,237=9.74, P=0.0020; Fig. 2C; Table S2), and lower thermal sensitivity (F_1,332=5.48, P=0.0198; Fig. 2B; Table S2). Variable winter temperatures also resulted in lower thermal sensitivity in overwintering larval Lepidoptera (Williams et al., 2012b). Thermal regimes approximated soil temperatures below snow (constant) or with snow removed (variable) at a high elevation site in the Eastern Sierra Nevada mountains (Roberts et al., 2021), so it is likely that natural variation in snow cover would elicit a plastic response in metabolic rate of these beetles.

Beetle lipid stores decreased linearly through time (F_3,102=2.87, P=0.040; Table S3), and did not differ between constant and variable acclimation treatments (F_1,104=0.008, P=0.93, Fig. 3A,B; Table S3). All models predicted energy use accurately, with estimates falling within the 25th to 75th percentile of empirically measured beetle lipid stores, with the exception of the constant acclimation treatment in May (Fig. 3A). This illustrates that the base model used in Roberts et al. (2021) predicts energy use accurately under laboratory, as well as field, conditions. Neither the dynamic model nor the acclimation model, incorporating population-level plasticity in response to acclimation, improved the accuracy of model estimates of whole winter energy use relative to the base model used in Roberts et al. (2021) (F_2,84=0.33, P=0.72; Fig. 3C,D; Table S4). Thus, incorporating population-level plasticity in metabolic rate (as a result of acclimation or developmental plasticity due to diapause progression) did not improve population-level estimates of energy use from February to May. This suggests that ecologically relevant magnitudes of plasticity will not compromise population-level estimates of energy use.

The mismatch between predictions of the dynamic model and empirical energy reserve estimates was striking. The dynamic model predicted accelerating rates of energy use throughout diapause, as a result of the observed increase in metabolic intensity, but we did not see any increase in rates of empirical lipid depletion throughout diapause (if anything, rates of lipid depletion slowed as diapause progressed, particularly in the constant conditions). This may have occurred as a consequence of a physiological shift that decoupled oxygen consumption and lipid use (e.g. switching to an alternative fuel source). This underlines the hazards in using thermal performance curves measured over short time scales to estimate long-term responses (Kingsolver and Buckley, 2017; Sinclair et al., 2016). It is possible that our ability to detect differences in energy use that correlate with changes in the metabolic rate–temperature curves would be improved by an individual-based approach, which would require measuring metabolic rate and energy reserves repeatedly on individual organisms (e.g. by qNMR; Guglielmo et al., 2011), but these experiments are technically challenging with small insects at present.

The magnitude of plasticity that we observed between months or acclimation treatments was dwarfed by the magnitude of variation in individual energy stores. The largest magnitude of plasticity in our experiment was developmental (presumably due to diapause progression) and, in the constant acclimation treatment, the difference in predicted energy use resulting from developmental plasticity was 0.04 mg (base versus dynamic model; Fig. 3A). However, there was large individual variation in observed lipid stores (s.d.=0.15 mg), swamping the magnitude of plasticity. Based on observed levels of variability, the magnitude of difference between model estimates would have to exceed ∼0.18 mg of lipid in order to detect a statistically significant effect (power analysis; power=0.8, n=15, α=0.05). To detect the observed predicted difference in energy use between the base and dynamic model in the constant environment (0.04 mg lipid), we would need a sample size of 112 per treatment (or 7.5 times what was used in this study). The mean (±s.d.) mass of beetles was 15.7±3.9 mg, and mean lipid content was 0.42±0.21 mg, making the standard deviation 25% of the mean. This is well within the 20–38% variance seen in across four other species of overwintering beetles (Lehmann et al., 2020; Watanabe and Tanaka, 2000) and the 20–41% variance in lepidopterans (Marshall and Sinclair, 2012; Williams et al., 2015, 2012b, 2011), indicating that individual variability is worth considering across insect taxa. The large degree of variation among individuals in lipid stores is likely to have more profound impacts on energetic stress in winter than the magnitude of plasticity that we observed. Individual lipid stores are vital in determining an organism's fasting endurance in winter (Trondrud et al., 2021), and may provide an interesting avenue for future research.

The degree of metabolic plasticity that we observed may be more impactful in natural thermal environments, which exhibit greater variation in winter energetic demands than did our two acclimation treatments (i.e. as a result of a wider range of temperatures or longer overwintering periods). Energy use estimates from models incorporating developmental plasticity (increase in metabolic intensity) would be particularly sensitive to warm temperatures because of the non-linear relationship of metabolic rate with temperature (Jensen's inequality; Ruel and Ayres, 1999) combined with the observed increase in metabolic intensity. This could have large impacts on overwinter energy use, where warm temperatures at the end of winter can heavily influence energy expenditure (Roberts et al., 2021). Adding the developmental plasticity to model estimates of winter energy use from Roberts et al. (2021) led to an average increase of 24.3% or 0.149 mg lipid, and when including thermal acclimation there was an 11.0% or 0.064 mg increase. As an ∼0.18 mg lipid difference is needed to detect an effect given intraspecific variability in lipid stores, acclimation plasticity will have negligible impacts, but developmental plasticity may come close to causing detectable differences in energy use in more ecologically relevant thermal regimes.

The energetic costs of winter are a critical and understudied aspect of understanding biotic responses to climate change (Fitzpatrick et al., 2020; Sinclair, 2015; Williams et al., 2012a). Metabolic rate–temperature relationships in dormancy are plastic, changing as a result of developmental plasticity and thermal acclimation. We developed a novel method to incorporate these types of plasticity into energy use models in a generalizable way, but this did not improve our model estimates relative to empirical lipid quantifications. We instead found that individual variation in lipid stores far overshadowed the effects of plasticity on energy use. This suggests that estimating winter energy use based on a single metabolic rate–temperature relationship may give us reasonable energy use estimates, and that individual variability in lipid stores is a critical component that is often overlooked, but may provide valuable insight that will help us to understand the energetic impacts of environmental change.

We would like to thank two anonymous reviewers and the editors of Journal of Experimental Biology for their constructive feedback that helped improve the manuscript. We would also like to thank Andre Szejner-Sigal, Lisa Treidel and the rest of the Williams’ Lab for helpful discussion and feedback on ideas; Carol Blanchette and the staff at Sierra Nevada Aquatic Research Laboratory for help with fieldwork logistics; and Nathan Rank, Elizabeth Dahlhoff and Jonathon Stillman for support on experimental design and analysis.