Study system & organisms
The study was conducted in the Hengill valley, Iceland12-17 (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 (Fig. S1). The streams have been heated in this way for centuries30 and are otherwise similar in their physical and chemical properties12,17, providing an ideal space-for-time-substitution in which to measure species responses after chronic exposure to different temperatures6,31. Fieldwork was performed from 30th April 2015 to 5th July 2017, with some supplementary experiments carried out from 3rd to 17th June 2018. Stream temperatures were logged every four hours using Maxim Integrated DS1921G Thermochron iButtons submerged in each stream (Fig. S2). 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 generations6,32).
Invertebrates were collected from nine streams spanning a temperature gradient of 5‑20 °C across the entire study system (Fig. S1-S2). Organisms were stored in containers within their ‘home stream’ until the end of each collection day, when they were transported within one hour to the University of Iceland and then transferred into 2 L aquaria filled with water from the main river in Hengill, the Hengladalsá. The aquaria were continuously aerated in temperature-controlled chambers set to the home-stream temperature of the organisms, which were maintained without food for at least 24 hours to standardise their digestive state prior to metabolic measurements33.
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 rate3,11. Before each experiment, individual organisms were confined in glass chambers and slowly adjusted to the (acute) experimental temperature over a 15-minute period to avoid a shock response. 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 after the acclimatisation period with gas-tight stoppers, 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 chamber34. Three 30-second 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. We converted the units of this rate (µmol O2 h-1) to energetic equivalents (J h-1) using atomic weight (1 mol O2 = 31.9988 g), density (1.429 g L‑1), and a standard conversion35 (1 ml O2 = 20.1 J). Organisms generally exhibited normal activity during experiments, thus these measurements can be classified as routine metabolic rates11, which reflect normal energy expenditure in field conditions, excluding the rates associated with activities such as mating, predator attack, or prey escape. Oxygen concentrations were never allowed to decline below 70% to minimise stress and avoid oxygen limitation. In total, oxygen consumption rates were measured for 1,819 individuals, none of which were ever reused in another experiment, thus every data point in the analysis corresponds to a single new individual.
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 microscope36. A linear dimension was precisely measured for every individual using an eyepiece graticule and converted to dry body mass using established length-weight relationships (Table S1).
Statistical analysis
All statistical analyses were conducted in R 4.0.2 (see Supplementary Methods for full details of statistical R code). According to the Metabolic Theory of Ecology3 (MTE), metabolic rate, I, depends on body mass and temperature as:
$$I={I}_{0}{M}^{b}{e}^{{E}_{A}{T}_{A}}$$
1
where I0 is the intercept, M is dry body mass (mg), b is an allometric exponent, EA is the activation energy (eV), and TA is a standardised Arrhenius temperature:
$${T}_{A}=\frac{{T}_{acute}-{T}_{0}}{k{T}_{acute}{T}_{0}}$$
2
Here, Tacute is an acute temperature exposure (K), T0 sets the intercept of the relationship at 283.15 K (i.e., 10°C), and k is the Boltzmann constant (8.618 × 10− 5 eV K1). 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 dataset3. Following these analyses, we excluded populations where n < 10 individuals, r2 < 0.5, and p > 0.05 for any term in the model (see Table S2, Fig. S4-S5). This excluded any poor quality species-level data and resulted in 1,359 individuals from 44 populations for further analysis. Note that we find the same overall conclusion if we analyse the entire dataset (Table S3, Fig. S6).
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:
$${ln}I={ln}{I}_{0}+b{ln}M+{E}_{A}{T}_{A}+{E}_{C}{T}_{C}+{b}_{A}{ln}M{T}_{A}+{b}_{C}{ln}M{T}_{C}+{E}_{AC}{T}_{A}{T}_{C}+{b}_{AC}{ln}M{T}_{A}{T}_{C}$$
3
Here, TC is a standardised Arrhenius temperature with Tchronic as a chronic temperature exposure (K) substituted for Tacute 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 structure37:
$$random=\tilde1+{ln}M+{T}_{A}+{T}_{C}|species$$
4
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 Table S4). 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 structure37 (i.e., Eq. 3). The optimal model was identified as follows:
$${ln}I={ln}{I}_{0}+b{ln}M+{E}_{A}{T}_{A}+{E}_{C}{T}_{C}+{b}_{C}{ln}M{T}_{C}+{E}_{AC}{T}_{A}{T}_{C}$$
5
Note that while the model with an additional interaction between ln(M) and TA performed similarly (ΔAIC = 0.2; see Table S5), 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 model37. 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 Supplementary Methods).
Modelling ecosystem-level energy fluxes
We used a recently proposed approach for inferring energy fluxes through trophic links24 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:
$${G}_{i} = {L}_{i}, i = 1, 2, \dots , n$$
6
where Gi and Li are the energy gain and loss rates [J h-1], 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:
\({G}_{i}=\sum _{k \in {\text{R}}_{i}}{e}_{ki}{w}_{ki}{F}_{ki}\) , and (7)
$${L}_{i}={Z}_{i}+\sum _{j \in {\text{C}}_{i}}{w}_{ij}{F}_{ij}$$
8
Here, for the ith species, Ri and Ci are the sets of its resource and consumer species respectively, and Zi 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, Fij is the maximum population-level feeding rate, eij is the assimilation efficiency (expressed as a proportion), and wij 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:
$${G}_{i}-{L}_{i}=\sum _{k \in {\text{R}\text{e}\text{s}}_{i}}{e}_{ki}{w}_{ki}{F}_{ki}-\left({Z}_{i}+\sum _{j \in {\text{C}\text{o}\text{n}}_{i}}{w}_{ij}{F}_{ij}\right)=0$$
9
which can be solved iteratively to obtain the unknown fluxes \({F}_{ij, i\ne j}\) of all consumer species, provided all the Zi’s, eij’s, and wij’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 interactions17; (2) biomasses for every species in each food web, characterised by 13,185 individual body mass measurements16; (3) assimilation efficiencies (eij’s) based on an established temperature-dependence and resource type (i.e., plant, detritus, or invertebrate)38; (4) preferences (wij’s) depending on resource biomasses; and (5) metabolic rates estimated using Equations 1 and 5 (assuming that I approximates Z). We treated TA in Equations 1 and 5 as the short-term temperature of the streams during food web sampling16,17 and TC in Eq. 5 as the long-term average temperature of the streams measured over the current study (Fig. S2). It is important to note that the energy balance assumption (Eq. 6) implies that Zi 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 publications16,17. 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 relationships16. The population biomass of each species in each stream was calculated as the total abundance [individuals m-2] 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 incomplete17.
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 stream24. To evaluate this, we compared F to whole-ecosystem respiration rates measured in the same study streams14. The ecosystem respiration estimates were based on a modified open-system oxygen change method using two stations corrected for lateral inflows39,40. 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-hour 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 model14. 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 TA in Eq. 1 and to both TA and TC 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 (ΔFH = FH* – FH), detritivores (ΔFD = FD* – FD), and predators (ΔFP = FP* – FP) in each stream.
Data and Code Availability:
The data and R code that support the findings of this study are available from the corresponding author upon reasonable request.
Additional References:
30 Arnason, B., Theodorsson, P., Björnsson, S. & Saemundsson, K. Hengill, a high temperature thermal area in Iceland. Bulletin of Volcanology 33, 245–259 (1969).
31 O'Gorman, E. J. et al. Climate change and geothermal ecosystems: natural laboratories, sentinel systems, and future refugia. Global Change Biology 20, 3291–3299 (2014).
32 Johansson, M. P., Quintela, M. & Laurila, A. Genetic divergence and isolation by thermal environment in geothermal populations of an aquatic invertebrate. Journal of evolutionary biology 29, 1701–1712 (2016).
33 Vucic-Pestic, O., Ehnes, R. B., Rall, B. C. & Brose, U. Warming up the system: higher predator feeding rates but lower energetic efficiencies. Global Change Biology 17, 1301–1310 (2011).
34 Brodersen, K. P., Pedersen, O., Walker, I. R. & Jensen, M. T. Respiration of midges (Diptera; Chironomidae) in British Columbian lakes: oxy-regulation, temperature and their role as palaeo-indicators. Freshwater Biology 53, 593–602 (2008).
35 Peters, R. H. The ecological implications of body size. (Cambridge University Press, 1983).
36 Brooks, S. J., Langdon, P. G. & Heiri, O. The identification and use of Palaearctic Chironomidae larvae in palaeoecology. (Quaternary Research Association, 2007).
37 Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. in Mixed effects models and extensions in ecology with R 101–142 (Springer, 2009).
38 Lang, B., Ehnes, R. B., Brose, U. & Rall, B. C. Temperature and consumer type dependencies of energy flows in natural communities. Oikos 126, 1717–1725 (2017).
39 Demars, B. O. L., Manson, J. R., Olafsson, J. S., Gislason, G. M. & Friberg, N. Stream hydraulics and temperature determine the metabolism of geothermal Icelandic streams. Knowledge and Management of Aquatic Ecosystems 402, 1–17 (2011).
40 Odum, H. T. Primary production in flowing waters. Limnology and Oceanography 1, 102–117 (1956).