For insects, life in water is challenging because oxygen supply is typically low compared with in air. Oxygen limitation may occur when oxygen levels or water flows are low or when warm temperatures stimulate metabolic demand for oxygen. A potential mechanism for mitigating oxygen shortages is behavior – moving to cooler, more oxygenated or faster flowing microhabitats. Whether stream insects can make meaningful choices, however, depends on: (i) how temperature, oxygen and flow vary at microspatial scales and (ii) the ability of insects to sense and exploit that variation. To assess the extent of microspatial variation in conditions, we measured temperature, oxygen saturation and flow velocity within riffles of two streams in Montana, USA. In the lab, we then examined preferences of nymphs of the stonefly Pteronarcys californica to experimental gradients based on field-measured values. Temperature and oxygen level varied only slightly within stream riffles. By contrast, flow velocity was highly heterogeneous, often varying by more than 125 cm s−1 within riffles and 44 cm s−1 around individual cobbles. Exploiting micro-variation in flow may thus be the most reliable option for altering rates of oxygen transport. In support of this prediction, P. californica showed little ability to exploit gradients in temperature and oxygen but readily exploited micro-variation in flow – consistently choosing higher flows when conditions were warm or hypoxic. These behaviors may help stream insects mitigate low-oxygen stress from climate change and other anthropogenic disturbances.

Life in water is shaped by the scarcity of oxygen (Hutchinson, 1981; Lancaster and Downes, 2013). Water contains ∼30 times less oxygen than air (Jones, 1972). Water is also dense and viscous, which slows rates of diffusion and makes respiratory ventilation energetically expensive (Denny, 1993; Verberk and Atkinson, 2013). For water-breathing ectotherms such as aquatic insects, the challenge is to extract sufficient oxygen to fuel metabolism, despite low levels of environmental oxygen supply (Verberk et al., 2016a,b; Woods and Moran, 2020).

This oxygen problem is complicated further by variation in temperature and water velocity (Woods, 1999; Jacobsen et al., 2003; Verberk et al., 2011; Frakes et al., 2021; Verberk et al., 2016a). Warmer water modestly increases rates of oxygen supply – oxygen diffuses more rapidly in warmer water, which more than offsets declining concentrations – but raises organismal demand even more (Woods, 1999; Woods and Moran, 2008; Verberk et al., 2011; Verberk and Atkinson, 2013). At high temperatures, oxygen demand can exceed the capacity of organisms to obtain sufficient oxygen, resulting in declines in organismal performance and potentially death (Pörtner and Knust, 2007; Verberk et al., 2016a,b; Harrison et al., 2018). Higher flows, by contrast, increase rates of oxygen supply by thinning the boundary layers enveloping insect respiratory surfaces (e.g. gills or thin segments of the cuticle) (Hynes, 1970). Water velocity is zero (i.e. no-slip condition) at the cuticle–water interface but increases with distance from the body surface until it approaches the velocity of the free-stream flow (Vogel, 1994). When boundary layers are thin, as in fast-flowing and turbulent water, oxygen diffuses over shorter distances between the water column and respiratory tissues, increasing supply rates (Pinder and Feder, 1990; Hall and Ulseth, 2020). Thus, aquatic ectotherms can tolerate deeper hypoxia and warmer temperatures when water velocities are sufficiently high (Frakes et al., 2021).

Taken together, climate change and other anthropogenic disturbances, such as eutrophication and dams, are raising water temperatures, depressing oxygen concentrations and lowering water velocities (at least seasonally) in many rivers and streams. These changes threaten aquatic communities in part by depressing the ratio of oxygen supply to demand (Jacobsen, 2008; Deutsch, et al., 2015; Verberk et al., 2016a). Predicting how aquatic communities will respond, however, will require a better understanding of the mechanisms that aquatic insects use to mitigate oxygen limitation. A broad literature base now addresses the mechanisms by which ectotherm thermal performances, thermal limits and acclimation responses depend on oxygen and flow (Verberk et al., 2011, 2016a,b; Rubalcaba et al., 2020; Collins et al., 2021; Frakes et al., 2021). Much less of the literature has focused on behavioral mitigation of the oxygen problem (but see Verberk and Bilton, 2015), even though behavior is well known to provide small terrestrial and intertidal ectotherms with the means to thermoregulate (Ebersole et al., 2001; Helmuth et al., 2005; Scheffers et al., 2014; Woods et al., 2015; Birrell et al., 2020; Pincebourde and Woods, 2020). Whether aquatic insects can mitigate oxygen shortages by moving among microclimates will depend on (i) whether diverse conditions are available locally and (ii) whether insects actually exploit that diversity. We expect that these two conditions are related to one another – in the sense that aquatic insects should evolve behavioral mechanisms to exploit microclimatic mosaics only when lineages are exposed regularly to biologically meaningful gradients in nature. Assessing both the availability of microclimates and the capacity for individuals to exploit them will be critical to understanding how aquatic insects will respond to future climate change, and has so far received little attention (e.g. Helmuth et al., 2005; Pincebourde and Woods, 2020).

Whether temperature, oxygen and flow are likely to show strong local gradients depends on complex interactions among physical, hydrogeomorphic and chemical factors. Because water has a high heat capacity, equilibrates slowly with local gas partial pressures, and is typically well mixed (in lotic systems), we predict that microclimatic diversity in temperature and oxygen will be slight (Denny, 1993). Warm microclimates can form near warm springs or near shallow, slow-moving stream margins, where solar radiation is high, and cool microclimates can form in association with groundwater inputs, side channels, springs, deep-water pools and shading (Mosley, 1983; Matthews and Berg, 1997; Clark et al., 1999; Ebersole et al., 2001, 2003). However, the spatial grain at which these differences are manifested appears too coarse for individual insects to exploit (e.g. Ebersole et al., 2001). Moreover, cold patches are often associated with inputs of hypoxic groundwater, which would at least partially counteract the benefits of colder water (Matthews and Berg, 1997; Torgersen et al., 1999; Elliott, 2000). Small gradients in temperature and oxygen may also form in stream riffles in summer as a consequence of patterns of exchange between in-stream and hyporheic flows, with slightly cooler, more oxygenated water at riffle heads and slightly colder, less oxygenated water at the tail (from upwelling of hyporheic water) (Davy-Bowker et al., 2006).

By contrast, microclimatic diversity in flow velocities is often very high in lotic systems. In streams, flow velocities often vary from ∼0 cm s−1 within the substrate to >100 cm s−1 in the free stream environment (White, 1990; Comiti et al., 2007), and individuals may often be able to alter local flows past their bodies by moving just a few centimeters. Indeed, gradients in flow velocity strongly drive stream insect microdistributions (Mérigoux and Dolédec, 2004). Some aquatic insects choose more exposed surfaces and perform more respiratory movements when flows are low or temperatures high, presumably to increase rates of oxygen supply (Kovalak, 1976, 1979; Wiley and Kohler, 1980; Genkai-Kato et al., 2005). Higher flows can also increase the upper temperature tolerances and lower oxygen limits of stream insects (Frakes et al., 2021).

Here, we examined potential behaviors by which aquatic insects can mitigate oxygen shortages. We first measured the spatial diversity of temperature, oxygen saturation and flow velocity available to aquatic insects within riffles at small spatial scales in two streams in western Montana, USA. We then established experimental gradients of each factor (temperature, oxygen, flow) in laboratory flow chambers and measured how strongly nymphs of the giant salmonfly, the stonefly Pteronarcys californica, chose among them. Broadly, we expected that when nymphs are subjected to oxygen deficiencies, they should move to ameliorating microclimates (lower temperatures, higher concentrations of oxygen, or higher flows). However, we also predicted that nymphs would make stronger choices among levels of the factors, e.g. flow, that show strong spatial gradients in nature.

Field microclimate measurements

In July 2019, we sampled water temperatures, oxygen saturation and flow velocities on the Blackfoot River and Rock Creek in western Montana, USA (see Fig. 1A for schematic of sampling methodology). The Blackfoot River and Rock Creek have similar hydrology and geomorphology. Both are large, cold-water, cobbled mountain streams fed by snowmelt, precipitation and groundwater. The Blackfoot River originates in the Flathead Mountains and flows ∼120 km southwest through open sage-prairie, rangeland and canyons into the Clark Fork River near Missoula, Montana. From its west fork, Rock Creek flows ∼100 km north through the Sapphire Mountains, cutting through rangeland and narrow canyons. It joins the Clark Fork River 40 km upstream of the Blackfoot.

We sampled at 9 sites on the Blackfoot River and 10 on Rock Creek. Sites were separated at relatively even intervals (∼10 km) across the elevation gradients of both rivers (1000–1630 m and 1090–2060 m, respectively). At each site, microclimates were measured surrounding 15 large cobbles (15–50 cm wide) within a 10×5 m quadrant in the middle of riffles. Cobbles were chosen haphazardly every meter along three transects within the quadrant. Transects were separated by 5 m (i.e. one each at 0, 5 and 10 m within the quadrant). We measured water temperature and oxygen saturation under the bottom surface toward the rear and under the bottom surface toward the front of each cobble and in the free stream ∼3/4 of the way up from the stream bed to the water surface. We measured flow velocity on the top surface, at the leading edge and at the tail edge of each cobble as well as in the free stream ∼3/4 of the way up from the stream bed to the water surface. Measurements therefore represent different microclimates around and among cobbles. To measure temperature and oxygen saturation, we used temperature and fiber-optic oxygen probes (TSUB21 and OXROB3, Pyroscience, Aachen, Germany, respectively) connected to an optical temperature and oxygen meter (FSO2-C2, Pyroscience). The sensors were fitted into a hollow metal rod, which allowed us to insert the tips into the substrate. At each location, the tips were held in place for 30 s. Flows were measured using a portable velocity meter (FP111 Flow Probe, Global Water, Davison, MI, USA).

Laboratory choice experiments

In the laboratory, we measured the preferences of Pteronarcys californica Newport 1848 (salmonflies) nymphs for different temperature, oxygen and flow conditions (see Fig. 1B,C for schematic of experimental setup and Fig. S1 for photos of a salmonfly nymph and the flow chambers). Salmonflies were collected from the Blackfoot River or Rock Creek with a kick screen (91×91 cm with 1 mm mesh openings) and returned to the University of Montana, where they were transferred to 20-liter buckets filled with de-chlorinated tap water and held in a temperature-controlled incubator (I-66LLC8, Percival Scientific, Perry, IA, USA). Nymphs were held at ∼10°C and fed cottonwood (Populus spp.) leaves from the Blackfoot River or Rock Creek. Leaves were soaked in stream water (i.e. conditioned) for 2 weeks before being given to the nymphs, which softens them and helps bacterial and fungal biofilms establish on their surfaces, from which salmonflies derive much of their nutrition (Eggert and Wallace, 2007). Water in the buckets was oxygenated and stirred by air directed through air stones. Nymphs were held for a minimum of 2 days and up to 3 weeks before experimentation. During this pre-experimental period, there was no mortality and we observed vigorous feeding and frequent molting.

Oxygen choice experiments

We performed oxygen choice experiments during the spring and summer of 2019 by measuring where nymphs distributed themselves between two static levels of oxygen saturation (hypoxic versus normoxic) within an experimental aquarium (53×39×16.5 cm). The oxygen gradient was sustained by bubbling air or nitrogen gas through air stones at each end of the aquarium. We placed stainless steel coils connected to cooling recirculating water baths to control the temperature (1160A, VWR Scientific, Radnor, PA, USA, and AACH10HP, Active Aqua, Petaluma, CA, USA) and small water pumps to provide flow (PL-118, 200 l h­−1, PULACO, Guangzhou, China). Mesh walls divided the air stones, coils and pumps from the compartment in which the nymphs were kept, while still allowing temperature- and oxygen-controlled water to circulate. Plexiglas walls divided the central chamber into five lanes and kept the nymphs separate. Median Plexiglas walls divided the aquarium (and its five lanes) in half, preventing the normoxic and hypoxic water from mixing. However, nymphs were able to move between the hypoxic and normoxic sides by crawling through a 1-cm diameter hole cut in the bottom-center of each median wall. Cobbles and conditioned cottonwood leaves were placed on both sides of all lanes to provide food and substrate for the nymphs. High-grit sandpaper (P150, Gator, Fairborn, OH, USA) glued to the floor of the chamber provided footing.

During oxygen experiments, we measured choices made by 25 nymphs – 13 at 5°C and 12 at 16°C. At the onset of each trial, one nymph was placed on a random side (left or right) of each lane and left to acclimate overnight, during which oxygen levels were held at ∼100% of air saturation (∼19.0 kPa oxygen at our elevation, 1000 m.a.s.l.) on both sides of the aquarium. On the following morning, the oxygen level on a randomly chosen side of the chamber was decreased to ∼30%, 50% or 70% of saturation (∼5.7, 9.5 and 13.3 kPa, respectively) by bubbling in nitrogen gas, or kept at ∼100% of saturation (∼19 kPa) (control). Nymphs were left in the chamber for 2–3 h, after which we recorded their locations (hypoxic versus normoxic side) and reset the oxygen level on each side of the aquarium to 100% of saturation (19 kPa). Nymphs were then acclimated until the next morning, being free to eat and move between sides of the chamber. Each nymph experienced every oxygen-choice combination over the course of 5 days. Afterward, each nymph was lighted patted dry with paper towels (Kimwipes, Kimtech, Roswell, GA, USA) and weighed. New nymphs were then put in the chamber to experience each choice.

Temperature choice experiments

During the summer and fall of 2020, we measured how salmonfly nymphs distributed themselves between two static temperatures in an experimental aquarium (49×68×18.5 cm). The arena consisted of two insulated containers (24 Can Party Stackers, Coleman, Chicago, IL, USA) connected via half-inch PVC pipes, providing a route for movement between differing temperatures on either side. Containers were divided into five lanes via mesh walls. Additional mesh walls separated the experimental lanes from the air stones, stainless steel coils and water pumps that controlled the temperature, oxygen and flow conditions of each side of the chamber. High-grit sandpaper (Gator, P150) glued to the floor of the chamber provided footing for the nymphs.

We measured temperature preferences of 51 nymphs: 32 in normoxia (100% of air saturation; 19 kPa) and 19 in hypoxia (60–95% of air saturation; 11.4–18.05 kPa). At the onset of each trial, one nymph was placed on a random side of each lane of the aquarium and left to acclimate for 6–12 h. During this period, oxygen levels were ∼100% of saturation (19 kPa) on both sides of the aquarium and temperatures were either 15°C or 20°C. For treatments with a 15°C reference temperature, nymphs were given a choice the following morning between 15°C versus 15 (control), 20 or 25°C. These occurred at either normoxia (n=22) or hypoxia (n=19). For 20°C reference trials, nymphs were given a choice of 20°C versus 20 (control), 15, 25 and 30°C at normoxia only (10 nymphs). Nymphs were left in the chamber for 6–12 h, after which we recorded the location (warm or cool side) of each nymph and reset the temperature on each side of the aquarium to the reference temperature (15 or 20°C). Nymphs were then acclimated until the next morning. Each nymph was assigned to either the 15 or the 20°C reference experiment, but within that experiment they experienced every temperature-choice combination. Once nymphs experienced each choice trial, we weighed them and placed new nymphs in the chamber for a new round of trials.

Flow choice experiments

In the flow experiments, nymphs were subjected to ramps of decreasing flow, and we measured the water velocities at which they relocated to higher flows in an adjacent chamber. We conducted these experiments in winter 2021. The flow chamber (15×17×7.5 cm) was made of glued Plexiglas and divided lengthwise into two lanes by a median wall with a small gap (∼0.5 cm high) allowing nymphs to move between sides. High-grit sandpaper (Gator, P150) glued to the floor of the chamber provided footing for the nymphs. Variable flows on the two sides were generated by water pumps (SFBP1-G500, SeaFlo, Xiamen, China) positioned inline in closed loops (PVC) connected to either side of each lane of the chamber via milled nylon end caps (designed in Autodesk based on designs from Mike Nishizaki and milled at The Friday Harbor Laboratories by Adam Summers). Mesh walls inserted between the flow chamber and the end caps prevented nymphs from escaping into the PVC loops. Calibrated flow velocities on each side were controlled by a two-way variable DC power supply (HY3005F-3,RSR, Manheim, PA, USA) connected to the pumps.

We established how flow velocity varied as a function of supply voltage by filming neutrally buoyant microscopic glass particles (8–12 μm diameter, 1.5 g cm–3, TSI, Shoreview, MN, USA) moving through the chamber with a cellphone camera (Samsung Galaxy S10e, Suwon-si, South Korea). A ruler was placed within the chamber in the field of view of the camera during each recording. We used MATLAB (function: DLTdv7) to measure how many pixels spanned 1 cm on the ruler, the average pixels the glass beads moved per frame, and the frame rate of the video (MATLAB ver. R2019a). We then calculated the velocity of the water via the following equation: 1/(pixels/cm)×pixels/frame×frames/s=cm s−1. Flow velocity values used in the analyses were estimated by interpolating these data via linear regression (Fig. S2, Table S1).

During each experiment, we controlled temperature by submerging the apparatus into a temperature-controlled water bath. Water temperatures were monitored throughout each experiment by a Type-T thermocouple and meter (600-1020, Barnant, London, Ontario, Canada). Oxygen levels were controlled by bubbling air or mixtures of air and N2 gas into both lanes via small drilled holes on top of the chamber. Oxygen levels were monitored with an oxygen optode (Pyroscience, OXROB3) connected to a meter (Pyroscience, FireSting-O2) inserted through an additional hole.

To determine how temperature and oxygen interact with flow to impact microclimate preferences of nymphs, we performed flow choice experiments at two temperatures (10 and 18°C) and two oxygen conditions (100 and 60% of air saturation; ∼19.0 and 11.4 kPa). Ten nymphs were used in each treatment (40 nymphs total). At the onset of the experiment, a nymph was put into a random side of the flow chamber. Preliminary experiments suggested that 9 cm s−1 was the maximum velocity in which nymphs could move freely within the chamber without losing their footing and being swept into the mesh at the end. Preliminary experiments also showed that most nymphs explored both sides of the chamber for ∼40 min before settling into a single lane for at least 2 h. Based on these observations, nymphs were allowed to acclimate within the chamber for 50 min with flow velocities on both sides set to 9 cm s−1. Subsequently, we recorded the location (left versus right side) of the nymph and began to ramp down the flow velocity on that side of the chamber at 10 min intervals (7, 5, 3, 2, 1.5, 1, 0.75, 0.5, 0.25 and 0 cm s−1), while the flow velocity on the other side was kept constant at 9 cm s−1. The location of the nymph was recorded at each interval.

Statistical analyses

Before analysis, we removed extreme outliers from the microclimate dataset. This was done post hoc by pooling the measurements of each variable and discarding values that were <Q1–3×IQR or >Q3+3×IQR, where Q1 is the lower quartile, Q3 is the upper quartile, and IQR is the interquartile range. We used linear mixed-effects models implemented in the R package ‘nlme’ (function: lme) (https://CRAN.R-project.org/package=nlme) to test whether significant variation in water temperature, oxygen saturation and flow velocity existed among cobble microclimates (i.e. bottom-fronts, bottom-backs, and tops of cobbles and free stream) and along the length of riffles in the Blackfoot River and Rock Creek. Because we were not interested in conditions at particular sites and cobbles per se, site and cobble identity were modeled as random effects. We included the duration of time elapsed from the first sample to each sample per site as a covariate to account for changing conditions throughout the sampling period at each site.

To determine the effect of factor – oxygen, temperature or flow – on the variability of measurements within riffles and around individual cobbles, we used ANOVA models (function: lm) (https://www.r-project.org/). We tested the effect of factor on the coefficients of variation (standard deviation/mean) of measurements for each variable at the whole site and individual cobble scales. For each model, river was included as a covariate.

After performing the oxygen choice trials, we discarded data for which oxygen levels on the normoxic side fell below 85% of saturation (16.15 kPa), which was due to equipment malfunction. To analyze effects of temperature and oxygen on salmonfly choices, we used zero-inflated negative binomial mixed models (function: glmm.zinb) in the R package ‘NBZIMM’ (Zhang and Yi, 2020). More standard linear mixed-effects models could not be used because they resulted in model singularity, likely because the data in many cases were non-normal. NBZIMM uses a non-parametric approach that resolved this issue. Because choices of individual nymphs were measured repeatedly, we used nymph identity as a random effect. In some trials, we were unable to maintain desired temperature and oxygen gradients precisely enough within chambers. To account for this variation, we included the difference between choice conditions (e.g. warm-side temperature–cold-side temperature) on both sides of the chambers as a covariate. In addition, because nymphs in the temperature choice experiment experienced different reference temperatures (i.e. 15°C and 20°C), reference temperature was also included as a covariate. For both experiments, preferences for the left versus right side of the chamber in control trials were analyzed by performing exact binomial tests (function: binom.test) (https://www.r-project.org/).

To test the effects of flow, water temperature and oxygen saturation on choices nymphs made in flow experiments, we used a mixed-effects logistic regression model, a type of generalized linear mixed-effects model, in the R package ‘lme4’ (function: glmer) (Bates et al., 2015). Because individual nymphs were measured repeatedly, we included nymph ID as a random effect. We also included mass as a covariate. In addition, we analyzed the effects of temperature and oxygen using linear regression models (function: lm) (https://www.r-project.org/). We first analyzed whether temperature and oxygen levels influenced the water velocity at which nymphs switched between chambers for the final time (i.e. moved to a higher flow). However, some nymphs switched to the high flow chamber early during the ramp but later moved back and forth near the end of the trial. This resulted in several low ‘last switch’ flow velocity values, which may have skewed the results of the first model. To account for this potential bias, we also analyzed total time nymphs spent on the high flow side of the chamber.

For all of the above investigations, we did not pre-specify a target effect size. None of these experiments has been performed previously, and we thus did not know what the variances would be within groups – for either the observational data collected on abiotic conditions within streams or the experimental data on nymph choices made in response to sets of conditions in the lab. Rather, we strived to maximize the number of replicates in each observation or experiment given existing limitations on time and resources.

Field microclimate measurements

Oxygen saturation varied little at small scales within riffles of the Blackfoot River and Rock Creek (Fig. 2, Table 1). The mean range within riffles was <6% of air saturation (1.14 kPa) in both streams. In the Blackfoot, oxygen varied with distance and time elapsed (F1,111=5.34, P=0.023 and F1,111=21.90, P<0.001), but not micro-position. In Rock Creek, oxygen varied significantly with micro-position and time elapsed (F2,272=3.63, P=0.028 and F1,137=4.14, P=0.044, respectively), but not distance. Differences between micro-positions were slight even when significantly different (i.e. <0.4% of air saturation; 0.08 kPa) and are likely of little relevance to invertebrates, as shown below.

Water temperature varied modestly within stream riffles (Fig. 2, Table 1). The mean temperature range within riffles was ∼1.7°C in both streams. Temperature varied significantly with distance along riffles (Blackfoot: F1,111=36.19, P<0.001 and Rock Creek: F1,137=14.00, P<0.001), even after accounting for the effect sampling time on micro-thermal variation, which was significant in Rock Creek but not the Blackfoot (F1,137=8.36, P=0.005 and F1,111=1.43, P=0.234, respectively). On average, temperatures were ∼0.7°C warmer at riffle heads than tails, though they varied more strongly at some sites (e.g. Blackfoot River site 2 where riffle head ∼3.2°C warmer than riffle tail). Micro-position surrounding cobbles also had a significant effect on temperature (F2,235=5.07, P=0.007 and F2,276=6.39, P=0.002 for Blackfoot and Rock, respectively). However, mean thermal differences between micro-positions were small, particularly between the front and rear of cobbles (i.e.<0.06°C), and are likely of little biological importance, as demonstrated below.

Flow velocity varied strongly within stream riffles (Fig. 2, Table 1). The mean range of flow velocities at sites in both rivers were >125 cm s−1. Although flows did not vary along the lengths of riffles (F1,113=0.49, P=0.486 and F1,113=0.81, P=0.370 for Blackfoot and Rock Creek, respectively), flow velocity depended strongly on micro-position in both rivers (Blackfoot: F3,361=232.83, P<0.001; Rock Creek: F3,432=337.42, P<0.001). The lowest flows – mean ∼20 cm s−1 but often near 0 cm s−1 – occurred at the front and rear of rocks. Flows were usually high on the tops of cobbles and higher still in the free stream, with a mean of ∼64 and ∼87 cm s−1 for both rivers, respectively.

Interactions among covariates were largely non-significant for temperature, oxygen and flow data in each river and are therefore excluded from Table 1. See Table S2 for all the statistics, including interactions, of each microclimate model.

Factor identity had a significant effect on the coefficients of variation of measurements made within stream riffles (F2,187=155.03, P<0.001) and around individual cobbles (F2,852=90.26, P<0.001) (Table S3). The effect of river was non-significant for each model; interactions were non-significant. Oxygen, temperature and flow measurements had mean coefficients of variation of 0.01, 0.03 and 0.68, respectively, at the site scale for each river. At the individual cobble scale, oxygen, temperature and flow measurements had mean coefficients of variation of <0.01, <0.01 and 0.32, respectively, for each river.

Laboratory choice experiments

Oxygen saturation had a significant effect on the microdistributions of salmonfly nymphs in the oxygen choices experiments (F1,43=16.78, P<0.001) (Fig. 3A, Table 2). Choices were strongest when nymphs were exposed to steep oxygen gradients. 100% and ∼87% of nymphs were found on the normoxic side when given a choice of 10–35% (1.9–6.65 kPa) and 35–60% (6.65–11.4 kPa) saturation versus normoxia, respectively. However, when presented with more modest oxygen gradients of 60-85% (11.4–16.15 kPa) of saturation versus normoxia, only 45% of nymphs were found on the normoxic side. Initial explorations of the data showed that differences between the oxygen levels on both sides, temperature (Fig. S3A), and nymph mass (range: 0.06–1.32 g) had no significant effect on oxygen choices, and these covariates were thus excluded from the final model. In addition, there was no significant difference between nymph movements towards the left versus right side of the chamber in the control trials (P=0.832).

Temperature had a near-significant effect on distributions of nymphs in the temperature choice experiments (F1,77=3.55, P=0.063) (Fig. 3B, Table 2). There was a significant interaction between temperature and oxygen (F1,77=7.19, P=0.009), such that nymphs made stronger choices based on temperature when oxygen levels were low (Fig. S3B). However, the effect of this interaction was small, with an interaction coefficient of −0.01. Nymph mass (range: 0.10–1.56 g), difference between the temperatures on both sides, and reference temperature had no significant effect and were excluded from the final model. There was also no significant difference between nymph movements towards the left versus right side of the chamber in the control trials (P=0.850).

In the flow choice experiments, flow had a significant effect on the locations of nymphs (z=−4.08, P<0.001); nearly all nymphs (n=37/40) moved to the high-flow side before the low-flow side reached 0 cm s−1 (Fig. 4, Table 3). The effect of temperature was significant (z=2.07, P=0.039), such that nymphs had a stronger preference for high flows when temperatures were warm. In addition, there was a significant interaction between oxygen and flow (z=−1.99, P=0.047), such that nymphs relocated to higher flows more readily when oxygen levels were low. Mass (range: 0.18–1.40 g) had no effect and was dropped from the final model. Supplemental analyses yielded similar results: oxygen saturation (F1,36=6.16, P=0.018) and temperature (F1,36=5.45, P=0.025) had a significant effect on the time nymphs spent on the high flow side (Fig. 5A, Table 4); oxygen saturation had a significant effect on the water velocity of the low flow side at the final time nymphs moved from the low to the high flow side of the chamber (F1,36=11.00, P=0.002); and temperature had a near significant effect (F1,36=3.48; P=0.070) (Fig. 5B, Table 4). No significant interactions were found. Mass (range: 0.18–1.40 g) had no effect and was dropped from both models.

For insects, aquatic life is challenging because water contains little oxygen (Woods and Moran, 2020; Harrison et al., 2018, 2012; Verberk et al., 2011). This basic problem of oxygen scarcity can be further exacerbated by warm temperatures and low flows, which increase oxygen demand and decrease oxygen supply, respectively (Verberk et al., 2016a,b; Frakes et al., 2021). For insects, one potential solution is to exploit mosaics of local microclimates as a way of obtaining combinations of temperature, oxygen and flow that promote high performance or that minimize low-oxygen stress (Birrell et al., 2020). Whether insects do so, however, depends on (1) whether biologically meaningful variation in conditions exists at spatial scales accessible by individuals and (2) whether individuals can in fact detect and choose to exploit local variation. Here, we demonstrate that flow velocities in mountain streams varied much more strongly at small scales than did either temperature or oxygen. Accordingly, and in support of our hypothesis, nymphs exploited laboratory gradients of flows far more readily than they did temperature or oxygen; nymphs made choices on temperature and oxygen only when presented with unrealistically strong gradients, which do not reflect field observations.

We measured variation in temperature, oxygen and flow at two spatial scales in Montana streams – around individual cobbles (scale of ∼5 to 50 cm) and along the length of riffles (scale of ∼5 to 10 m). Differences in temperature and oxygen were slight around individual cobbles: <0.06°C and<0.4% of air saturation (0.8 kPa), respectively. Weak gradients such as these offer little opportunity to shift levels of oxygen supply and demand (Verberk et al., 2011) and are probably not actively exploited as evidenced by the results from our choice experiments. Within riffles, oxygen and temperature varied more substantially: oxygen saturation varied irregularly by ∼6% of air saturation (1.14 kPa), and temperature varied by ∼1.7°C, with water near the tail of the riffle being consistently cooler (∼0.7°C) than water near the head. These patterns likely arise because water downwells at riffle heads and is partially replaced by upwelling groundwater at the tails (Davy-Bowker et al., 2006). At several sites, temperatures changed more strongly with distance, being ∼3–4°C cooler at riffle trails. Strong, but inconsistent cold patches have been recorded in other streams, and are likely important refugia for fish (Ebersole et al., 2001, 2003). For stream insects, moving more than a few meters likely raises risks of dislodgement or predation. Indeed, in a rare mark–recapture study of salmonfly nymphs, individuals moved an average of only 1.8 m in 3 months (Freilich, 1991). In our system, movements at this scale would allow individuals to modify local temperatures and oxygen saturations by ∼0.31°C and 0.11% of saturation (0.02 kPa), respectively. Drifting can of course lead to much greater movement and may be important for exploiting distant microclimates (Lancaster and Downes, 2013).

Flow velocities, by contrast, varied substantially at small spatial scales, often by more than 125 cm s−1 around a single cobble. Local flows available to insects are likely even more variable still, as the flow meter we used was too large to fit underneath cobbles, where the lowest flows occur. Therefore, by moving only a few centimeters, individual insects should be able to strongly alter their local flow environments, from <1 cm s−1 under cobbles to >60 cm s−1 on the tops of rocks, decreasing boundary layer thicknesses and increasing rates of oxygen supply (Denny, 1993). Moving greater distances to exploit high flows (e.g. by drifting) is likely ineffective for increasing oxygen supply, as we found that flow velocities did not vary significantly along the lengths of riffles.

In support of our hypothesis that insects should evolve to exploit strong gradients in their local environments, salmonfly nymphs exploited laboratory variation in flow much more readily than they did temperature and oxygen. For example, the temperatures and oxygen levels that induced nymphs to make choices represented gradients that were ∼1–2 orders of magnitude greater than those observed in the field. In contrast, nymphs made consistently strong choices in response to ecologically relevant gradients of flow. This suggests that although differences in design between experiments – i.e. two static choices for temperature and oxygen versus a ramp for flow – may have had some influence on the results, nymphs nonetheless exploit local flows more readily than temperature or oxygen. Similar responses to laboratory flows have been observed in mayflies (Wiley & Kohler, 1980). In addition, we show some of the first experimental evidence that aquatic insects move to areas of high flow more readily, and spend more time in them, when temperatures are high and oxygen levels low. Our results corroborate field observations of caddisflies and stoneflies, showing that individuals occur more frequently in high-flow microhabitats when bulk flows were low or temperatures high (Kovalak, 1976, 1979; Genkai-Kato et al., 2005). Taken together, these results support the idea that interactions among temperature, oxygen and flow strongly influence aquatic ectotherm metabolism, performance and survival (Pörtner and Knust, 2007; Rubalcaba et al., 2020; Verberk et al., 2016a,b; Harrison et al., 2018; Frakes et al., 2021). In addition, they suggest that aquatic insects exploit mosaics of flow to increase ratios of oxygen supply to demand in nature. Curiously, however, mass had no effect on the choices made by nymphs, even though larger nymphs typically have greater demand for oxygen (e.g. Malison et al., 2022). We suggest that compensatory mechanisms may allow nymphs to maintain adequate rates of oxygen supply, e.g. by growing disproportionately larger gills or by increasing the density of tracheal tubes in their tissues in later instars.

As predicted from the weak gradients in temperature and oxygen measured at the scale of individual cobbles, nymphs made weak choices based on temperature and only when presented with ecologically unrealistic gradients (i.e. 20 versus 28°C). Nymphs made stronger choices based on oxygen, but as in the temperature experiments, only when gradients were stronger than those found commonly in nature (i.e. 10–55% versus 100% of air saturation; 1.9–10.5 kPa versus 19.0 kPa). Counter to our predictions, oxygen choices were not affected by temperature, possibly because 15°C was not warm enough to raise metabolic demand significantly (Malison et al., 2022). In addition, although hypoxia led nymphs to make stronger temperature choices, the effect of this interaction was weak (interaction coefficient −0.01). This is perhaps due to the apparently weak abilities of nymphs to sense and respond to temperature and oxygen. However, they clearly do have at least some ability exploit these conditions, which may allow them to avoid stressful conditions that occur irregularly in nature – likely avoiding deep, hypoxic pools or inflows from hot springs. In more usual circumstances, salmonflies and other aquatic insects may rely, at least in part, on their strong acclimation abilities to buffer against warming and hypoxia (Gunderson and Stillman, 2015).

In the face of climate change and other anthropogenic disturbances, organisms may go extinct, shift their ranges, evolve new capacities or tolerances, or employ different forms of plasticity (e.g. Donelson et al., 2019). Our study shows that behavioral plasticity – i.e. moving among locally available flows – may allow stream insects to mitigate oxygen limitation arising from climate change, river dewatering and damming, and nutrient-driven eutrophication (Birrell et al., 2020). Aquatic insects may also rely on other behaviors, such as more frequent ventilatory movements, to increase ratios of oxygen supply to demand (Genkai-Kato et al., 2000; Frakes et al., 2021). Other forms of plasticity – i.e. morphological and physiology plasticity – will likely also be important (Angilletta, 2009). For example, aquatic insects such as stoneflies and caddisflies tend to grow larger gills (Wichard, 1974; Badcock et al., 1987) and tolerate high temperatures better following acclimation to high temperatures or hypoxia (Gunderson and Stillman, 2015; Malison et al., 2022). Indeed, aquatic taxa have approximately twice the acclimation abilities of terrestrial taxa, which may help buffer them from future functional hypoxia from droughts, warming and other disturbances (Gunderson and Stillman, 2015).

We thank Michael Nishizaki and Adam Summers for help on constructing the flow choice chamber and Alisha Shah and Alex Birrell for helping measure microclimates in the field. We thank Tim Twardoski, Sean Kellogg and Elani Borhegyi for assistance in collecting choice data. We acknowledge James Frakes and Jack Hanson for helpful edits on the manuscript. We also thank the Montana Water Center for funding this research through the 2019 Montana Water Center Student Fellowship and the University of Montana for providing facilities for experiments. The paper was significantly clarified by constructive comments from two anonymous reviewers.

Author contributions

Conceptualization: J.H.B., H.A.W.; Methodology: J.H.B., H.A.W.; Formal analysis: J.H.B.; Investigation: J.H.B.; Resources: H.A.W.; Writing - original draft: J.H.B.; Writing - review & editing: J.H.B., H.A.W.; Visualization: J.H.B.; Supervision: H.A.W.; Project administration: H.A.W.; Funding acquisition: J.H.B., H.A.W.

Funding

This project was funded by a Montana Water Center Student Fellowship to J.H.B.

Data availability

All relevant data can be found within the article and its supplementary information.

Angilletta
,
M.
(
2009
).
Thermal Adaptation: A Theoretical and Empirical Synthesis
.
Oxford
:
Oxford University Press
.
Badcock
,
R. M.
,
Bales
,
M. T.
and
Harrison
,
J. D.
(
1987
).
Observations on gill number and respiratory adaptation in caddis larvae
. In
Proceedings of the Fifth International Symposium on Trichoptera
(ed.
M.
Bournaud
and
H.
Tachet
), pp.
175
-
178
.
Dordrecht
:
Dr W. Junk Publishers
.
Bates
,
D.
,
Mächler
,
M.
,
Bolker
,
B.
and
Walker
,
S.
(
2015
).
Fitting linear mixed-effects models using lme4
.
J. Stat. Soft.
67
,
1
-
48
.
Birrell
,
J. H.
,
Shah
,
A. A.
,
Hotaling
,
S.
,
Giersch
,
J. J.
,
Williamson
,
C. E.
,
Jacobsen
,
D.
and
Woods
,
H. A.
(
2020
).
Insects in high-elevation streams: life in extreme environments imperiled by climate change
.
Glob. Chang. Biol.
26
,
6667
-
6684
.
Clark
,
E.
,
Webb
,
B. W.
and
Ladle
,
M.
(
1999
).
Microthermal gradients and ecological implications in Dorset rivers
.
Hydrol. Process.
13
,
423
-
438
.
Collins
,
M.
,
Truebano
,
M.
,
Verberk
,
W. C. E. P.
and
Spicer
,
J. I.
(
2021
).
Do aquatic ectotherms perform better under hypoxia after warm acclimation?
J. Exp. Biol.
224
,
232512
.
Comiti
,
F.
,
Mao
,
L.
,
Wilcox
,
A.
,
Wohl
,
E. E.
and
Lenzi
,
M. A.
(
2007
).
Field-derived relationships for flow velocity and resistance in high-gradient streams
.
J. Hydrol.
340
,
48
-
62
.
Davy-Bowker
,
J.
,
Sweeting
,
W.
,
Wright
,
N.
,
Clarke
,
R. T.
and
Arnott
,
S.
(
2006
).
The distribution of benthic and hyporheic macroinvertebrates from the heads and tails of riffles
.
Hydrobiologia
563
,
109
-
123
.
Denny
,
M. W.
(
1993
).
Air and Water: The Biology and Physics of Life's Media
.
Princeton
:
Princeton University Press
.
Deutsch
,
C.
,
Ferrel
,
A.
,
Seibel
,
B.
,
Portner
,
H.-O.
and
Huey
,
R. B.
(
2015
).
Climate change tightens a metabolic constraint on marine habitats
.
Science
348
,
1132
-
1135
.
Donelson
,
J. M.
,
Sunday
,
J. M.
,
Figueira
,
W. F.
,
Gaitán-Espitia
,
J. D.
,
Hobday
,
A. J.
,
Johnson
,
C. R.
,
Leis
,
J. M.
,
Ling
,
S. D.
,
Marshall
,
D.
,
Pandolfi
,
J. M.
et al. 
(
2019
).
Understanding interactions between plasticity, adaptation and range shifts in response to marine environmental change
.
Philos. Trans. R. Soc. B Biol. Sci.
374
,
20180186
.
Ebersole
,
J. L.
,
Liss
,
W. J.
,
Frissell
,
C. A.
(
2001
).
Relationship between stream temperature, thermal refugia and rainbow trout Oncorhynchus mykiss abundance in arid-land streams in the northwestern United States
.
Ecol. Freshw. Fish
10
,
1
-
10
.
Ebersole
,
J. L.
,
Liss
,
W. J.
and
Frissell
,
C. A.
(
2003
).
Cold water patches in warm streams: physicochemical characteristics and the influence of shading
.
J. Am. Water Resour. Assoc.
39
,
355
-
368
.
Eggert
,
S. L.
and
Wallace
,
J. B.
(
2007
).
Wood biofilm as a food resource for stream detritivores
.
Limnol. Oceanogr.
52
,
1239
-
1245
.
Elliott
,
J. M.
(
2000
).
Pools as refugia for brown trout during two summer droughts: trout responses to thermal and oxygen stress
.
J. Fish Biol.
56
,
938
-
948
.
Frakes
,
J. I.
,
Birrell
,
J. H.
,
Shah
,
A. A.
and
Woods
,
H. A.
(
2021
).
Flow increases tolerance of heat and hypoxia of an aquatic insect
.
Biol. Lett.
17
,
20210004
.
Freilich
,
J. E.
(
1991
).
Movement patterns and ecology of Pteronarcys nymphs (Plecoptera) observations of marked individuals in a Rocky Mountain stream
.
Freshw. Biol.
25
,
379
-
394
.
Genkai-Kato
,
M.
,
Nozaki
,
K.
,
Mitsuhashi
,
H.
,
Kohmatsu
,
Y.
,
Miyasaka
,
H.
and
Nakanishi
,
M.
(
2000
).
Push-up response of stonefly larvae in low-oxygen conditions
.
Ecol. Res.
15
,
175
-
179
.
Genkai-Kato
,
M.
,
Mitsuhashi
,
H.
,
Kohmatsu
,
Y.
,
Miyasaka
,
H.
,
Nozaki
,
K.
and
Nakanishi
,
M.
(
2005
).
A seasonal change in the distribution of a stream-dwelling stonefly nymph reflects oxygen supply and water flow
.
Ecol. Res.
20
,
223
-
226
.
Gunderson
,
A. R.
and
Stillman
,
J. H.
(
2015
).
Plasticity in thermal tolerance has limited potential to buffer ectotherms from global warming
.
Proc. R. Soc. B Biol. Sci.
282
,
20150401
.
Hall
,
R. O.
and
Ulseth
,
A. J.
(
2020
).
Gas exchange in streams and rivers
.
WIREs Water
7
,
e1391
.
Harrison
,
J. F.
,
Woods
,
H. A.
and
Roberts
,
S. P.
(
2012
).
Ecological and Environmental Physiology of Insects
.
Oxford
:
Oxford University Press
.
Harrison
,
J. F.
,
Greenlee
,
K. J.
and
Verberk
,
W. C. E. P.
(
2018
).
Functional hypoxia in insects: definition, assessment, and consequences for physiology, ecology, and evolution
.
Annu. Rev. Entomol.
63
,
303
-
325
.
Helmuth
,
B.
,
Kingsolver
,
J. G.
and
Carrington
,
E.
(
2005
).
Biophysics, physiological ecology, and climate change: does mechanism matter?
Annu. Rev. Physiol.
67
,
177
-
201
.
Hutchinson
,
G. E.
(
1981
).
Thoughts on aquatic insects
.
Bioscience
31
,
495
-
500
.
Hynes
,
H. B. N.
(
1970
).
The Ecology of Running Waters
.
Liverpool
:
Liverpool University Press
.
Jacobsen
,
D.
(
2008
).
Low oxygen pressure as a driving factor for the altitudinal decline in taxon richness of stream macroinvertebrates
.
Oecologia
154
,
795
-
807
.
Jacobsen
,
D.
,
Rostgaard
,
S.
and
Vásconez
,
J. J.
(
2003
).
Are macroinvertebrates in high altitude streams affected by oxygen deficiency?
Freshw. Biol.
48
,
2025
-
2032
.
Jones
,
J. D.
(
1972
).
Comparative Physiology of Respiration
.
London
:
Edward Arnold
.
Kovalak
,
W. P.
(
1976
).
Seasonal and diel changes in the positioning of Glossosoma nigrior Banks (Trichoptera: Glossosomatidae) on artificial substrates
.
Can. J. Zool.
54
,
1585
-
1594
.
Kovalak
,
W. P.
(
1979
).
Day-night changes in stream benthos density in relation to current velocity
.
Arch. Hydrobiol.
87
,
1
-
18
.
Lancaster
,
J.
and
Downes
,
B. J.
(
2013
).
Aquatic Entomology
.
Oxford
:
Oxford University Press
.
Malison
,
R. L.
,
Frakes
,
J. I.
,
Andreas
,
A. L.
,
Keller
,
P. R.
,
Hamant
,
E.
,
Shah
,
A. A.
and
Woods
,
H. A.
(
2022
).
Plasticity of salmonfly (Pteronarcys californica) respiratory phenotypes in response to changes in temperature and oxygen
.
J. Exp. Biol.
225
,
244253
.
Matthews
,
K. R.
and
Berg
,
N. H.
(
1997
).
Rainbow trout responses to water temperature and dissolved oxygen stress in two southern California stream pools
.
J. Fish Biol.
50
,
50
-
67
.
Mérigoux
,
S.
and
Dolédec
,
S.
(
2004
).
Hydraulic requirements of stream communities: a case study on invertebrates
.
Freshw. Biol.
49
,
600
-
613
.
Mosley
,
M. P.
(
1983
).
Variability of water temperatures in the braided Ashley and Rakaia rivers
.
N. Z. J. Mar. Freshw. Res.
17
,
331
-
342
.
Pincebourde
,
S.
and
Woods
,
H. A.
(
2020
).
There is plenty of room at the bottom: microclimates drive insect vulnerability to climate change
.
Curr. Opin. Insect Sci.
41
,
63
-
70
.
Pinder
,
A. W.
and
Feder
,
M. E.
(
1990
).
Effect of boundary layers on cutaneous gas exchange
.
J. Exp. Biol.
154
,
67
-
80
.
Pörtner
,
H. O.
and
Knust
,
R.
(
2007
).
Climate change affects marine fishes through the oxygen limitation of thermal tolerance
.
Science
315
,
95
-
97
.
Rubalcaba
,
J. G.
,
Verberk
,
W. C. E. P.
,
Jan Hendriks
,
A.
,
Saris
,
B.
and
Woods
,
H. A.
(
2020
).
Oxygen limitation may affect the temperature and size dependence of metabolism in aquatic ectotherms
.
Proc. Natl. Acad. Sci. USA
117
,
31963
-
31968
.
Scheffers
,
B. R.
,
Edwards
,
D. P.
,
Diesmos
,
A.
,
Williams
,
S. E.
and
Evans
,
T. A.
(
2014
).
Microhabitats reduce animal's exposure to climate extremes
.
Glob. Chang. Biol.
20
,
495
-
503
.
Torgersen
,
C. E.
,
Price
,
D. M.
,
Li
,
H. W.
and
Mcintosh
,
B. A.
(
1999
).
Multiscale thermal refugia and stream habitat associations of chinook salmon in northeastern Oregon
.
Ecol. Appl.
9
,
301
-
319
.
Verberk
,
W. C. E. P.
and
Atkinson
,
D.
(
2013
).
Why polar gigantism and Palaeozoic gigantism are not equivalent: effects of oxygen and temperature on the body size of ectotherms
.
Funct. Ecol.
27
,
1275
-
1285
.
Verberk
,
W. C. E. P.
and
Bilton
,
D. T.
(
2015
).
Oxygen-limited thermal tolerance is seen in a plastron-breathing insect and can be induced in a bimodal gas exchanger
.
J. Exp. Biol.
218
,
2083
-
2088
.
Verberk
,
W. C. E. P.
,
Bilton
,
D. T.
,
Calosi
,
P.
and
Spicer
,
J. I.
(
2011
).
Oxygen supply in aquatic ectotherms: partial pressure and solubility together explain biodiversity and size patterns
.
Ecology
92
,
1565
-
1572
.
Verberk
,
W. C. E. P.
,
Overgaard
,
J.
,
Ern
,
R.
,
Bayley
,
M.
,
Wang
,
T.
,
Boardman
,
L.
and
Terblanche
,
J. S.
(
2016a
).
Does oxygen limit thermal tolerance in arthropods? A critical review of current evidence
.
Comp. Biochem. Physiol. Part A Mol. Integr. Physiol.
192
,
64
-
78
.
Verberk
,
W. C. E. P.
,
Durance
,
I.
,
Vaughan
,
I. P.
and
Ormerod
,
S. J.
(
2016b
).
Field and laboratory studies reveal interacting effects of stream oxygenation and warming on aquatic ectotherms
.
Glob. Chang. Biol.
22
,
1769
-
1778
.
Vogel
,
S.
(
1994
).
Life in Moving Fluids: The Physical Biology of Flow
, 2nd edn.
Princeton
:
Princeton University Press
.
White
,
D. S.
(
1990
).
Biological relationships to convective flow patterns within stream beds
.
Hydrobiologia
196
,
149
-
158
.
Wichard
,
W.
(
1974
).
Zur morphologischen Anpassung von Tracheenkiemen bei Larven der Limnephilini Kol. (Insecta, Trichoptera)
.
Oecologia
15
,
159
-
167
.
Wiley
,
M. J.
and
Kohler
,
S. L.
(
1980
).
Positioning changes of mayfly nymphs due to behavioral regulation of oxygen consumption
.
Can. J. Zool.
58
,
618
-
622
.
Woods
,
H. A.
(
1999
).
Egg-mass size and cell size: effects of temperature on oxygen distribution
.
Am. Zool.
39
,
244
-
252
.
Woods
,
H. A.
and
Moran
,
A. L.
(
2008
).
Temperature-oxygen interactions in Antarctic nudibranch egg masses
.
J. Exp. Biol.
211
,
798
-
804
.
Woods
,
H. A.
and
Moran
,
A. L.
(
2020
).
Reconsidering the oxygen-temperature hypothesis of polar gigantism: successes, failures, and nuance
.
Integr. Comp. Biol.
60
,
1438
-
1453
.
Woods
,
H. A.
,
Dillon
,
M. E.
and
Pincebourde
,
S.
(
2015
).
The roles of microclimatic diversity and of behavior in mediating the responses of ectotherms to climate change
.
J. Therm. Biol.
54
,
86
-
97
.
Zhang
,
X.
and
Yi
,
N.
(
2020
).
NBZIMM: negative binomial and zero-inflated mixed models, with application to microbiome/metagenomics data analysis
.
BMC Bioinform.
21
,
488
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information