Anaerobic metabolism at thermal extremes: a metabolomic test of the oxygen limitation hypothesis in an aquatic insect.

Thermal limits in ectotherms may arise through a mismatch between supply and demand of oxygen. At higher temperatures, the ability of their cardiac and ventilatory activities to supply oxygen becomes insufficient to meet their elevated oxygen demand. Consequently, higher levels of oxygen in the environment are predicted to enhance tolerance of heat, whereas reductions in oxygen are expected to reduce thermal limits. Here, we extend previous research on thermal limits and oxygen limitation in aquatic insect larvae and directly test the hypothesis of increased anaerobic metabolism and lower energy status at thermal extremes. We quantified metabolite profiles in stonefly nymphs under varying temperatures and oxygen levels. Under normoxia, the concept of oxygen limitation applies to the insects studied. Shifts in the metabolome of heat-stressed stonefly nymphs clearly indicate the onset of anaerobic metabolism (e.g., accumulation of lactate, acetate, and alanine), a perturbation of the tricarboxylic acid cycle (e.g., accumulation of succinate and malate), and a decrease in energy status (e.g., ATP), with corresponding decreases in their ability to survive heat stress. These shifts were more pronounced under hypoxic conditions, and negated by hyperoxia, which also improved heat tolerance. Perturbations of metabolic pathways in response to either heat stress or hypoxia were found to be somewhat similar but not identical. Under hypoxia, energy status was greatly compromised at thermal extremes, but energy shortage and anaerobic metabolism could not be conclusively identified as the sole cause underlying thermal limits under hyperoxia. Metabolomics proved useful for suggesting a range of possible mechanisms to explore in future investigations, such as the involvement of leaking membranes or free radicals. In doing so, metabolomics provided a more complete picture of changes in metabolism under hypoxia and heat stress.

Introduction

Thermal limits in ectotherms may arise through a mismatch between supply and demand of oxygen (Winterstein 1905; Pörtner 2001). As temperatures increase, there is an increase in both the environmental availability of oxygen and organismal demand for oxygen, with the latter tending to increase more rapidly (Verberk et al. 2011). As a result, the ratio of oxygen supply relative to demand decreases and shortages of oxygen may arise in warmer waters. Consequently, cardiac and ventilatory activities of ectotherms may be insufficient to meet their elevated demand for oxygen at the higher temperatures (Frederich and Pörtner 2000). The resulting whole-animal drop in aerobic scope induces a shift from aerobic to anaerobic metabolism to maintain energy status (Pörtner 2006).

Evidence for temperature-dependent oxygen limitation has not been forthcoming for terrestrial insects (Klok et al. 2004; Stevens et al. 2010); most of the evidence for oxygen limitation at thermal extremes to date comes from a variety of marine taxa (Pörtner 2001). In tracheated arthropods, oxygen limitation may be unlikely due to the efficiency and plasticity of tracheal systems in supplying oxygen directly to metabolically active tissues. Recently, Verberk and Bilton (2011) found that better oxygenation of the medium improves heat tolerance in the aquatic nymphs of the stonefly Dinocras cephalotes (Curtis 1827). These nymphs have closed tracheal systems and rely on diffusion across their tegument for gas exchange and are hence limited in their regulatory capacity. Although these findings support the idea that a mismatch between external oxygen supply and internal oxygen demand can set thermal limits in aquatic insects they do not conclusively show that thermal limits are caused by a loss of aerobic scope of the whole animal. This requires demonstrating the onset of anaerobic metabolic pathways at thermal extremes and a modulating effect of oxygen on the onset of these pathways. Here, we extend previous research on thermal limits and oxygen limitation in aquatic insect larvae and investigate oxygen limitation directly by looking at the onset of anaerobic metabolism.

Generally, anaerobic metabolism is not intensively studied in insects (Müller et al. 2012) and specific anaerobic end products are unknown for many insects, including this species. Therefore, rather than measuring a small set of metabolites that are typically associated with anaerobic metabolism (e.g., lactate, succinate, and alanine) but risk missing the actual end products, we chose to use an untargeted metabolomic approach. Metabolomics is the now well-established technology concerned with the study of naturally occurring, low-molecular-weight organic metabolites within a biological specimen. Environmental metabolomics can be defined simply as the application of this technology to characterize the interactions of organisms with their environment (Malmendal et al. 2006; Bundy et al. 2009; Colinet et al. 2012; Viant and Sommer 2013). This approach has considerable potential for characterizing the responses of aquatic organisms to natural and anthropogenic stressors (Viant 2007). Here, we used metabolomics in stonefly nymphs for the first time to discover (1) which metabolites are accumulating and if they indicate the onset of anaerobic metabolism and (2) whether the pattern in metabolites supports the hypothesis of increased anaerobic metabolism and lower energy status at thermal extremes.

Methods

Maintenance of animals

The life-cycle of the stonefly D. cephalotes lasts 3 years, and most of this time is spent as aquatic nymphs. Aquatic perlid nymphs such as D. cephalotes are among the largest European freshwater carnivorous invertebrates, feeding primarily at night (Malmqvist and Sjöström 1980). Nymphs pass through 15–21 instars, and growth is slower for males than females, resulting in pronounced sexual dimorphism (Frutiger 1987). Maximum growth of nymphs occurs in spring-summer. Adult flight and oviposition occur chiefly in May and June.

Aquatic nymphs were collected from the River Dart, Devon, UK at the same site in (early) spring in 2010 and (late) spring in 2011. Higher air temperatures and less rainfall were recorded for SW England in 2011 compared with 2010 for the period prior to animal collection (January–April). Nymphs were housed individually in separate chambers that were placed in a flow-through aquarium (10 l min−1), fed with artificial pond water (ASTM 1980), buffered, and diluted to reflect the pH and conductivity of the field site (pH 6.4–6.6, 70–150 µS cm−1). Nymphs were fed chironomid larvae and maintained in the laboratory at 10 ± 1°C in a 12 L:12 D regime. We could maintain nymphs under these conditions for extended periods of time (months), and all nymphs were acclimated for at least 7 days to laboratory conditions before trials were undertaken to measure endpoints used to assess critical temperatures. Nymphs could not be sexed reliably, and therefore, we could not account for differences between sexes, which may have increased the variance in our results.

Oxygen consumption

Oxygen consumption was measured for each larvae at 5, 10, and 15°C using closed glass respiration chambers of 67.5–68.5 ml. The chambers were immersed in a temperature controlled bath (±0.1°C) and stirred using submersible magnetic stirrers to ensure mixing of water. Respiration chambers were fitted with a fine nylon mesh forming a false bottom to prevent contact between the larva and the magnetic stirrer bar. Individuals were allowed to acclimate for 60 min before the chambers were closed and left for at least 60 min (60–140 min). Oxygen content was measured before closing the chambers and after the incubation period using an oxygen electrode (1302 Oxygen Electrode, Strathkelvin Instruments) connected to a calibrated meter (Oxygen Meter 929, Strathkelvin Instruments). On average, larvae depleted oxygen levels to 94% of the initial values (minimum 82%), and oxygen consumption was expressed as μg O2 h–1 g fresh weight–1.

Critical temperatures

To assess critical temperatures, we employed the same methods as previously described (Verberk and Bilton 2011; Verberk and Calosi 2012). Briefly, individual nymphs were placed in five parallel flow-through chambers (70 × 70 × 30 mm; flow rate 0.016 l s−1), and water was supplied to these chambers by gravity from a 25-l header tank after having passed through a tubular counter-current heat exchanger. Water in the header tank was of the same composition as that used to maintain animals and was bubbled with a mixture of 20% O2 and 80% N2, obtained using a gas-mixing pump (Wösthoff, Bochum, Germany). Individuals were left resting for 1 h at the equilibration temperature of 10°C, after which temperature in the experimental chambers was increased to 0.25°C min−1, using a Grant R5 water bath with a GP200 pump unit (Grant Instrument Ltd, Cambridge, UK), connected to the heat exchanger. Temperatures were logged using a HH806AU digital thermometer (Omega Engineering Inc., Stamford, CT).

Critical maximal temperatures (CTmax) were recorded as the point at which animals no longer showed any body movement or muscular spasms. Animals that were transferred at this point to fully oxygenated water of 10°C recovered (i.e., they survived and performed equally well in subsequent trials). At lower temperatures, larvae initiated in repeated swimming behavior (at about 29°C; interpreted as attempts to escape experimental conditions) and fell upon their backs (at about 30°C), and an onset of spasms was observed (at about 31°C). CTmax was assessed at normoxia, hypoxia, and hyperoxia. Different levels of oxygenation were achieved by changing the O2-N2 gas mixture obtained using the gas-mixing pump (Wösthoff, Bochum, Germany). The gas mixture was adjusted 10 min after placing the animals in the small flow-through chambers, to allow for gradual exposure to hypoxic and hyperoxic conditions during the 1 h resting period. During this period, oxygen levels in the outflow water from the chambers were measured repeatedly, to verify that the oxygen levels had stabilized. Trials to assess CTmax were performed in 2010 and 2011. In 2010, mild hypoxic and hyperoxic conditions were used (14 and 36 kPa, respectively). In 2011, more extreme hypoxic and hyperoxic conditions were used (5 and 60 kPa, respectively).

Design of the study of metabolomics

At the end of the trial, when animals reached the critical temperature, they were freeze-clamped in liquid nitrogen. These samples constituted three treatments, one for each level of oxygenation (hypoxia, normoxia, and hyperoxia). Two additional treatments comprised (1) animals freeze-clamped in liquid nitrogen before commencement of a heating trial (i.e., control animals at 10°C) and (2) animals subjected to a warming trial under hyperoxic conditions, but freeze-clamped in liquid nitrogen after reaching the critical temperature observed under hypoxic conditions (i.e., before reaching their critical temperature). All animals thus frozen were stored at −85°C until further analysis. So in total there are five treatments: The first three treatments were all at the thermal limit and are therefore expected to be at a temperature at which oxygen limitation has set in. In this sense, they could be said to have been standardized for this effect, even though the absolute temperatures varied in accordance with the level of oxygenation. In the last two treatments, animals were freeze-clamped in liquid nitrogen before reaching their thermal limit, with the control animals experiencing no heat stress, whereas the animals under hyperoxia in treatment 4 did experience the same heat stress as the first treatment, but (presumably) before oxygen limitation set in. These different conditions across these treatments allow three analyses that provide insight in the metabolites that are indicative of heat stress, oxygen limitation, or both. The first comparison is between nymphs at their critical temperature (treatments 1–3) and control animals (treatment 5). The metabolites differentiating these two groups provide a test of the hypothesis that thermal tolerance is limited by a shortage of oxygen and energy (as indicated by anaerobic end products and low ATP levels). The second comparison is between nymphs at their critical temperature but at different levels of oxygen (treatments 1–3). This allows a test of the hypothesis that metabolic profiles of animals at their critical temperature are comparable irrespective of the oxygen conditions as they should all be experiencing oxygen limitation. The third comparison between nymphs under hypoxia and hyperoxia, controlling for temperature (treatments 1 and 4), provides a test as to whether the increased heat tolerance under hyperoxia is associated with an alleviation or postponement of anaerobic metabolism.

Hypoxia at CTmax.

Normoxia at CTmax.

Hyperoxia at CTmax.

Hyperoxia at CThypoxia (treatment 1).

Control (acclimation at 10°C, 20 kPa).

Extraction of metabolites

Metabolites were extracted from each biological sample using a two-step methanol/water/chloroform protocol (Wu 2008). In brief, per milligram of tissue, 1 µl methanol (high-performance liquid chromatography [HPLC] grade) and 0.2 µl water (HPLC grade) were added to each sample (one freeze-clamped stonefly nymph) prior to homogenization, using a Precellys-24 ceramic bead-based homogenizer (Stretton Scientific Ltd, UK). Thereafter, per milligram of tissue, 0.5 µl water (HPLC grade) and 1 µl chloroform (pesticide analysis grade) were added to each sample, which was vortexed and centrifuged (4000 rcf, 10 min), yielding an upper (polar) and lower (non-polar) fraction for each sample. Each polar fraction was collected into a 1.5-ml Eppendorf tube, and the non-polar fraction into a 1.8-ml glass vial. Of the polar extract, 5 µl was taken for mass spectrometry and 50 µl for NMR spectroscopy. All samples were dried using a centrifugal concentrator (Thermo Savant, Holbrook, NY) and stored at −80°C prior to analysis.

FT-ICR mass spectrometry and spectral processing

Immediately prior to mass spectrometric analysis, polar fractions were resuspended in 40 µl methanol/water (4:1) with 20 mM ammonium acetate for negative ion analysis (the non-polar samples were not analyzed), and 5 µl each was taken to pool a quality control sample. For initial optimization of concentration, test samples were further diluted 1:1, 1:3, 1:7, and 1:15; following this, the main samples were not further diluted. Samples were centrifuged (10 min, 22,000 rcf, 4°C) to remove any particulate matter. Fourier transform ion cyclotron resonance (FT-ICR) direct infusion mass spectrometry was performed using an LTQ FT Ultra (Thermo Fisher Scientific, Bremen, Germany) equipped with a chip-based direct infusion nanoelectrospray ion source (Triversa, Advion Biosciences, Ithaca, NY). Each sample was analyzed in triplicate from a 96-well plate using the selected ion monitoring (SIM)-stitching method as adapted for the LTQ FT Ultra (Weber et al. 2011) from m/z 70 to 590 in both ion modes. Raw mass spectral data were processed using the SIM-stitching algorithm (Southam et al. 2007; Payne et al. 2009). Missing values were filled into the resulting data matrix using the KNN algorithm (Hrydziuszko and Viant 2012) in an in-house R script, whereas PQN normalization (Dieterle et al. 2006) was performed using again an in-house Matlab script. This normalized matrix was transformed using the generalized logarithm (Parsons et al. 2007) to stabilize the variance across the thousands of measured peaks prior to analysis with multivariate statistics (see below).

NMR spectroscopy and spectral processing

The dried polar extracts were resuspended in 650 µl sodium phosphate buffer in D2O (0.1 M, pH 7.0) containing 0.5 mM sodium 3-trimethylsilyl-2,2,3,3-d4-propionate, which serves as an internal chemical shift standard. The samples were centrifuged, and then 600 µl of each supernatant was transferred to a 5-mm NMR tube and analyzed immediately using a DRX-500 NMR spectrometer (Bruker Biospin, Coventry, UK), equipped with a cryoprobe and operated at 500.18 MHz (at 300 K). One-dimensional (1D) 1H NMR spectra were obtained as described by Hines et al. (2007). Spectra were obtained using excitation sculpting for water suppression (Hwang and Shaka 1995) and using an 8.4 μs (60°) pulse, 6 kHz spectral width, and a 2.5 s relaxation delay with water presaturation. A total of 64 transients were collected into 16,348 data points, requiring a 4.5-min acquisition time. Datasets were zero-filled to 32,768 points, before line-broadenings of 0.5 Hz were applied prior to Fourier transformation.

To maximize metabolite discrimination two-dimensional (2D) 1H J-resolved (JRES) NMR spectra were also acquired (Viant 2003), being processed according to Hines et al. (2007). 2D JRES spectra were acquired for each sample using 16 transients per increment, for 16 increments, which were collected into 16,000 data points with spectral widths of 6 kHz in F2 (chemical shift axis) and 50 kHz in F1 (spin-spin coupling constant axis). A 4.0-s relaxation delay was employed resulting in a total acquisition time of 24 min. Datasets were zero-filled in F1, the F2 dimension was then multiplied by a SEM window function using 0.5 Hz line broadening while the F1 dimension was multiplied by a sine-bell window function, all prior to Fourier transformation. 2D JRES spectra were tilted by 45°, symmetrized about F1, and calibrated using TopSpin (Bruker Biospin). Data were exported as the 1D skyline projections of 2D JRES spectra and converted to a format for multivariate statistical analyses using NMRLab (Günther et al. 2000).

Statistical analysis and identification of metabolites

Principal components analysis (PCA) was used initially to assess the overall metabolic differences between the sample groups in an unbiased manner, using the PLS_Toolbox (version 5.5.1, Eigenvector Research, Manson, WA) within Matlab (version 7.8; The MathsWorks, Natick, MA). Supervized multivariate analyses were performed using partial least-squares discriminant and regression approaches (PLS-DA and PLS-R, respectively) using the PLS Toolbox, with internal cross-validation (venetian blinds) and permutation testing (1000 tests each) using in-house Matlab scripts. Univariate statistical analyses were used to confirm the significance of changes in individual mass-spectral signals. Specifically, analysis of variance was conducted using in-house R scripts (with a false discovery rate of 5% to correct for multiple hypothesis testing) (Benjamini and Hochberg 1995). For FT-ICR mass-spectrometry data, identification of putative empirical formulae was achieved using MI-Pack software (Weber and Viant 2010), the KEGG database (http://www.genome.jp/kegg/download/), and calculations in Xcalibur 2.0.7 (Thermo Fisher Scientific). For NMR data, peaks were identified and integrated using in-house peak-picking algorithms (J. Byrne, C. Ludwig, J. M. Easton, S. He, and M. R. Viant, in preparation).

Results

Ecophysiology of stonefly nymphs

Oxygen consumption by stonefly nymphs increased with increasing temperatures and body size, with temperature effects being stronger in larger bodied animals (F1,207 = 27.16; P < 0.0001; Fig. 1). Overall, oxygen consumption was comparable for nymphs collected in 2010 and 2011 (F1,207 = 1.51; P = 0.221; Supplementary Table S1), although slightly stronger responses to temperature were found for nymphs collected in 2011 (F1,207 = 5.537; P = 0.0196; Supplementary Table S1). The critical thermal maxima (CTmax) observed by behavioral responses was consistently influenced by the level of oxygen (F1,77 = 31.22; P < 0.0001), suggesting oxygen limitation at thermal extremes (Fig. 2 and Supplementary Table S2). Interannual differences in CTmax were found (F1,77 = 9.88; P = 0.0024), with nymphs collected in 2011 displaying lower CTmax at normoxia (36.1°C in 2010 versus 35.0°C in 2011). In addition, there was a modest effect of body size, with larger sized individuals displaying lower heat tolerance (β = −3.53, t1,76 = −2.052, P = 0.044). Nevertheless, treatment effects of the level of oxygen predominated; relative to normoxia, hypoxia decreased CTmax by 2.9°C (at 14 kPa) and 7.8°C (at 5 kPa), whereas hyperoxia increased CTmax by 1.5°C (at 36 kPa) and 1.9°C (at 60 kPa).

Fig. 1

Oxygen consumption of stonefly nymphs at 5, 10, and 15°C shown separately for small (<260 mg) and large nymphs (>260 mg). Values in brackets refer to the numbers of measurements. Note that mass specific rates of oxygen consumption are shown; absolute rates of oxygen consumption were higher for larger animals.

Fig. 2

Critical thermal maxima (CTmax) at different levels of oxygen for stonefly nymphs collected in 2010 and 2011. Values in brackets refer to the numbers of individuals measured. Differences relative to normoxia are indicated. Note that the hypoxia and hyperoxia levels differ for 2010 and 2011.

Metabolomics of stonefly nymphs

Peaks from the FT-ICR mass spectra of homogenates of whole organisms were annotated with empirical formula(e) and putative metabolite names using the MI-Pack software. Of the >3000 peaks detected, 99.6% were assigned to one or more empirical formula(e), representing the most extensive putative annotation of the stonefly metabolome to date. Transformation mapping in conjunction with the KEGG database was then used to attempt to assign putative metabolite names to these empirical formulae (Weber and Viant 2010). In total, 196 peaks could be assigned to one or more compound names listed in the KEGG database. The metabolomics data presented here are based primarily on those derived from the FT-ICR mass spectra, and where NMR data were available for the same metabolites, the measurements were compared across the two analytical approaches. Both methods yielded highly similar results (Supplementary Fig. S1), which is interesting from a methodological point of view. In addition, the NMR data provided greater confidence in the identification of the metabolites, in particular when a signal in the mass spectra could originate from multiple putative metabolite names having the same empirical formula. For example, identification of lactate, succinate, fumarate, tyrosine, and alanine (Supplementary Table S3) is based on the fact that these metabolites were also identified in the NMR spectra.

Preliminary data analysis using PCA showed that differences in metabolic profiles between the nymphs collected in 2010 and 2011 were subordinate to the differences observed across the different treatment (Supplementary Fig. S2). Therefore, the remainder of the analyses focuses on comparing the various treatments to elucidate metabolic responses to temperature and oxygen.

Comparing the metabolites of nymphs at 10°C with those at their thermal limit

To determine whether our metabolomics approach could distinguish between the metabolome of nymphs at their critical temperature (treatments 1–3) and that of control animals (treatment 5), we conducted supervized multivariate classification (PLS-DA). The optimal PLS model comprised five latent variables (LVs), based on the minimization of the cross-validated classification errors, and accounted for 45.6% of the total variance in the metabolic dataset. The metabolic profile of stonefly nymphs reaching the critical thermal maxima (CTmax; >33°C), and control animals kept at 10°C were clearly distinguishable (Fig. 3). The cross-validated classification error was 15.9% with a significance of P < 0.001 when compared with 1000 random permutations. Furthermore, the metabolites important in separating these two groups did, to some extent, also differentiate within nymphs at their thermal limit according to oxygen level; nymphs under hyperoxia tended to be more similar in metabolic profile to control animals, and nymphs under hypoxia were least similar to control animals (Supplementary Table S3). The metabolites differentiating between both groups are associated with the generation of energy, involving sugar metabolism and anaerobic metabolism (e.g., lactate and succinate; see Fig. 4 and Supplementary Table S3). In addition, ATP levels were also lower in animals at their thermal limit. Uridine levels were elevated in tandem with lactate (Spearman rank correlation across all individuals: r = 0.734, P < 0.001). Furthermore, ATP levels were strongly related to levels of l-arginine phosphate (Spearman rank correlation across all individuals: r = 0.405, P < 0.001).

Fig. 3

PLS-DA scores plot, showing a clear separation in metabolic profile between nymphs at 10°C and those at their thermal limit. All nymphs at their thermal limit (treatments 1–3) were treated as one class for the PLS-DA calculation (cross-validated class error rate of 15.9%, P < 0.001 [1000 random permutations]). The best model had five LVs in total, but for clarity only the first two axes are shown. The percentage of data variance explained by each latent variable is indicated. Ellipses are drawn by eye to improve visibility.

Fig. 4

Levels of individual metabolites relative to control animals (expressed as a percentage ± standard error) for treatments 1–4. All metabolite data are based on FT-ICR mass spectrometry except for isocitrate, oxaloacetate, and acetate, which are based on NMR data. Nonoverlapping 95% confidence intervals are indicated by an asterisk. Metabolites are color coded, indicating metabolites involved in glycolysis (blue), the TCA cycle (red), energetic and sugars (green), amino acid and nucleotide metabolism (yellow), and antioxidant responses (purple). Note the difference in scale on the x axis for positive values (+500%) and negative values (−100%).

Comparing the metabolism of nymphs reaching their thermal limit at different levels of oxygen

To assess whether there were consistent differences in the metabolite profile of animals reaching their critical thermal maxima but at different oxygen levels (treatments 1–3), we conducted supervized multivariate regression (PLS-R). Perturbations in the metabolic pathways were highly predictive of the oxygen conditions that nymphs experienced during the experiments. The model with four LVs accounted for 43.3% of the variance in the total metabolic dataset (x-block). Oxygen values (y) predicted from the metabolic profiles accurately matched those used in the experiment (x) (Fig. 5; partial least squares regression; y = 0.981x; R2 = 0.944, Q2 [cross-validated] = 0.605, P < 0.001 [1000 random permutations]). This suggests that the metabolic responses to heat stress in the various oxygen treatments were distinct. The metabolites differentiating between oxygen conditions were again associated with the generation of energy, with anaerobic metabolism being more pronounced under hypoxia (treatment 1) (e.g., lactate and succinate) (see Fig. 4 and Supplementary Table S3). In addition, metabolism of amino acids was perturbed (e.g., alanine and l-arginine phosphate) (see Fig. 4 and Supplementary Table S3).

Fig. 5

Plot of measured versus PLS-R-predicted oxygen levels, showing a high correlation between observed and predicted oxygen levels (partial least squares regression; y = 0.981x; R2 = 0.944, Q2 [cross-validated] = 0.605, P < 0.001 [1000 random permutations]). Plot is based on nymphs at their critical thermal maxima only (treatments 1–3).

Comparing the metabolism of nymphs under hypoxia and hyperoxia; controlling for temperature

To test whether the observed improvement of heat tolerance under hyperoxic conditions (Fig. 2) is associated with alleviation or postponement of anaerobic metabolism, we compared the metabolic profile of nymphs reaching their thermal maximum under hypoxic conditions (treatment 1) with that of nymphs at the same temperature but under hyperoxic conditions (treatment 4). The metabolic profiles were clearly different for both treatments as demonstrated by PLS-DA (Fig. 6). The optimal PLS model comprised three LVs, based on the minimization of cross-validated classification errors, and accounted for 36.7% of the total variance in the metabolic dataset (x-block). The cross-validated classification error was 14.4% with a significance of P < 0.001 when compared with 1000 random permutations. So there were consistent differences between hypoxia and hyperoxia in the perturbations in metabolic pathways in response to heat stress. Again, the metabolites differentiating between the two groups are associated with the generation of energy, involving both anaerobic metabolism (e.g., lactate and alanine) (see Fig. 4 and Supplementary Table S3) and amino-acid metabolism (e.g., tyrosine).

Fig. 6

PLS-DA scores plot, showing a clear separation in metabolic profile between nymphs under hypoxia (treatment 1) and hyperoxia (treatment 4) irrespective of the temperature (33.1°C at 14 kPa and 27.3°C at 5 kPa) (cross-validated class error rate of 14.3%, P < 0.001 [1000 random permutations]). The best model had three LVs in total, but for clarity only the first two axes are shown. The percentage of data variance explained by each latent variable is indicated. Ellipses are drawn by eye to improve visibility.

Discussion

The respiratory system of terrestrial arthropods consists of an open tracheal system, a highly branched system of tubes extending throughout the body and supplying oxygen directly to metabolically active tissues. The air-filled trachea opens to the external atmosphere and can be ventilated, which greatly facilitates gas exchange (Dejours 1981; Socha et al. 2008). Consequently, insects with open trachea can regulate oxygen intake and maintain aerobic scope while avoiding oxygen toxicity under widely varying conditions (Hetz and Bradley 2005; Harrison et al. 2006). These considerations could explain why evidence for oxygen limitation at thermal maxima has not been forthcoming in several tracheated arthropods (Klok et al. 2004; Stevens et al. 2010). In this study, we studied nymphs with closed tracheal systems and which rely on diffusion across their tegument for gas exchange. We extend previous research showing that results for terrestrial insects cannot be generalized to all insect taxa and life stages (Verberk and Bilton 2011).

In the stonefly studied here, hyperoxia improved tolerance to heat whereas hypoxia reduced heat tolerance (Fig. 1). More extreme levels of either hypoxia or hyperoxia did correspondingly change heat tolerance to a greater extent, although the effects were not symmetrical: Hypoxia had larger effects than hyperoxia (see also Verberk and Calosi 2012). In addition, heat tolerance differed among years by 1.1°C under normoxia, showing that upper thermal limits are somewhat plastic in these aquatic nymphs, perhaps to a greater degree than in terrestrial insects (Overgaard et al. 2011; Hoffmann et al. 2013). Nymphs were collected from the same field site, and field conditions did not differ dramatically between years (in fact conditions prior to collection were warmer and drier in 2011, contrasting with the lower upper thermal limits observed), suggesting that a greater plasticity in upper thermal limits might well be observed under more challenging conditions.

Importantly, we directly investigated the link between thermal limits and the expected onset of anaerobic metabolism when oxygen limitation sets in. Higher levels of lactate and succinate clearly indicate the onset of anaerobic metabolism at thermal maxima in stonefly nymphs under normoxic conditions (Fig. 4). Higher levels of alanine, produced together with lactate from pyruvate under conditions of limited oxygen, further support the onset of anaerobic metabolism. The pyruvate generated from glucose during glycolysis fits with the increased sugar-based metabolism we observed (e.g., pentitol phosphates). To reduce the accumulation of pyruvate, it can be transformed into acetylpyruvate using acetate (in addition to the transformation into lactate and alanine mentioned above). Although these specific end products were not known for D. cephalotes, they corroborate previous research on responses to hypoxia in Drosophila (Feala et al. 2007) and correspond to known by-products of anaerobic metabolism in other terrestrial insects (Hoback and Stanley 2001). Insects typically rely on the glycerol-3-phosphate shuttle to regenerate NAD+, which is important for maintaining glycolysis (Gilmour 1961), thus explaining the higher levels of glycerol and glycerol-3-phosphate. Our findings in Dinocras show good correspondence to those by Malmendaal et al. (2006) who also include alanine, acetate, and glucose as important markers for heat stress in Drosophila.

The higher dependency on glycolysis under oxygen-limited conditions is also evident from an overall disruption of the tricarboxylic acid (TCA) cycle with strong increases in oxaloacetate, malate, and fumarate, whereas levels of isocitrate decreased. Interestingly, the TCA cycle intermediate 2-oxoglutarate (alfa-ketoglutaric acid) was unaffected by the treatments, but was strongly related to the body size of individuals, with larger animals having higher levels (β = 0.843, t1,72 = 6.49, P < 0.001; Supplementary Fig. S3). Michaud et al. (2008) studied in heat stress in the Antarctic midge and suggested that increased 2-oxoglutarate levels were indicative of a perturbed TCA cycle. In addition, 2-oxoglutarate is involved in the formation of alanine from pyruvate in various metazoans, including sipunculid and polychaete worms (Müller et al. 2012). Larger animals also consistently showed reduced levels of acetylpyruvate (β = −0.365, t1,69 = 2.88, P = 0.0054). This suggests that metabolic pathways are size dependent, with the TCA cycle being perhaps more important relative to glycolysis in larger animals. Such size-dependent differences in metabolic pathways may underlie the size-dependent differences in both heat tolerance (Supplementary Table S2) and thermal sensitivity of metabolism found (Fig. 1; see also Verberk and Bilton 2011).

The metabolic shifts toward anaerobic metabolism described above in response to oxygen limitation were not sufficient to compensate for the decrease in energy production via aerobic metabolism. ATP levels indicate a lower energy status in nymphs at their thermal maxima (treatments 1–3). In addition, Sacktor and Cochran (1957) found that in addition to hydrolysis of ATP to ADP, the hydrolysis of UTP to release the stored chemical energy involved splitting all three phosphates, resulting in free uridine. Such generation of free uridine under conditions of ATP shortage could explain both the high levels of uridine under conditions at which energy status is compromised and the tight association between uridine and lactate accumulation. In addition, comparing nymphs at their thermal maxima (treatment 1–3), arginine phosphate levels declined in parallel with the shift toward anaerobic metabolism. Arginine phosphate is a well-established marker of cellular energy status, donating a phosphate to ADP when the ATP pool is depleted (Uda et al. 2006), explaining the positive correlation observed between ATP and arginine phosphate.

These results clearly show that under normoxia, the concept of oxygen limitation applies to the insects studied. The onset of anaerobic metabolism and the decrease in energy status, which eventually limit survival of heat stress, can all be traced back to shifts in the metabolome of the stonefly nymphs. Comparing only nymphs at their thermal limit (treatments 1–3), we found an attenuation of these shifts in metabolites under hyperoxia and a potentiation under hypoxia. As nymphs reached their critical thermal maxima at different temperatures under different oxygen conditions (Fig. 2), this indicates that perturbations of metabolic pathways in response to heat stress and hypoxia may be similar, but not identical (in which case changes in metabolism due to heat or hypoxia would be interchangeable and organisms at their thermal limit would have a very comparable metabolite profile). Possibly animals warmed under hypoxic conditions started to go into anaerobic metabolism quite some time before reaching their critical thermal maximum. Furthermore, when comparing the nymphs under hypoxia (treatment 1) and nymphs at the same temperature but under hyperoxic conditions (treatment 4), it is clear that hyperoxia rescued these animals and negated the severe energy shortage and anaerobic metabolism observed under hypoxia. Thus, the hypothesis that heat tolerance is limited by insufficient energy as a consequence of oxygen limitation is strongly supported at normoxic and especially hypoxic conditions. However, for nymphs at their thermal limit under hyperoxia, anaerobic metabolism and associated changes were not very distinct. Shortage of energy resulting from oxygen limitation could still be a factor involved in limiting heat tolerance, as there was a trend to lower ATP levels, higher alanine levels, and a general disruption of the TCA cycle (trend to higher fumarate and malate levels and lower isocitrate and oxaloacetate levels). Yet, the results allow an exploration to see whether other factors are involved in setting the level of heat tolerance.

Alternative mechanisms underpinning heat tolerance are related to membrane properties and protein denaturation (Feder and Hoffmann 1999; Pörtner 2001). The consequences of heat on membranes are increased fluidity and membrane leakage, which will promote proton leaks and lead to an uncoupling of mitochondria. This will lower ATP production. As no strong decline in ATP was observed under hyperoxia, even though this conformed to the highest temperature experienced by nymphs across all treatments, membrane fluidity is probably not the main driver. Nevertheless, we did find some evidence for disruption of membranes. Phosphocholine is an intermediate in the synthesis of phosphatidylcholine in tissues, a major component of biological membranes. Phosphocholine levels were elevated in all animals at their thermal limit, suggesting some generic membrane response to heat stress. Subsequent experiments would be necessary to confirm the importance of reorganization of membranes. Denaturation of proteins is another potentially important mechanism underlying heat tolerance (Somero 1995). In our study, we did see disruption of amino acid metabolism and nucleotide metabolism, but this was most likely due to active catabolism of proteins to generate energy, rather than protein denaturation. Thermal limits in terrestrial insects are usually higher, making protein denaturation more likely (see also Verberk and Bilton 2011; Verberk and Calosi 2012).

Specifically for what may cause heat coma under hyperoxia, we also have to consider the effects of radical oxygen species, which may be more prevalent during metabolic disruption. Levels of 2-hydroxybutyrate were elevated relative to control animals, especially under hypoxia and normoxia (treatments 1 and 2). Higher levels could indicate an upregulation of protection against the action radical oxygen species, as 2-hydroxybutyrate is released as a byproduct when cystathionine is cleaved to cysteine that is incorporated into the antioxidant glutathione. This fits with the strong disruption of the TCA cycle under hypoxia and normoxia. However, we found no compelling evidence for damage by radical oxygen species that set levels of heat tolerance levels under hyperoxic conditions.

Finally, there is the possibility that neuronal motor control is impacted in the nymphs before the onset of anaerobic metabolism and the associated organism-wide shortage of energy. Dawson-Scully et al. (2010) found that activation of the Protein kinase G (PKG) pathway increased survival in adult Drosophila, but decreased locomotory function. For aquatic damselfly nymphs, heat hardening was previously shown to increase survival but decrease locomotory control, possibly via activation of the PKG pathway (Verberk and Calosi 2012). Therefore, this pathway could uncouple the behavioral responses observed (Fig. 2) and the biochemical responses in the metabolome.

In conclusion, this study provides broad support for the hypothesis of oxygen limitation; anaerobic metabolites accumulated in nymphs at their thermal limit under hypoxic and normoxic conditions. The metabolomics approach showed that energy shortage and anaerobic metabolism could not be conclusively identified as the cause underlying thermal limits under hyperoxia. Metabolomics proved useful for suggesting a range of possible mechanisms to explore in future investigations, such as the involvement of leaking membranes or free radicals. In doing so, metabolomics provided a more complete picture of changes in metabolism under hypoxia and heat stress.