Metabolic plasticity can amplify ecosystem responses to global warming

Study system & organisms

The study was conducted in the Hengill valley, Iceland^{13,14,15,16,17,18} (N 64°03; W 21°18), which contains many streams of different temperature due to geothermal heating of the bedrock or soils surrounding the springs (Supplementary Fig. 1). The streams have been heated in this way for centuries³³ and are otherwise similar in their physical and chemical properties^13,18, providing an ideal space-for-time substitution in which to measure species responses after chronic exposure to different temperatures^6,34. Fieldwork was performed in the summers of 2015–2018, between May and July. Stream temperatures were logged every 4 h using Maxim Integrated DS1921G Thermochron iButtons submerged in each stream (Supplementary Fig. 2). The average stream temperature over this study period was used as a measure of chronic temperature exposure, encompassing at least the lifetime of every invertebrate species under investigation (and potentially multiple generations^6,35).

Invertebrates were collected from nine streams spanning a temperature gradient of 5–20 °C across the entire study system (Supplementary Figs. 1–2). The streams exhibit some differences in the annual variability of their thermal regimes, but there are examples of both cold and warm streams that have high (IS12 and IS2) and low (IS13 and IS8) variability throughout the year. Our main finding is also robust to the inclusion of stream temperature variability as a random effect in our modelling framework (Supplementary Table 4; Supplementary Fig. 7). Note that we present temperature data from 15 streams in Supplementary Fig. 2, but it was not logistically feasible to study acute thermal responses of invertebrates collected from all of them, thus we focused on a subset of nine streams that best spanned the temperature gradient. The remaining six streams were included in other studies from the system, quantifying the biomass of all the constituent species¹⁷, describing food web structure¹⁸, and measuring whole-stream respiration¹⁵ (described in detail below).

Individual organisms were stored in containers within their ‘home stream’ until the end of each collection day, when they were transported within 1 h to the University of Iceland and then transferred into 2 L aquaria filled with water from the main river in Hengill, the Hengladalsá. The water was passed through a 125 µm sieve to ensure no organisms or filamentous algae entered the aquaria, and thus limiting the potential food available to the study organisms. The aquaria were continuously aerated in temperature-controlled chambers set to the home stream temperature of the organisms during sampling, which were maintained without food for at least 24 h to standardise their digestive state prior to metabolic measurements³⁶. While we did not observe any cannibalism or organisms feeding on dead bodies in the laboratory, we cannot rule out the possibility that organisms fed on fine algal or detrital particles in the water, thus increasing variability in our metabolic measurements due to differences in digestive state.

Quantifying metabolic rates

Experiments were carried out to determine the effects of body mass, acute temperature exposure (5, 10, 15, 20 and 25 °C), and chronic temperature exposure (i.e., average stream temperature) on oxygen consumption rates as a measure of metabolic rate^3,12. Before each experiment, individual organisms were confined in glass chambers in a temperature-controlled water bath and slowly adjusted to the (acute) experimental temperature over a 15 min period to avoid a shock response. Glass chambers ranged in volume from 0.8–5 ml and scaled with the size of the organism. The glass chambers were filled with water from the Hengladalsá, which was filtered through a 0.45 µm Whatman membrane after aeration to 100% oxygen saturation. A magnetic stir bar was placed at the bottom of each chamber and separated from the organism by a mesh screen. In each experiment, one individual organism was placed in each of seven chambers and the eighth chamber was used as an animal-free control to correct for potential sensor drift. The chambers were sealed with gas-tight stoppers after the 15 min acclimatisation period, ensuring there was no headspace or air bubbles.

Oxygen consumption by individual organisms was measured using an oxygen microelectrode (MicroRespiration, Unisense, Denmark), fitted through a capillary in the gas-tight stopper of each chamber³⁷. A total of 330 s measurement periods were recorded for each individual, where dissolved oxygen was measured every second. Oxygen consumption rate was calculated as the slope of the linear regression through all the data points from a single chamber, corrected for differences in chamber volume and the background rate measured from the control chamber (which was never >5% of the measured metabolic rates). We converted the units of this rate (µmol O₂ h⁻¹) to energetic equivalents (J h⁻¹) using atomic weight (1 mol O₂ = 31.9988 g), density (1.429 g L⁻¹), and a standard conversion³⁸ (1 ml O₂ = 20.1 J). Organisms generally exhibited some activity during experiments, thus these measurements can be classified as routine metabolic rates¹², which are more reflective of energy expenditure in field conditions. Nevertheless, activity levels were minimal due to the space constraints of the chambers (volume equal to 5–100 times the mass of the measured organism), indicating that the measured rates were likely to be closer to resting metabolic rates. Oxygen concentrations were never allowed to decline below 70% to minimise stress and avoid oxygen limitation. The system was cleaned with bleach at the end of each measurement day to avoid accumulation of microbial organisms on the insides of glass chambers and the water bath. In total, oxygen consumption rates were measured for 1819 individuals, none of which were ever reused in another experiment, thus every data point in the analysis corresponds to a single new individual (see below for details of how this dataset was curated to the final analysed subset of 1359 individuals based on quality-control procedures).

Following each experiment, individuals were preserved in 70% ethanol and later identified to species level under a dissecting microscope, except for Chironomidae, which were identified by examining head capsules under a compound microscope³⁹. A linear dimension was precisely measured for every individual using an eyepiece graticule and converted to dry body mass using established length-weight relationships (Supplementary Table 1).

Statistical analysis

All statistical analyses were conducted in R 4.0.2 (see the Supplementary Note for full details of statistical R code). According to the Metabolic Theory of Ecology³ (MTE), metabolic rate, I, depends on body mass and temperature as:

where I₀ is the intercept, M is dry body mass (mg), b is an allometric exponent, E_A is the activation energy (eV), and T_A is a standardised Arrhenius temperature:

Here, T_acute is an acute temperature exposure (K), T₀ sets the intercept of the relationship at 283.15 K (i.e., 10 °C), and k is the Boltzmann constant (8.618 × 10⁻⁵ eV K⁻¹). We performed a multiple linear regression (‘lm’ function in the ‘stats’ package) on the natural logarithm of Eq. (1) to explore the main effects of temperature and body mass on the metabolic rate of each population (i.e., species × stream combination) in our dataset³. Following these analyses, we excluded populations where n < 10 individuals, r² < 0.5, and p > 0.05 for any term in the model (see Supplementary Table 5, Supplementary Figs. 8–9). This excluded any poor quality species-level data and resulted in 1359 individuals from 44 populations for further analysis. Note that we find the same overall conclusion if we analyse the entire dataset (Supplementary Table 6, Supplementary Fig. 10).

To determine whether chronic temperature exposure alters the size- and acute temperature-dependence of metabolic rate, we added a term for chronic temperature exposure to Eq. 1. We began our analysis by considering the natural logarithm of all possible combinations of the main and interactive effects in this model:

Here, T_C is a standardised Arrhenius temperature with T_chronic as a chronic temperature exposure (K) substituted for T_acute in Eq. (2). To determine the optimal random effects structure for this model, we compared a generalised least squares model of Eq. 3 with linear mixed-effects models (‘gls’ and ‘lme’ functions in the ‘nlme’ package) containing all possible subsets of the following random effects structure⁴⁰:

Here, we are accounting for the possibility that metabolic rate could be different for each species (i.e., a random intercept) and that the effect of body mass, acute temperature exposure, or chronic temperature exposure on metabolic rate could also be different for each species (i.e., random slopes).

The full random structure (Eq. 4) was identified as the best model using Akaike Information Criterion (ΔAIC > 31.2; see Supplementary Table 7). We used this random structure in subsequent analyses, set ‘method = ‘ML’’ in the ‘lme’ function, and performed AIC comparison on all possible combinations of the fixed-effect structure⁴⁰ (i.e., Equation 3). The optimal model was identified as follows:

Note that while the model with an additional interaction between ln(M) and T_A performed similarly (ΔAIC = 0.2; see Supplementary Table 8), that term was not significant (t = −1.645; p = 0.1002). We set ‘method = ‘REML’’ before extracting model summaries and partial residuals from the best-fitting model⁴⁰. Note that the models were always fitted to the raw metabolic rate data, with residuals only extracted for a visual representation of the best-fitting models, excluding the noise explained by the random effect of species identity (see R code in the Supplementary Note).

Exploration of spatial autocorrelation

A Mantel test (‘mantel’ function in the ‘vegan’ R package) was used to test for spatial autocorrelation in the temperature gradient, by comparing pairwise temperature difference between streams to the pairwise distance between streams. Pairwise distances were calculated from GPS coordinates taken at the confluence of each stream with the main river and the ‘earth.dist’ function in the ‘fossil’ R package. This analysis revealed no significant relationship between pairwise temperature and pairwise distance between sites (Mantel r = −0.1293, p = 0.780).

In addition, we explored for spatial autocorrelation in the residuals of our optimal model (Table 1a) by generating an empirical semivariogram cloud, illustrating the squared difference between all pairwise residual data points as a function of the distance between the two points. We also calculated Moran’s I as a measure of spatial autocorrelation in the model residuals. The semivariogram indicated no clear patterns in the residuals as a function of the distance between data points (Supplementary Fig. 3) and there was no statistical evidence for spatial autocorrelation in the model residuals (Moran’s I = 0.1187, p = 0.453).

Exploration of phylogenetic structure

To examine the influence of evolutionary relatedness on metabolic rate measurements, we reconstructed a time-calibrated phylogeny of the 16 species in our final dataset (Supplementary Table 1). To this end, we combined: (i) nucleotide sequences of the 5′ region of the cytochrome c oxidase subunit I gene (COI-5P) from the Barcode of Life Data System database⁴¹; (ii) tree topology information from the Open Tree of Life⁴² (OTL; v. 13.4); and (iii) previously reported divergence time estimates between pairs of genera from the TimeTree database⁴³. More precisely, we were able to obtain COI-5P nucleotide sequences for 15 out of 16 species (Supplementary Table 2), which we aligned using the G-INS-i algorithm of MAFFT⁴⁴ (v. 7.490). To constrain the topology of our phylogeny based on the results of previous studies, we queried the OTL via the ‘rotl’ R package⁴⁵ (v. 3.0.11). This yielded topological information for all 16 species. Finally, we manually queried the TimeTree database to obtain node age estimates. We only used three such estimates that (a) were based on more than five previous studies and (b) did not force any tree branches to have a length of zero.

We next used MrBayes⁴⁶ (v. 3.2.7a) to obtain a time-calibrated phylogeny based on the sequence alignment, the OTL topology, and the node ages from TimeTree. For this, we first determined the most appropriate nucleotide substitution model using ModelTest-NG⁴⁷ (v. 0.1.7). This was the General Time-Reversible model with Gamma-distributed rate variation across sites and a proportion of invariant sites. To allow branches of the phylogeny to differ in their rate of sequence evolution, we specified the Independent Gamma Rates model⁴⁸ and used a normal distribution with a mean of 0.00003 and a standard deviation of 0.00001 as the prior for the mean clock rate. Finally, we executed four MrBayes runs with two chains per run for 100 million generations, sampling from the posterior distribution every 500 generations. Samples from the first ten million generations were treated as burn-in and were discarded. We examined the remaining samples to ensure that the four MrBayes runs had converged on statistically indistinguishable posterior distributions (i.e., all potential scale reduction factor values were below 1.1) and the parameter space was sufficiently explored (i.e., all effective sample size values were higher than 200). We summarised the sampled trees into a single time-calibrated phylogeny by calculating the median age estimate for each node (Supplementary Fig. 4).

To investigate the influence of evolutionary and acclimatory processes on metabolic rate, we first estimated the phylogenetic heritability of metabolic rate, i.e., the extent to which closely related species have more similar trait values than species chosen at random⁴⁹. This metric takes values from 0 (trait values are independent of the phylogeny) to 1 (trait values evolve similarly to a random walk in the parameter space), with intermediate values indicating deviations from a pure random walk. To estimate phylogenetic heritability, we fitted a generalised linear mixed-effects model using the ‘MCMCglmm’ R package⁵⁰ (v. 2.32). We set the natural logarithm of metabolic rate as the response variable and only an intercept as a fixed effect. We also specified a phylogenetic species-level random effect on the intercept, using the phylogenetic variance-covariance matrix obtained from our time-calibrated phylogeny. We used the default (normal) prior for the fixed effect, an uninformative Cauchy prior for the random effect, and an uninformative inverse Gamma prior for the residual variance. We then executed four independent runs for 500,000 MCMC generations each, with parameter samples being obtained every 50 generations after the first 50,000. We verified that sufficient convergence was reached, based on potential scale reduction factor and effective sample size values, as described earlier. Phylogenetic heritability was calculated as the ratio of the variance captured by the species-level random effect to the sum of the random and residual variances. The mean posterior phylogenetic heritability estimate of the natural logarithm of metabolic rate was 0.48. This means that nearly half (48%) of the variation can be explained by the evolution of metabolic rate along the phylogeny (Supplementary Fig. 4), with the other half arising from other sources including (but not necessarily limited to) acclimation and measurement error.

To describe the remaining unexplained variation, we fitted a series of models using MCMCglmm in R with all possible combinations of log body mass, acute temperature exposure, and chronic temperature exposure (fixed effects, as in Eq. (3) of the main text) and species-level random effects on the intercept and slopes (as in Eq. (4) of the main text). Furthermore, we specified both phylogenetic and non-phylogenetic variants of each model to understand if such a correction is warranted when the fixed effects are included. We determined the most appropriate model based on the Deviance Information Criterion⁵¹ (DIC). The optimal model (ΔDIC > 19; Supplementary Table 3; Supplementary Fig. 5) was found to include the full random effects structure (Eq. 4), the main effects of log body mass, acute temperature exposure, and chronic temperature exposure, the interaction between log body mass and chronic temperature exposure, and the interaction between acute temperature exposure and chronic temperature exposure (as for Eq. 5 in the main text), i.e., the same optimal model as that containing only species-level, rather than phylogenetic, information (Table 1a; Fig. 1). We calculated the marginal and conditional coefficients of determination to report the amounts of variance explained by the fixed and random effects, or left unexplained⁵². We found that the unexplained variation dropped from 52% to 8%, indicating that metabolic rate is strongly influenced by acclimatory processes in addition to evolutionary processes (see above).

It should be noted, however, that a definitive empirical quantification of the relative strength of evolutionary and acclimatory processes would require population genetics (to determine evolutionary divergent populations among streams), transcriptomics (to identify the expression of genes associated with thermal adaptation), and exhaustive common garden experiments (to disentangle acclimation from adaptation in all populations). Such an undertaking was logistically unfeasible in this study, but should be a focus for follow-up research on this topic.

Modelling ecosystem-level energy fluxes

We used a recently proposed approach for inferring energy fluxes through trophic links²⁵ to predict the effects of climate warming on ecosystem-level energy fluxes. We began by assuming that each stream ecosystem is at energetic steady state, i.e., for all n consumer species in the system:

where G_i and L_i are the energy gain and loss rates [J h⁻¹], respectively, of the ith species in that stream. All basal species are implicitly assumed to be at energy balance. The two terms in Eq. (6) can be specified in a general way as

Here, for the ith species, R_i and C_i are the sets of its resource and consumer species respectively, and Z_i is its population-level energy loss rate stemming from mortality and metabolic expenditure on various activities realised over the timescale of the system’s dynamics. For the jth species feeding on the ith species, F_ij is the maximum population-level feeding rate, e_ij is the assimilation efficiency (expressed as a proportion), and w_ij is the consumer’s preference for that species (all preferences for a given consumer sum to 1). Thus, the effective flux through a trophic link is ({e}_{{ki}}{w}_{{ki}}{F}_{{ki}}). Next, assuming the energy balance condition in Eq. 6 holds for all species, there are n linear equations (corresponding to the n consumer species) of the form:

which can be solved iteratively to obtain the unknown fluxes ({F}_{{ij},{i}ne j}) of all consumer species, provided all the Z_i’s, e_ij’s, and w_ij’s are known.

For this, we used the ‘fluxing’ function in the ‘fluxweb’ R package, parameterised with: (1) binary predation matrices for 14 stream food webs, characterised by 49,324 directly observed feeding interactions¹⁸; (2) biomasses for every species in each food web, characterised by 13,185 individual body mass measurements¹⁷; (3) assimilation efficiencies (e_ij’s) based on an established temperature-dependence and resource type (i.e., plant, detritus, or invertebrate)⁵³; (4) preferences (w_ij’s) depending on resource biomasses; and (5) metabolic rates estimated using Eqs. (1) and (5) (assuming that I approximates Z). We treated T_A in Eqs. (1) and (5) as the short-term temperature of the streams during food web sampling^17,18 and T_C in Eq. (5) as the long-term average temperature of the streams measured over the current study (Supplementary Fig. 2). It is important to note that the energy balance assumption (Eq. 6) implies that Z_i in Eq. (8) is a combination of basal, routine, and active metabolic rates, stemming from the combination of activities realised over the timescale of the system’s dynamics. Therefore, our use of routine metabolic rate I is an underestimate of Z, which in turn means that the fluxes (which must balance the losses) are an underestimate.

Biomass and food web data were sampled in August 2008, with extensive protocols described in previous publications^17,18. Briefly, this involved three stone scrapes per stream for benthic diatoms, five Surber samples per stream for macroinvertebrates, and three-run depletion electrofishing for fish. All individuals in the samples were identified to species level where possible and counted. Linear dimensions were measured for at least ten individuals of each species in each stream, with body masses estimated from length-weight relationships¹⁷. The population biomass of each species in each stream was calculated as the total abundance [individuals m⁻²] multiplied by the mean body mass [mg dry weight]. Food web links were largely assembled from gut content analysis of individual organisms collected from the streams (>87% of all links in the database), but additional links were added from the literature when yield-effort curves indicated that the diet of a consumer species was incomplete¹⁸.

Validation of the ecosystem flux model using field data

To test whether our model of energy fluxes through trophic links was empirically meaningful, we calculated the sum of all energy fluxes through each stream food web to get the total energy flux, F (i.e., the sum of all ({e}_{{ki}}{w}_{{ki}}{F}_{{ki}})’s in Eq. 7). This quantity is a measure of multitrophic functioning and is expected to be positively correlated with the total respiration of each stream²⁵. To evaluate this, we compared F to whole-ecosystem respiration rates measured in the same study streams¹⁵. The ecosystem respiration estimates were based on a modified open-system oxygen change method using two stations corrected for lateral inflows^54,55. Essentially, this was an in-stream mass balance of oxygen inflows and outflows along stream reaches (17–51 m long). Oxygen concentrations were measured during 24- to 48 h periods from 6th to 16th August 2008, i.e., the exact same time period during which biomass and food web data were sampled to parameterise the energy flux model¹⁵. Dissolved oxygen concentrations were measured every minute with optic oxygen sensors (TROLL9500 Professional, In-Situ Inc. and Universal Controller SC100, Hach Lange GMBF). Hourly ecosystem respiration was calculated from the net metabolism at night, i.e., when no primary production occurs due to lack of sunlight.

Modelling the consequences of metabolic plasticity for global warming impacts on ecosystem-level energy flux

In addition to total energy flux, F, we also calculated a modified total energy flux, F^*, for each food web after considering a global warming scenario, where we added 2 °C to T_A in Eq. (1) and to both T_A and T_C in Eq. (5). We calculated the change in total energy flux as a result of the global warming scenario as ΔF = F^* – F. We tested whether the (statistically optimal) model with metabolic plasticity (Eq. 5) predicted a greater ΔF across the 14 empirical stream food webs from the Hengill system than the model without metabolic plasticity using paired Wilcoxon tests (since the data did not conform to homogeneity of variance). To determine whether our results were consistent for all major trophic groupings in the system, we repeated the analysis after calculating the change in energy flux to herbivores (ΔF_H = F_H^* – F_H), detritivores (ΔF_D = F_D^* – F_D), and predators (ΔF_P = F_P^* – F_P) in each stream.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com