Plasticity of salmonfly (Pteronarcys californica) respiratory phenotypes in response to changes in temperature and oxygen

Climate change is altering the abiotic characteristics and increasing the temperature of mountain streams around the globe (Webb and Walling, 1992; Ashizawa and Cole, 1994; Hari et al., 2006; Moatar and Gailhard, 2006; Pekarova et al., 2008; Isaak et al., 2012). Altered conditions affect mountain-dwelling taxa directly and also, in some cases, indirectly via interactions with other biotic and abiotic factors (Birrell et al., 2020; Shah et al., 2020; Frakes et al., 2021). For example, rising temperatures can stimulate metabolic demand for oxygen more than they raise rates of supply, resulting in oxygen insufficiency (Verberk and Bilton, 2013).

Whether populations will persist in any particular location is difficult to predict. One possibility is extinction (Rosset and Oertli, 2011; Poff et al., 2012); however, if they persist, other potential responses include evolution (if it is rapid enough to keep up with the pace of climate change; Diffenbaugh et al., 2018; Wilczek et al., 2014), shifts in ranges including uphill movement (if suitable new habitats exist; Giersch et al., 2017; Parmesan, 2006; Shah et al., 2020), and physiological and behavioral plasticity (Chevin and Hoffmann, 2017). Adaptive plasticity is especially interesting, as it may occur rapidly within individuals and may allow populations to persist in their current locations, while buying time for longer-term ecological or evolutionary adjustments (Chevin and Hoffmann, 2017). For example, stonefly nymphs acclimated to warmer water typically have higher critical thermal maxima (CT_max) than those acclimated to cooler water (J.I.F., unpublished data; Shah et al., 2017), which may allow warm-acclimated populations to better survive transient exposure to high-temperature extremes. More broadly, although plasticity may alter the thermal limits of aquatic insect respiratory physiology, its quantitative importance is unknown.

Like most aquatic insects, stonefly nymphs (Order Plecoptera) support metabolism with oxygen extracted from water. Obtaining sufficient amounts can be difficult because oxygen is sparingly soluble, and water is dense and viscous enough to make rapid ventilation of respiratory surfaces potentially impractical and costly (Woods and Moran, 2020). Recent experiments support the idea that the morphology and physiology of respiratory surfaces influence upper temperature tolerance (Verberk and Bilton, 2013). The gills, which occur in many species, are thin-walled extensions of the cuticle containing many small tracheal tubes connected to larger trunks in the thorax and abdomen (Wichard and Komnick, 1974). Tracheal gills are key sites of gas exchange, but whether gills show plasticity is largely unknown for Plecoptera and other aquatic insects. By contrast, terrestrial insects held in experimental hypoxia are well known to develop larger tracheal trunks (Wigglesworth, 1954; Locke, 1958; Loudon, 1989; Henry and Harrison, 2004; VandenBrooks et al., 2012; Harrison, 2015) and greater tracheolar branching (Jarecki et al., 1999). Plasticity in gill morphology may allow nymphs to extract oxygen more readily from warmer or more hypoxic water and consequently to tolerate higher temperatures occurring with climate change.

Here, we investigated the degree to which plasticity shapes respiratory and metabolic phenotypes in an aquatic insect. Because stonefly nymphs take up oxygen across specialized respiratory surfaces – the gills – we expect that those surfaces will show substantial, and possibly adaptive, plasticity. To test this prediction, we examined the performance of stonefly nymphs during and after acclimation to multi-week exposures to different experimental combinations of temperature and oxygen. We measured rates of growth and survival, and then examined plasticity in several underlying respiratory phenotypes: (1) the upper CT_max of individuals, (2) the oxygen sensitivity of their metabolic rates and (3) the morphology of tracheal gills.

We measured responses of Pteronarcys californica Newport, 1848 nymphs to manipulations of temperature and oxygen in an experiment that included four treatments (2×2; Table 1). Briefly, insects were weighed and then held at 9°C or 16°C and subjected to normoxia (19.5 kPa, ∼100% saturation) or hypoxia (11.7 kPa, ∼60% saturation) for 6–8 weeks of acclimation. Subsequently, individuals were removed and weighed for growth. Half of the individuals from each treatment group were then used for metabolic rate measurements (n=45) and half for CT_max experiments and gill morphology analysis (n=44; Table 1).

All P. californica nymphs (0.04–1.24 g) were collected from lower Rock Creek, a relatively pristine stream that flows north out of the Sapphire and Pintler mountain ranges (46°41′52.908″N, 113°40′10.272″W; Fig. 1). Flows in Rock Creek range from ∼150 to 2000 cfs (∼4–57 m³ s⁻¹) and water temperatures range from ∼0 to 20°C (United States Geological Survey, gauge station 17010202). Like other small, turbulent streams, Rock Creek is well oxygenated, yet the risk of lower oxygen levels remains with climate change. Insect collection took place in December 2019 near Clinton, MT, USA, from riffle habitats using a screen barrier net (91×91 cm) and insects were transported to the University of Montana (1000 m above sea level) in aerated stream water. They were then held together in buckets and fed conditioned leaves (leaves that had been in-stream long enough to be soft and colonized with biofilms) from Rock Creek and placed within incubators (Percival Scientific, I-66LLC8, Perry, IA, USA) at 12°C overnight prior to the 6–8 week acclimation period. The light cycle was set to 12 h:12 h light:dark both within the incubators and during the acclimation period.

Tupperware chambers (280 ml) were used to hold nymphs in a 2×2 design, with high and low temperatures (16 and 9°C) crossed with normoxia and hypoxia (∼60% of air saturation). Nymphs in the high-temperature treatments were held for 50 days and those in the low-temperature treatments for 64 days prior to removal for measurements. Normoxia was imposed by bubbling small streams of air into chambers, and hypoxia (approximately 60% of air saturation) was imposed by bubbling streams of air mixed with nitrogen gas. Nitrogen flow was controlled by a mass-flow controller (UFC-1100, Unit Instruments) connected to controlling electronics (MFC-4, Sable Systems). The nitrogen stream was mixed, via a T-connector, with a stream of air driven by an aquarium pump, and the mixed stream was then distributed to the chambers via a custom-built manifold. Levels of oxygen in normoxic and hypoxic chambers were measured daily using a Firesting O₂ probe (PyroScience, Aachen, Germany). Values in hypoxic chambers were always within 5% of the target value (60% of air saturation).

Prior to placement in the chambers, 96 nymphs were weighed individually (alive but gently blotted dry). They were then separated into groups of three, including three individuals of distinguishable sizes (small, medium and large), which were placed together in a single chamber with three conditioned leaves and one pebble (32 total chambers). Half of the chambers were held at 16°C (16 chambers, 48 individuals) and the other half at 9°C (16 chambers, 48 individuals). Then, within each temperature treatment, half of the chambers were held in normoxia (19.5 kPa, eight chambers, 24 individuals) and half in hypoxia (11.7 kPa, eight chambers, 24 individuals).

During the acclimation period, conditioned leaves and fresh stream water from Rock Creek were placed in each chamber every 10 days. Before these routine changes, temperature and oxygen levels were verified using a thermocouple and the Firesting O₂ probe, respectively. At the end of the experiment, individuals were weighed again, then half from each treatment (n=12) were transported to the Flathead Lake Biological Station (FLBS) for metabolic rate measurements (see below). The other half (n=12 from each treatment) were transferred to a separate incubator at the University of Montana set to an intermediate temperature of 12°C and held without food for 3 days, before CT_max measurements were conducted to ensure nymphs were not actively digesting during the CT_max experiment (see below; Shah et al., 2017). Finally, these nymphs were fixed in glutaraldehyde for microscopy analysis of gills (see below).

Nymphs were monitored every 10 days for survival, and dead individuals were removed when found. Growth rate was estimated as final minus initial mass divided by the number of days in treatments (50 days in the warm treatment and 64 days in the cool).

All analyses were run in R (https://www.r-project.org/). We analyzed survival using logistic regression implemented as generalized linear models (glm, family ‘binomial’). We analyzed growth rates using linear mixed-effects models to account for the blocking of individuals into individual containers. Growth rate was analyzed as a function of initial mass, oxygen treatment (normoxic or hypoxic) and temperature (warm or cool), including all interactions.

Tests of CT_max were conducted in a Plexiglas water bath (55×32×14 cm) submerged in an insulated container that held ∼18 l of water (Fig. S1). At the start of each CT_max test, cool water (12°C) was added to the water bath and individuals were transferred into partially submerged, numbered mesh chambers (12 in each run), in which they were allowed to adjust to the new conditions for 2 min. Temperature was then ramped at +0.3°C min⁻¹ using a 500 W titanium heating rod controlled by a temperature controller (Dwyer 16B-22, Dwyer Instruments, Michigan City, IN, USA; Shah et al., 2017). Normoxia was maintained throughout the tests by bubbling in air with an air stone and 2.5 W pump (Imagitarium, Petco Brand, www.petco.com). Gradients in temperature or oxygen were minimized by circulating the water with a pump (Rio Plus 600, Rio Production, Taipei, Taiwan, www.rio-pump.com). Temperature was monitored and recorded with an Omega digital thermometer (HH 802W, Omega Engineering, Norwalk, CT, USA) and platinum probe. Once the water temperature reached ∼15°C, individuals were gently flipped onto their backs with forceps once per 1°C increase in temperature, and we noted the temperature at which they lost their righting response (Lutterschmidt and Hutchison, 1997). Individuals displaying loss of righting response were immediately removed and placed into cooler (12°C) water for recovery. Afterward, they were transferred to marked cups and preserved in glutaraldehyde for later analyses of gill morphology. CT_max data were analyzed using linear mixed-effects models in R, with body mass, acclimation temperature and acclimation oxygen level (normoxia or hypoxia) as predictors.

Metabolic rates were measured for 87 nymphs (i.e. approximately half of the individuals from each acclimation treatment) using intermittent respirometry (Table 1). Nymphs were transported to FLBS in the experimental Tupperware chambers and, on arrival, transferred to perforated Tupperware chambers, which were placed in a 12.5°C constant-temperature water bath (stream water cooled by a refrigerator, ActiveAqua Hydro-Culture, Active Aqua, www.activeaquahydroponics.com). Individuals from all treatments were held at this common temperature for 2.5–5.5 days prior to measurements. Holding times varied because only six individuals could be measured at once in the respirometry system. Thus, some individuals were held longer than others. To avoid confounding acclimation treatment with differences in holding time, we measured metabolic rates for individuals from each acclimation treatment during each trial.

The respirometry system consisted of two Plexiglas boxes housing eight horizontal glass mini-chambers and a water bath reservoir to control temperature (Loligo Systems, Aarhus, Denmark). Three individual stoneflies were placed into three of the mini-chambers in each box, leaving a blank mini-chamber control in one box and an open chamber for oxygen sensing and regulation in the second box (Fig. 2A). Small (2.2 ml) and large (17 ml) mini-chambers were used, depending on the size of nymph being tested. Filtered (0.2 µm) water, held constant at 12.5°C by a refrigerated unit, circulated continuously between the water bath and the Plexiglas boxes. Witrox 4 oxygen meters monitored dissolved oxygen concentration, which was measured in each mini-chamber with polymer optical bare-tip fibers pointed at 2 mm sensor spots (thin planar oxygen mini-sensors mounted onto the inside of each mini-chamber, PreSens Precision Sensing). Sensor spots were calibrated for a zero baseline using a sodium sulfite solution (10 g per 500 ml water) and 100% saturation baseline in bubbled recirculated water. During trials, oxygen concentration and temperature were recorded every second using AutoResp software (Loligo Systems). Initial oxygen saturation in normoxia was 11–12 mg l⁻¹ for all experiments. Ambient barometric pressure ranged from 90.2 to 92.4 kPa, or 18.9–19.4 kPa of oxygen, at the laboratory (elevation of 892 m).

The mini-chambers were connected to flush pumps controlled by the program AutoResp for intermittent measurements. Each intermittent cycle consisted of a flush period (90 s) during which water was pulled from the box into the mini-chambers, followed by a waiting period (30 s) that allowed oxygen levels in the chambers to stabilize prior to measurement, and finally a closed period (600 s) during which pumps were off and oxygen was depleted by stoneflies (from which metabolic rates were later calculated). Measurement windows of 600 s for each cycle provided data with high R² values, and allowed us to make replicate measurements (four for each oxygen level).

To test how sensitive metabolic rates were to hypoxia, we measured metabolic rates at progressively decreasing oxygen levels (in the same order for every stonefly). We had four replicate measurement periods (from which to calculate metabolic rates) at each oxygen level of normoxia (∼10 mg l⁻¹), and progressively increasing levels of hypoxia (6, 4 and ∼2 mg l⁻¹). Oxygen levels were lowered by bubbling nitrogen gas into the water bath reservoir and raised by bubbling air. At the end of each trial, we removed the individuals and documented survival. Each trial took ∼5 h, and we completed one trial per day.

During metabolic trials, stonefly behavior was noted once per minute on a scale from 0 to 5 indicating the degree of movement (0, upside down and still; 1, right-side up and still; 2, moving legs or antennae; 3, slow crawl; 4, fast crawl; and 5, rapid or frantic movement). We also documented push-ups (i.e. the body is raised and lowered or moved side to side repeatedly) if they occurred, because they are a strong indicator of oxygen insufficiency in P. californica and other stoneflies (Hynes, 1976; Genkai-Kato et al., 2000).

Metabolic rates were estimated from raw data on changes in oxygen concentration following methods used by Malison et al. (2020a). Oxygen consumption was calculated as the slope of the linear regression of oxygen concentration with time, and data were corrected by subtracting the slope of the control chamber from each experimental chamber. Mean mass-specific metabolic rates (MSMRs) were calculated for each individual for each oxygen level from the replicate 10 min windows. Negative metabolic rates resulted when control chambers had larger slopes than measurement chambers, which we consider to be experimental error. This occurred most commonly when respiration rates of stoneflies were very low.

We used linear mixed-effects models to analyze how acclimation treatments influenced MSMRs. Because individual metabolic rates were measured repeatedly, we included individual as a random effect. Because we predicted similar metabolic rates at higher experimental dissolved oxygen (DO) levels, we first tested the effect of acclimation temperature and oxygen level on metabolic rates for the three highest oxygen levels. A second model tested whether MSMRs differed between the three highest experimental oxygen levels (all above 4 mg l⁻¹) compared with the acute hypoxia level (nominal 2 mg l⁻¹). A simple linear model was then used to further analyze treatment effects within the nominal 2 mg l⁻¹ level of experimental DO. Next, we investigated whether behavior changed with experimental oxygen level using a linear mixed-effects model with acclimation temperature, acclimation oxygen and experimental oxygen as the main effects. We then included behavior as a main effect along with acclimation temperature, acclimation oxygen and experimental oxygen to investigate how all variables influenced MSMR. We primarily focused analyses on MSMRs to avoid overparameterization of models; however, because mass is an important variable that may influence metabolic rates, we ran additional models for absolute metabolic rates that included mass, acclimation oxygen, acclimation temperature and experimental oxygen level as main effects.

The morphology of individual gill filaments was quantified using a combination of light microscopy and scanning electron microscopy (SEM). Three live P. californica individuals (small, medium and large body sizes) were collected from each treatment group after CT_max measurements. Prior to imaging, individuals were rinsed with clean water and their gill tufts were gently brushed with a fine paintbrush to separate clumps of filaments and remove sediment.

To preserve the integrity of gill tufts, individuals were immobilized ventral-side up using rubber bands threaded through holes drilled in a plastic Petri dish filled with water. Legs were kept out of the image with forceps or mounting pins. Gill tufts were imaged using a stereomicroscope (Nikon SMZ 1500, Nikon Metrology, Brighton, MI, USA) illuminated with fiber-optic light. Images used for measurements were taken at 60× magnification.

Stoneflies were processed for SEM using a typical glutaraldehyde fixation technique. Briefly, individuals were placed into a 3.75% glutaraldehyde solution for a minimum of 24 h at room temperature, then water was removed using an ethanol dehydration series (30%, 50%, 70%, 90% and 100% ethanol). Samples were incubated in anhydrous ethanol overnight and dried using a critical point dryer (CPD-030, Oerlikon Balzers, Schaumburg, IL, USA). The thorax and gills were then mounted on an aluminium stub and sputter coated (Desk V Standard, Denton) prior to imaging. Gills were imaged at a low accelerating voltage (3–5 kV) using a S-4700 Cold-Field Emission SEM (Hitachi, Schaumburg, IL, USA).

On each stonefly, we measured the length and width of a total of 60–150 individual gill filaments. Gill filaments were selected from six gill tufts per individual, three on each side from the upper, middle and lower thorax (Fig. 2B,C). Single gill tufts generally contained between 60 and 120 individual gill filaments. However, many individual filaments within a tuft were difficult to distinguish. To maintain accuracy, only sufficiently visible and in-focus tufts were measured and analyzed, resulting in 60–150 gill filaments per stonefly. Filaments that were highly curved or that were anomalously short or long (e.g. short and stubby gills in the early growth stage, or filaments that clearly were damaged) were excluded. The lengths and widths of gill filaments were extracted from images using the Fiji image-processing application (Schindelin et al., 2012). Gill filament length was measured from the tip to its attachment point at the tuft. Filament thickness was measured at the midpoint between tip and base.

We then calculated mean filament width and filament length for each individual, as well as mean filament area (by treating each filament as a cylinder and calculating its surface area). Finally, on 10 individuals, viewed via a stereomicroscope, we counted all filaments on each of the thoracic segments. Individuals contain additional, uncounted filaments on some abdominal segments, but our approach captured much of the area and allowed analysis of size scaling and effects of experimental treatments.

Gill filament width, length and number were analyzed separately using linear models, with body mass and experimental exposure to oxygen and temperature as predictors (including their interactions). Because not all parameters were available for each individual, we were unable to calculate total filament area per individual, Nevertheless, we were able to combine the scaling coefficients for each parameter to estimate the scaling coefficient for total gill area on thoracic segments using the equation for a cylinder: total area≈number of gill filaments×mean filament length×width².

Survival was high for the duration of the experiment in all acclimation treatments; only four individuals died in the warm temperature treatment (Fig. 3). Preliminary analyses indicated that body mass did not predict survival, so it was dropped from further analyses. A logistic regression with acclimation temperature and acclimation oxygen as predictors indicated that temperature was significant (χ²=5.72, P=0.017) but that the main effect of oxygen and the interaction term were not (χ²=0, P=1).

Growth rates depended in complex ways on body mass, acclimation temperature and acclimation oxygen level (Table 2, Fig. 4; Fig. S2, Table S1). Overall, the smallest individuals consistently gained mass, whereas the largest individuals primarily lost mass. This pattern was modified by oxygen availability both as a main effect (those in hypoxia grew less) and in interaction with mass (large individuals in hypoxia lost disproportionately large amounts of mass) (Table 2). Temperature did not have a main effect, but it did interact significantly with both mass and oxygen. These interactions reflect that growth rates of individuals in the warm treatments were disproportionately depressed by large body size and hypoxia.

The mean value of CT_max across treatments was 31.4°C, and most values lay between 30 and 33°C (Fig. 5). Linear modeling of CT_max as a function of initial stonefly mass and acclimation oxygen and temperature levels indicated that mass was not a significant predictor, so it was dropped from further analyses. A simplified model with acclimation oxygen and acclimation temperature suggested that CT_max was significantly higher in individuals acclimated in hypoxia (F=5.99, P=0.019). Although those from the warm treatment had higher CT_max, the main effect of acclimation temperature was not significant (F=1.86, P=0.18) nor was the interaction term between acclimation oxygen and acclimation temperature (F=0.01, P=0.91).

Pteronarcys californica from all acclimation treatments had similar MSMRs at the three highest levels of experimental DO (4, 6 and 10.5 mg l⁻¹; P>0.12; Figs 6 and 7). Because MSMRs were similar at higher experimental DO levels, we compared MSMRs from the three higher experimental DO levels with MSMRs measured in acute hypoxia. The main effect of experimental DO level (2 mg l⁻¹ versus high DO) was significant (F=10.91, P=0.001), indicating that individuals had lower metabolic rates at 2 mg l⁻¹ than at higher DO levels. The interaction of acclimation temperature and experimental DO level was significant (F=7.67, P=0.007), illustrating that the high-temperature acclimation group had higher MSMRs at 2 mg l⁻¹ oxygen and lower MSMRs in normoxia. Within the 2 mg l⁻¹ experimental oxygen level only, the effect of acclimation temperature was significant (F=6.12, P=0.018). Mass was not included in the main analysis to avoid overparameterization of the models. In additional models for absolute metabolic rates (that included mass, acclimation oxygen, acclimation temperature and experimental oxygen level as main effects), mass had a highly significant effect (F=50.29, P<0.0001) and acclimation oxygen level was also significant (F=4.424, P=0.0434). There were many significant interactions, including mass×acclimation oxygen level, acclimation oxygen×experimental oxygen level, mass×acclimation oxygen×acclimation temperature, mass×acclimation oxygen×experimental oxygen, and mass×acclimation oxygen×acclimation temperature×experimental oxygen level (F>5.686, P<0.0198; Fig. S5).

Stonefly behavioral scores were similar across experimental DO levels and for oxygen acclimation levels (F=0.96, P=0.331, and F=0.08, P=0.78, respectively; Fig. S3), but differed significantly by acclimation temperature (F=4.87, P=0.034), reflecting that individuals that acclimated to the warm temperature were more active during measurements. Our final model tested the main effects of behavior, acclimation temperature, acclimation oxygen and experimental oxygen on MSMR. Behavior significantly influenced MSMRs (F=15.48, P=0.0001; Fig. S4) and the interaction of experimental oxygen level and acclimation temperature was significant (F=11.20, P=0.001).

The morphology of gill filaments was affected in complex ways by body size and treatments. Body size was the most important predictor of filament morphology (Fig. 8; Table 2); larger individuals had more and larger filaments. For filament width and number, acclimation to experimental factors (temperature or oxygen) resulted in no significant shifts. By contrast, filament length was significantly affected by both acclimation oxygen and temperature. Filaments were shorter in individuals acclimated to the warm treatment and in individuals acclimated to the hypoxic treatment. Thus, individuals in cold normoxia had, for their body sizes, the longest gill filaments (Fig. 8).

Because of missing data, we were unable to calculate total thoracic gill area for each individual. Nevertheless, by combining scaling coefficients for length (0.152), width (0.168) and number (0.393) (using the equation: total gill surface area≈2π×radius×length×number of filaments), we estimate that total thoracic gill area scaled with a coefficient of 0.713.

How will aquatic insects respond to climate change? This question has generated a robust literature over the past several decades (Rank and Dahlhoff, 2002; Parmesan, 2006; Rosset and Oertli, 2011; Poff et al., 2012; Wilczek et al., 2014; Giersch et al., 2017; Diffenbaugh et al., 2018; Birrell et al., 2020), which has examined possible responses, including extinction, range shifts, rapid evolution and plasticity. Although plastic responses may play an important role, we have little empirical data on which traits change or the magnitude of those changes (Sgrò et al., 2016). In this study, we explored the potential roles of plasticity in respiratory and metabolic phenotypes in an aquatic insect. Because P. californica nymphs take up oxygen across the gills, we expected that gills would show plasticity. We subjected juvenile stoneflies (nymphs) to ecologically relevant combinations of temperature (9 or 16°C) and oxygen (air saturated: 19.5 kPa at our elevation in the Rocky Mountains; or hypoxic: 11.7 kPa or 60% of air saturation), which are commonly experienced in nature, and then measured performance (survival and growth) and a set of phenotypes linked to respiratory physiology. Despite relatively small differences in temperature and modest changes in oxygen levels, nymphs showed clear differences in growth among treatments. Larger individuals grew less (or lost mass) compared with smaller individuals, and these effects were magnified for individuals in the warm water treatment and for individuals in hypoxia. Lower growth rates from the largest individuals may have occurred because final instar nymphs go through a period of dormancy, with little or no feeding and growth, prior to emerging as adults (Townsend and Pritchard, 1998). However, survival overall was quite high, with only four mortalities occurring in the warm treatment.

An ongoing discussion centers on two questions about upper critical thermal limits (CT_max): how much scope do aquatic insects have for plasticity in CT_max, and are upper limits set, or at least influenced by, oxygen insufficiency? Evidence for the first question generally shows that aquatic insects held at higher (but not stressful) temperatures exhibit higher CT_max, by about 1–4°C (Heiman and Knight, 1972; Calosi et al., 2008; Dallas and Rivers-Moore, 2012; Shah et al., 2017; Verberk et al., 2018; but, for a counterexample, see Treanor et al., 2013). In general, species adapted to higher temperatures may have reduced ability to acclimate (the tolerance–plasticity tradeoff hypothesis; Stillman, 2003; van Heerwaarden and Kellermann, 2020), although the opposite pattern has been found among species of Deronectes beetles (Verberk et al., 2018). More broadly, temperature-driven plasticity in upper thermal limits is observed widely across other taxa besides insects, and also in terrestrial and marine habitats (Gunderson and Stillman, 2015). Gunderson and Stillman (2015) estimated that, compared with terrestrial taxa, aquatic taxa have approximately twice the capacity to acclimate to changes in temperature. They suggested that lower levels of local thermal variation available to aquatic taxa force them to rely more strongly on physiological acclimation and less on behavioral thermoregulation. This is likely also true for P. californica in rivers in western Montana, in which temperatures often vary by only 1–2°C across 10 m stretches (Jackson H. Birrell, The University of Montana, unpublished data). Thus, in nature, individual P. californica have restricted opportunities for behavioral thermoregulation.

The second question – whether oxygen plays roles in setting upper limits – has also received increasing attention (Verberk et al., 2016). Theoretically, oxygen could play a role if, with rising temperatures, metabolic demand for oxygen outstrips the rate of oxygen supply from the environment (Woods, 1999; Verberk et al., 2011; Pörtner et al., 2017; Birrell et al., 2020). Most tests of this idea have measured whether CT_max changes with experimentally altered levels of oxygen. Typically, raising oxygen levels (hyperoxia) has little effect, whereas lowering oxygen levels (hypoxia) strongly depresses CT_max. For example, upper thermal limits to performance of diving beetles and damselflies fell by 4–8°C at 5 kPa oxygen compared with 21 kPa (Verberk and Calosi, 2012; Verberk et al., 2018). In P. californica, Frakes et al. (2021) found that CT_max fell continuously with declining levels of oxygen, reaching values that were more than 10°C lower at oxygen levels below 25% of air saturation (∼5 kPa). This same study also showed that CT_max was significantly lower in still water than in moderately flowing water, which the authors interpreted as indicating that low- or no-flow conditions impeded delivery of oxygen to the insects and contributed to oxygen insufficiency.

Few studies have examined whether aquatic insects ‘acclimate’ to altered levels of oxygen (as opposed to examining the acute effects of altered oxygen). In this study, we held nymphs for weeks at factorial combinations of temperature and oxygen, allowing us to ask whether acclimation to oxygen is possible. Interestingly, nymphs acclimated to the hypoxia treatments showed significantly higher CT_max, by about 1°C, compared with those acclimated to normoxia. Contrary to a number of other studies, individuals held at the warmer temperature did not show significantly higher CT_max, although they trended in that direction (J.I.F., unpublished data; Shah et al., 2017). Hypoxia-driven increases in CT_max may reflect acclimation of the respiratory system in hypoxia to promote better capture and transport of oxygen (see ‘Gill Morphology’, below).

Rates of metabolism integrate a broad range of influences, including body size, current biotic circumstances such as feeding history and infection status, and position in the life-history trajectory, as well as abiotic factors including temperature, water flow and oxygen availability. Of these factors, temperature and body size may have the most explanatory power (Gillooly et al., 2001).

Metabolic thermal performance curves can be used to determine how stressful different temperatures are to aquatic insects (Shah et al., 2021). Compared with the metabolic rates of some other benthic taxa, metabolic rates of P. californica nymphs in western Montana are less variable in response to warming temperatures, on average only ∼3.5 times higher at peak metabolic rates compared with rates at the coldest temperature (versus 3.5–7 times higher for other taxa; R.L.M., Brian K. Hand, Jack A. Stanford and Gordon Luikart, unpublished data). Metabolic rates of stoneflies decline with declining oxygen concentrations and, for P. californica, decline sharply at concentrations below ∼3.6 mg l⁻¹ (this study and Malison et al., 2020a,b). The main effects of acclimation temperature and acclimation oxygen did not significantly influence metabolic rates, but the interaction did. Individuals acclimated to the warm treatment had lower metabolic rates in normoxia and elevated metabolic rates in hypoxia. This pattern suggests that warm-acclimated individuals had increased ability to obtain oxygen in acute hypoxia, especially compared with P. californica nymphs from our previous studies (Malison et al. 2020a,b), which could not effectively extract oxygen at concentrations lower than ∼3.6 mg l⁻¹. It is also possible, however, that the increased MSMR in hypoxia indicates higher stress for individuals acclimated to warmer temperatures. It is not clear why individuals acclimated to the warm treatment had lower metabolic rates in normoxia, but it may be an acute effect of colder temperature as the metabolic rates were measured at 12°C (below the acclimation temperature of 16°C).

Insects distribute oxygen via the tracheal system, a ramified set of air-filled tubes of which the finest branches supply essentially all tissues with oxygen. Terrestrial insects typically regulate oxygen and carbon dioxide levels in their tracheal systems using short-term control over spiracles – valved openings on thoracic and abdominal segments (Förster and Hetz, 2010). Aquatic insects, by contrast, typically have closed (apneustic) spiracles, meaning that they cannot regulate oxygen via this mechanism. Adaptive plasticity over longer durations has been documented in terrestrial insects, with greater tracheal volume and increased tracheal branching occurring when individuals develop in hypoxia (Loudon, 1989; Jarecki et al., 1999; Henry and Harrison, 2004; Harrison et al., 2006; VandenBrooks et al., 2012). Aquatic insects may instead need to rely more on behavioral regulation in moderate hypoxia – by moving among flow gradients, actively ventilating their surfaces or gills by flapping the gills, movements of various kinds (including ‘push-ups’ exhibited by stoneflies), or movements of the whole body (Knight and Gaufin, 1963; Bäumer et al., 2000; Genkai-Kato et al., 2000; Van Der Geest, 2007). In deep hypoxia, they may use more complex physiological mechanisms (Harrison et al., 2018; Hoback and Stanley, 2001; Malison et al., 2020a,b), including changes in the number or morphology of tracheal gill filaments (Wichard and Komnick, 1974).

Our data, however, do not support this idea – there were no significant effects of acclimation treatment on filament width or number (although sample sizes were small, which decreases our power to detect effects). In addition, filament length was significantly shorter in the warmer acclimation temperature and in hypoxia. Shorter filament lengths will result in shorter diffusion distances between the site of oxygen uptake and the rest of the tracheal system. However, shorter filaments in conjunction with lack of systematic changes in width or number implies that total gill surface area was lower in the conditions (warmer temperatures, hypoxia) when metabolic demand for oxygen was higher and environmental supply lower. We suggest potential adaptive plasticity may have been constrained by the length of our treatments, which, although long by many experimental standards, still did not permit individuals to molt. Even longer experimental exposures may clarify the patterns and functional consequences of gill plasticity, if any.

Finally, we note that total thoracic gill area showed strong size scaling, with an estimated scaling coefficient of 0.713. Thus, total gill area scales in approximately the same way that metabolic rate scales with body size.

Taken together, these results suggest that temperature and oxygen interact to affect performance of stonefly nymphs in combinations of conditions that we expect to occur more frequently as climate change progresses. In small, turbulent streams, oxygen is not generally limiting. In larger, less turbulent systems, by contrast, oxygen levels can be depressed during periods of high temperature and low flow. In addition, photosynthesis and respiration strongly influence concentrations of oxygen in slow-moving rivers with abundant macrophytes (Allan and Castillo, 2007), and drawdowns can be exacerbated by anthropogenic loading of nitrogen and phosphorus (Chambers et al., 2006). Our study suggests that plasticity in traits at several levels of organization may partially ameliorate the severity of climate-driven risks.

Understanding how natural populations will respond to climate change is a pressing problem (Seebacher et al., 2015). Physiological plasticity may be a key mechanism of population resilience, but assessing its potential is difficult because it can occur across many different traits simultaneously at different levels of biological organization. In addition, plasticity may play out over time scales longer than those used in typical experiments (including ours), thus limiting inferences about responses in natural systems. Our study illustrates the power of a multi-trait approach over a relatively long experimental time frame. Our results suggest that acclimation may allow stoneflies to persist in warmer streams longer than expected.

Additional studies on more species over time scales long enough to allow insects to molt one or more times will provide even more ecologically relevant information on the degree to which acclimation will allow persistence in the face of climate change. We also suggest that future studies employ multi-trait approaches. Although such studies require more time and resources, they will likely reveal responses that are more relevant to the long-term effects of climate change in the wild.

We thank Jackson Birrell for help in collecting the samples, Jim Driver for support and the use of the University of Montana EMtrix electron microscopy facility, and October Moynahan for incubator support at the University of Montana Environmental Control for Organismal Research facility.