Thermal optima in the hypoxia tolerance of marine ectotherms: Physiological causes and biogeographic consequences

The minimum O2 needed to fuel the demand of aquatic animals is commonly observed to increase with temperature, driven by accelerating metabolism. However, recent measurements of critical O2 thresholds (“Pcrit”) reveal more complex patterns, including those with a minimum at an intermediate thermal “optimum”. To discern the prevalence, physiological drivers, and biogeographic manifestations of such curves, we analyze new experimental and biogeographic data using a general dynamic model of aquatic water breathers. The model simulates the transfer of oxygen from ambient water through a boundary layer and into animal tissues driven by temperature-dependent rates of metabolism, diffusive gas exchange, and ventilatory and circulatory systems with O2-protein binding. We find that a thermal optimum in Pcrit can arise even when all physiological rates increase steadily with temperature. This occurs when O2 supply at low temperatures is limited by a process that is more temperature sensitive than metabolism, but becomes limited by a less sensitive process at warmer temperatures. Analysis of published species respiratory traits suggests that this scenario is not uncommon in marine biota, with ventilation and circulation limiting supply under cold conditions and diffusion limiting supply at high temperatures. Using occurrence data, we show that species with these physiological traits inhabit lowest O2 waters near the optimal temperature for hypoxia tolerance and are restricted to higher O2 at temperatures above and below this optimum. Our results imply that hypoxia tolerance can decline under both cold and warm conditions and thus may influence both poleward and equatorward species range limits.


Introduction
Climate change is raising temperatures throughout the upper ocean, while decreasing its oxygen content.These trends are among the most robustly observed and well-understood aspects of global ocean change [1].They also pose a major challenge for marine ectotherms, whose metabolic rates rise exponentially with temperature [2,3], requiring a concomitant increase in O 2 supply to maintain aerobic energy balance that is at odds with the ocean's declining global O 2 inventory [4,5].The temperature-dependent hypoxia tolerance of marine species already limits their geographic distributions, most commonly at the equatorward (warm) and/or deep (low O 2 ) range edge of species distributions [6][7][8][9], yielding a simple physiological mechanism for understanding species responses to climate change [10,11].
The minimum environmental O 2 at which an organism can sustain its resting metabolism is typically reported as a critical pressure (P crit ) and remains the most common measure of hypoxia tolerance, despite challenges associated with its experimental determination and interpretation [12][13][14][15].In most studied species, P crit increases with temperature, implying that their O 2 demand accelerates faster with warming than their rate of O 2 supply [6].Some species show a decrease in P crit as temperatures rise, implying that supply accelerates faster than demand, although this has rarely been observed [16,17].In recent experiments with higher than typical sample size over a broad temperature range, species exhibited both a decline in P crit as temperatures rise from the coldest water, followed by an increase from further warming, resulting in a minimum P crit , and thus a maximum hypoxia tolerance, at an intermediate optimum temperature [18,19].While the individual processes of supply and demand all tend to increase steadily with temperature [20,21], these bowl-shaped P crit curves require that the ratio of these rates exhibits a more complex relationship to temperature.
Thermal optima for hypoxia tolerance have been posited [22,23], but empirical support has been limited, in part because P crit is commonly estimated at too few temperatures, typically focused at the warm end of species tolerance ranges most relevant to anthropogenic climate change.These data limitations prevent the development of mechanistic and quantitative models capable of evaluating the role of hypoxia tolerance at the cold edge of range limits, and the associated implications of climate change, especially for populations not living near a species' warm range limit, or exposed to ocean cooling.
To examine the prevalence of thermal optima in hypoxia tolerance, diagnose the physiological conditions under which it can arise, and evaluate its relevance to species biogeography, we combined new laboratory experiments from species across multiple phyla, a dynamic model of O 2 supply in aquatic ectotherms, and species biogeographic distribution data.Among all studied species, we find complex P crit behavior across a broad range of temperatures.A general model of aquatic water breathers demonstrates the conditions under which thermal optima can emerge from the multistep nature of the O 2 supply chain, and analysis of published laboratory data suggests that marine species commonly meet those conditions.
The behavior of the dynamic model can be reproduced with a generalized model of temperature-dependent O 2 supply to demand ratios, termed "the Metabolic Index" [6], which captures a wide range of observed P crit curves.Finally, we present evidence that P crit curves with a thermal optimum are also reflected in the biogeography of marine species and thus may explain the cold limit in such species geographic ranges.

Laboratory observations
To evaluate the prevalence of nonexponential P crit curves, we measured P crit and metabolic rates across a wide temperature spectrum at high sample sizes for 4 previously unmeasured aquatic invertebrate species (Fig 1

Funding:
The study was made possible by grants to C.A.D. from the National Oceanic and Atmospheric Administration (NOAA NA18NOS4780167), the California SeaGrant and Ocean Protection Council, and the National Science Foundation (NSF OCE-1737282).E.A.S. is funded by the National Science Foundation (NSF EAR-1922966) and an Environmental Ventures Project grant from the Stanford Woods Institute.The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests:
The authors have declared that no competing interests exist.
we conducted P crit trials for the freshwater oligochaete worm Tubifex tubifex (n = 133), the outer shelf/upper slope sea urchin Lytechinus pictus (n = 39) from California, and the Atlantic intertidal anemone Nematostella vectensis (n = 107).We also used P crit measurements of groups of the squid Doryteuthis opalescens (n = 14 groups of 15 to 30 individuals), which is exposed to strong gradients of temperature and O 2 in the California Current System [24].To robustly estimate mean species P crit at each temperature, we performed the experiments on a much larger population of individuals (10s to >100) than in most previous studies.Full P crit and metabolic rate data are provided in the Supporting information (Fig A in S1 Text and S1 Data).Unprocessed O 2 drawdown curves can be found in https://doi.org/10.6084/m9.figshare.

24236257.v1.
Measurements of all species reveal complex P crit patterns across the studied temperature range (Fig 1A) even as metabolic rates increase exponentially with warming (Fig 1B).Despite the considerable variability and uncertainty in individual P crit estimates, the large number of measured individuals reveal a common pattern for all species.At intermediate temperatures, the mean P crit across individuals is lower than at warmer and colder temperatures, indicating a thermal optimum in resting hypoxia tolerance.The location, depth, and width of this optimum varies across species.For example, T. tubifex shows a deep bowl with a significant minimum in P crit at intermediate temperatures (Dunnett's test, comparing the 2 intermediate temperatures to the 2 colder and warmer temperatures, p = 0.00039 and p = 0.048, respectively).N. vectensis features a weaker but significant drop in P crit at 20˚C relative to colder temperatures (t test, p = 0.031), followed by a sharp rise at warmer temperatures.The broad bowl shapes seen in L. pictus and D. opalescens include thermal minima in P crit but are not statistically significant (Dunnett's test, p > 0.05).
These new respirometry data combined with published data [18,19] indicate that thermal optima are found in multiple phyla and across multiple modes of oxygen supply (e.g., gills and a blood vascular system in squid versus cutaneous respiration in anemones) and may therefore represent a widespread pattern.These patterns cannot be explained by the measured temperature sensitivity of metabolic O 2 demand alone, which rises throughout the temperature range (Fig 1B ), thus implying a key role for O 2 supply mechanisms.

Dynamic model
To explore the conditions that lead to a thermal optimum in hypoxia tolerance, we develop a dynamic model of O 2 supply and demand in water-breathing animals.The model simulates the transfer of O 2 from the environment to the metabolizing tissues of an organism across a range of temperatures using a system of coupled nonlinear ordinary differential equations (ODEs).To make the model generally applicable to aquatic animals, we include all the potential pools and fluxes of O 2 , including external ventilation of water from the ambient fluid to the boundary layer at the exchange surface, the molecular O 2 diffusion across that surface, and internal flux of O 2 to metabolizing body tissues, which may be mediated by a circulatory system (Fig 2A).In addition to dissolved O 2 , the model also tracks the concentration of the bound and unbound forms of an oxygen-transport protein such as hemoglobin or hemocyanin (denoted HxO and Hx, respectively), which bind and release molecular O 2 according to the associated chemical equilibrium.This is captured by the pO 2 at half-saturation (denoted P 50 ) and the enthalpy (ΔH) of the binding reaction, which governs the temperature dependence of that equilibrium [25] (S1 Text).
Each of the 3 O 2 supply processes (ventilation, diffusion, and circulation) is described by a rate S i that is represented as the product of the pO 2 difference between the respective compartments, Δ i pO 2 , and a temperature-dependent rate coefficient âi ðTÞ that characterizes the kinetics of that process: The 3 rate coefficients-flow rates of ventilated water and circulated blood and the diffusivity of O 2 -each vary exponentially with temperature (Arrhenius function), as does the metabolic rate, but with distinct temperature sensitivities.The resulting 8 parameters (3 supply rate coefficients, the metabolic rate, and the temperature sensitivity of each) along with the 3 chemical parameters (P 50 , ΔH, and total Hx + HxO concentration) represent a set of traits that determine a model organism's hypoxia tolerance and its variation with temperature.The welldocumented trait variations in real animals (e.g., overall O 2 supply capacity [26] and adaptation to hypoxia [27], gill surface area [28,29], and blood properties [30,31]) are simulated by scaling these parameters in the model.Our analysis aims to discern how such biological traits govern the shape of the resulting P crit curves with respect to temperature.
Model simulations resemble standard closed system respirometry experiments used to determine P crit values (Fig 2B) [32], in which O 2 is depleted from the ambient water as it gets transferred to metabolizing tissues.Both the O 2 concentrations and the fraction of O 2 -bound protein (HxO) decline in all compartments.Once O 2 levels in the tissue compartment can no longer support resting metabolism, consumption slows down with the onset of hypoxemia, allowing P crit to be diagnosed from the rate of environmental O 2 depletion using breakpoint analysis (Materials and methods; full model in S1 Text).
Simulations across a range of temperatures yield the P crit curve, which integrates the contribution of all traits to a single metric of hypoxia tolerance.Across a wide range of model parameters centered on the most common traits observed in marine organisms [7], the P crit curves exhibit an overall rise with temperature, driven by the increase in metabolic rate.Both the chemical properties (Hx + HxO, ΔH, and P 50 ) as well as the rate coefficients of supply and demand ( âi ) have qualitative impacts on the P crit curves that are intuitive.For example, a higher concentration of total oxygen transport protein (Hx + HxO) acts to lower the P crit curve across all temperatures (Fig 2C ), enhancing the tolerance to hypoxia.An equivalent effect can The shape of the P crit curves also exhibits a sensitivity to parameters that is quantitatively less intuitive, as the fractional change in P crit is not always the same across the full temperature range (Fig 2 ).For example, a given increase in the ventilation rate does not lower the P crit by the same fraction at all temperatures but instead has a larger impact under cold conditions than under warm conditions (Fig 2D).In other words, the P crit curves resulting from a multistep supply chain can depart from simple exponential relationships with temperature, even when each single supply process accelerates exponentially with warming.We conclude that the well-known nonlinearities in blood-O 2 binding are not the essential cause of this behavior, because the variation due to biophysical properties is similar to that induced by variations in blood chemistry and complex P crit curves are observed in organisms without O 2 -binding proteins (e.g., N. vectensis in Fig 1; [17]).Instead, we focus our analysis on the mechanisms by which the linear combination of biophysical transfer processes in a multistep O 2 supply chain leads to the complex patterns observed in P crit curves.
The origins of nonexponential P crit curves can be demonstrated quantitatively in a model with a supply chain consisting only of ventilation and diffusive gas exchange (Fig 3).In isolation, each step yields a simple (exponential) P crit curve with a slope depending on the temperature sensitivities of supply and demand.The curve is increasing if metabolic demand accelerates faster with temperature than supply (shown for diffusion;  The additive nature of the P crit curve resulting from a linear supply chain can also be derived analytically from the system of model ODEs for more than 2 supply steps (S1 Text).Conceptually, this property can be thought of as analogous to an electrical circuit in which a fixed voltage is applied to a series of resistors.Just like the total voltage can be obtained as the sum of the individual voltage drops across each resistor, the total P crit curve of a multistep supply chain can be obtained as the sum of the pO 2 drops that drive each individual supply process.
A bowl-shaped P crit curve can emerge if the supply chain includes processes that are both more and less sensitive to temperature changes than metabolism.In Fig 3C, the P crit curve rises under warm conditions because a large pO 2 gradient is required to drive sufficient diffusion at high temperatures.This is due to the fact that diffusion accelerates slower than metabolism with warming.On the other hand, the curve also remains flat or even reverses under cold conditions because a large pO 2 gradient is required to provide sufficient O 2 transport via ventilation at low temperatures, since this process has a higher temperature sensitivity than metabolism and thus slows down more strongly in cold water.
Because the critical pO 2 differences required to drive the individual supply steps are not the same, the total P crit curve is not equally sensitive to changes in the biologically controlled rate coefficients at all temperatures.In the example above, the change in P crit at high temperatures due to a change in ventilation rate might be small or even negligible, while its response to a change in diffusivity might be substantial, even for the same relative increase in the biologically controlled parameter.More generally, a change in the coefficient of any supply process that accelerates faster with warming than metabolism will have the largest impact on P crit under cold conditions, as in the case of ventilation.On the other hand, such an increase has the largest impact on P crit under warm conditions for a supply process that accelerates slower than metabolism, such as diffusion.This relationship is particularly important for processes under immediate biological control like ventilation and circulation and has implications for understanding their temperature sensitivity.Incurring the energetic costs of accelerating heart rate or ventilation across the entire temperature range may not be beneficial if P crit is instead much more sensitive to changes in diffusion at high temperatures.We illustrate this in a model variant with a ventilation rate that has a high temperature sensitivity at low temperatures but reaches an upper limit under warm conditions, as, for example, observed in [33][34][35].The resulting change in P crit at high temperatures compared to a simple exponential ventilation rate is minimal (Fig C in S1 Text).In this scenario, increasing the ventilation rate throughout the warm side of the temperature range barely impacts hypoxia tolerance, because O 2 supply is largely determined by diffusion.

Evidence from physiology
To determine whether the physiological conditions for a thermal optimum in hypoxia tolerance are common among marine biota, we compiled experimental data on the temperature dependence of ventilation and circulation rates of aquatic water breathers (Materials and methods).The compilation covers 58 data sets from 35 species, including 21 chordates, 9 arthropods, 3 annelids, and 2 mollusks.Estimates of the temperature sensitivity parameters (E V , E C ) that control the slope of the exponential relationship with temperature are obtained by fitting Arrhenius functions to the data.The results show an increase in ventilatory and circulatory activity with temperature in almost all species (i.e., E V , E C > 0).The estimated sensitivities range from −0.14 eV to 0.9 eV and have a mean of 0.39 eV (±0.22 eV SD) that is much greater than the theoretical sensitivity of diffusion (E D = 0.06 eV), which is the O 2 supply process that is present in all organisms.Since these estimates include ventilation and circulation frequencies (e.g., heart rates) and stroke volumes (e.g., liters per heartbeat) in addition to actual volumetric flow rates (i.e., liters per min), they represent a lower bound on the sensitivity of (volumetric) ventilation and circulation rates as considered in the dynamic model.A higher sensitivity (0.49 eV ± 0.21 eV SD) is obtained if only volumetric rates (n = 12) are considered (Fig 4A).
The biophysical responses to temperature can be compared to existing estimates of the temperature sensitivity of metabolism (E M ) with a mean of 0.71 eV [±0.46 eV SD] from a diverse set of 186 species [7].If the traits of O 2 supply and demand are considered independent, the estimated frequency distributions in Fig 4A predict that the conditions for thermal optima (i.e., E D < E M < E V , E C ) are met in about 23% of species after accounting for the effect of species) fall between the theoretical prediction for the sensitivity of diffusion (vertical orange line) and published estimates for the sensitivity of metabolic rates (n = 186 species; data from [7]) on average but with significant overlap, indicating the widespread potential for thermal optima in P crit .Lines show kernel density estimates of the trait frequency distributions.(B) Example for the estimation of E V from published data across the full inhabited temperature range of the polychaete worm Nereis succinea (solid black; [33]) as well as under cold and warm conditions only (blue dashed and red dotted lines, respectively).(C) Estimated temperature sensitivities E V , E C of volumetric ventilation and circulation rates from published data under cold (blue) and warm (red) conditions as illustrated in panel (B).The data underlying panels (A) and (C) can be found in S1 Data.The data underlying panel (B) can be found in [33].
https://doi.org/10.1371/journal.pbio.3002443.g004decreasing solubility with temperature (S1 Text).However, this estimate likely represents a lower bound for the occurrence of bowl-shaped P crit curves, for 2 reasons.First, supply sensitivities exceed that of demand in 7 of 17 species for which both estimates are available (S1 Data and Fig F in S1 Text).Thus, about 40% of species with paired E V,C and E M data meet this condition for having a thermal optimum in hypoxia tolerance.
Second, the thermal increase in ventilation and circulation rates is stronger on the cold side of the inhabited temperature range for many species (Fig 4B ), where it is more likely to limit metabolism.We estimated the temperature sensitivity of volumetric flow rates in both warm and cold temperature ranges for all species with sufficient data (n = 10; Fig 4C).On the cold side, the mean E V and E C are 0.69 eV, significantly exceeding the warm side with a mean of 0.07 eV (p = 0.009).The difference remains significant if stroke volumes and frequencies are also considered (n = 40 from 25 species).Thus, the acceleration of ventilation and circulation rates slows down in the warmer part of the temperature range on average.This behavior is consistent with biophysical processes conferring little additional hypoxia tolerance at high temperatures for the associated energetic cost, with P crit being most sensitive to diffusion under warm conditions, as illustrated in the model variant for ventilation (Fig C in S1 Text).It also means that interspecies estimates of E V , E C for comparison to E M (Fig 4A) are likely lower than they would be if measured at the cold edge, where these biophysical processes limit metabolism.Taken together, the interspecies distributions of E V , E C , and E M (Fig 4A ), and the higher temperature sensitivity of biophysical rates (circulation and ventilation) in colder waters (Fig 4B and 4C), including the reduction of these sensitivities in warmer waters, suggest that the conditions for thermal optima are commonly found among the traits of marine species.

A Metabolic Index with thermal optima
We generalized the Metabolic Index introduced by Deutsch and colleagues [6] to account for the occurrence of complex shaped P crit curves.The index is defined as the ratio of O 2 supply to resting demand of an aquatic water breather in its environment, and it has been applied to understand how species biogeography is shaped by climate [6,10,11,18,36].However, the original formulation assumed that P crit varies exponentially with temperature.Our generalized version is able to reproduce the full range of behaviors exhibited by the dynamic model, requiring only 5 parameters, which can be calibrated from experimental data through a single equation.
The generalized Metabolic Index can be derived analytically from the dynamical model (S1 Text) but can also be developed heuristically from the electrical circuit analogy, noted above.The metabolic demand requires an equal rate of O 2 flow (a "current") to sustain it.Every step in the supply chain acts like a resistor with a biologically determined and temperature-dependent resistance whose inverse (a "conductance") is given by the rate coefficient âi of the process.The flow at every step is the product of the conductance and the pO 2 difference ΔpO 2 (a "voltage") between the compartments, whichAU : Pleasecheckandconfirmthattheeditsto}Theflowateverys drives the current (Eq 1).Each step is associated with such a required voltage drop-a pO 2 difference-determined by its single step conductance.The steps in the supply chain, from ventilation to diffusion and circulation, act as resistors wired in series, and the P crit of the composite supply chain is the minimum "voltage" required to achieve an O 2 supply matching demand and is the sum of the minimum pO 2 differences of the single supply steps, as illustrated in Fig 3C .The temperature-dependent rate coefficient (or "conductance") âi of a single supply step can be expressed as a i � RðT; E S i Þ, where α i denotes the value of the coefficient at reference temperature, which is scaled by an exponential (Arrhenius) function R with temperature sensitivity E S i [eV].More generally, in a chain with n supply steps in series, the total conductance of the chain is the reciprocal of the sum of single step resistances.When divided by metabolic demand α M �R(T, E M ), the resulting expression for the generalized supply-to-demand ratio F is where the α i represent the supply rate coefficients at reference temperature, and [eV] denote the differences between the sensitivities of metabolic demand and the supply processes.The dependence of supply and demand on body mass B is reflected in the allometric exponent � as in the original index [6].
The condition for the existence of a bowl-shaped P crit curve, i.e., supply steps having temperature sensitivities both less than and greater than that of metabolic demand, thus reads E S i < E M < E S j for any 2 supply steps i and j.Eq 2 can include any number of supply processes.However, we find that P crit curves generated by the full model (n = 3) can still be reproduced by curves assuming only 2 steps (Fig D in S1 Text).The generalized Metabolic Index in Eq 2 can reproduce P crit curves that include the Hx/HxO system (Fig E in S1 Text), because the effects of the chemical blood component on the P crit curve are qualitatively the same as those of the biophysical parameters (details in S1 Text).Adding more exponential curves also does not change the qualitative range of possible P crit curves beyond those of concave-up bowl shapes (see Figs D and E in S1 Text).In such cases, the parameters can no longer be associated with single steps in the supply chain but instead capture the combined properties of the processes that limit the O 2 supply toward the cold and warm ends of the temperature range, respectively.The asymmetry of the bowl-shaped curves can vary among species, depending on the values of the temperature sensitivities.

Connecting physiology to biogeography
The Metabolic Index framework establishes a direct link between physiological traits and species distributions, as the range boundaries of a diverse set of species align more strongly with a lower threshold of O 2 supply to demand predicted by the index than with either temperature or pO 2 alone [7].The generalized formulation has the potential to further improve this description of species habitats, especially at the cold edges of a species distribution.
To examine whether the thermal optima in physiological hypoxia tolerance are reflected in a species' biogeography, we investigate state-space habitats of biogeographic occurrence data from the Ocean Biodiversity Information System [37] (see Materials and methods).For the species presented in Fig 1, the aerobic habitat conditions (T and pO 2 ) are poorly represented in large-scale datasets (Fig G in S1 Text).However, for 2 additional species with physiological traits suggesting a thermal optimum, the starry flounder Platichthys stellatus and the shrimp Oplophorus spinosus, adequate occurrence and environmental data are available.In Platichthys stellatus, estimates from published experimental results yield a temperature sensitivity E M = 0.68 eV for metabolism and E V = 0.9 eV [38] for the ventilation rate, indicating a bowl-shaped O 2 limitation.For Oplophorus spinosus, critical O 2 pressures have been measured and display a minimum at intermediate temperatures [39], such that Eq 2 can be fit directly.
For both species, the environmental conditions in occupied habitats reveal a clear minimum in inhabited pO 2 at intermediate temperatures, consistent with the measured and predicted physiological thresholds (Fig 5).In contrast, the minimum inhabited temperatures of each species are inconsistent with a model based on a lower threshold value of temperature that is independent of O 2 .Instead, minimum temperatures decrease to lower values as oxygen levels increase.Similar patterns are also observed in other species for which laboratory experiments on metabolic and O 2 supply traits indicate thermal optima and for which sufficient occurrence data are available (details in S1 Text and Fig H in S1 Text; [18]).
In all these cases, the generalized Metabolic Index reveals how the reversal in hypoxia tolerance at low temperatures results from physiological traits, and how this bidirectionality is reflected in biogeographic ranges.In particular, it suggests that O 2 limitation is the mechanism that restricts habitat toward the cold edges of species distributions.

Discussion
The dynamic model of temperature-dependent hypoxia reveals that a series of biophysical O 2 supply steps can give rise to thermal optima in hypoxia tolerance as observed in new high- resolution respirometry data.This occurs when the supply chain includes at least 2 processes such that one accelerates with temperature more slowly than metabolic demand, and another accelerates more rapidly.In this case, the process with a lower temperature sensitivity drives an increase in P crit under warm conditions, while the more sensitive process leads to a reversal with higher P crit in cold waters.A generalized Metabolic Index adequately captures these complex patterns in a single metric based on mechanistic principles.
Our analysis of available physiological evidence suggests that such bidirectional effects of temperature on hypoxia tolerance may not be uncommon in aquatic animals across taxonomic groups.Currently available estimates of the temperature sensitivity of ventilation and circulation rates in aquatic ectotherms imply the existence of thermal optima in a significant fraction of species because many estimates of E V /E C exceed observations of E M , especially in cold waters where biophysical supply mechanisms appear most important.Sampling the involved traits across a broader range of the taxonomic, morphological, and ecological diversity is a key step toward further advancing and testing this framework and its implications, as there are only a few teleost and crustacean species for which all required physiological estimates are available.Our observation that the acceleration of ventilation and circulation rates slows down under warm conditions (Fig 4B and 4C and Fig C in S1 Text) also needs to be investigated more closely, as we only obtained 10 datasets that allowed this estimation.However, the causes and consequences of this phenomenon can also be investigated further theoretically by incorporating the metabolic cost incurred due to ventilation and circulation rates in the dynamic model.In particular, such a cost analysis should include the temperature-dependent viscosity of water and blood, which are not included here.
Our results highlight new opportunities to diagnose the role of O 2 limitation in limiting the geographic range and body size of marine species.Patterns of body size under cold conditions in polar fish has recently been attributed to limitations on O 2 supply [40].The generalized Metabolic Index provides a framework to test this body size hypothesis, given additional measurements of P crit across the range of body size and temperatures among these taxa [41].In the environment, species minimum inhabited pO 2 levels are elevated above physiological thresholds by the additional O 2 requirements needed to sustain metabolic activity above the resting state, equivalent to the sustained metabolic scope [7,42], with the potential to influence T-pO 2 relationships based on physiology alone.However, in case studies presented here, species' biogeographic state space strongly reflected thermal optima in physiological hypoxia tolerance (Fig 5).In contrast to the sparsity of detailed physiological measurements, global occurrence data are available for a much larger number and diversity of marine species (e.g., in OBIS).Future work should aim to clarify and quantify such potentially synergistic physiological and ecological mechanisms driving the dependence of thermal tolerance on O 2 supply and demand.
Oxygen limitation of aerobic metabolism at low temperature also has broad implications for marine ecosystems and their response to climate change [43].Marine species richness is generally observed to decline toward the poles and equator from a peak in midlatitudes [44,45] and is often cited as being driven by gradients in ocean temperature, with coldest and warmest waters taken to inhibit diversity.Our results indicate that long-term aerobic energy constraints on viable habitat in cold water could be a physiological cause of this poleward diversity decline, just as aerobic constraints in warm water can explain the equatorial diversity dip [43].At the same time, warming at species' poleward range limits would relieve such aerobic constraints, allowing species to disperse toward, and establish in, higher latitudes.This mechanism could thus potentially explain widespread poleward range shifts of marine species seen in response to recent anthropogenic warming [45][46][47].On longer timescales, O 2 limitation at species' cold edge habitat limits provides a novel mechanism for driving habitat loss and extinction risk during periods of global cooling [48,49].

Ethics statement
The studied invertebrates were collected, reared, and handled according to all applicable permits and ethical standards of treatment.In particular, collection of the cephalopod D. opalescencs was carried out under California Department of Fish and Wildlife (CDFW) Scientific Collection Permit SC-13563, and husbandry and experiments were conducted under Institutional Animal Care and Use Committee (IACUC) #10643.

Laboratory measurements
Critical O 2 levels were measured following standard closed system respirometry protocols [17,50] for individuals of T. tubifex (n = 132), N. vectensis (n = 107), and L. pictus (n = 40).For the social squid D. opalescens, we measured critical O 2 levels for 14 groups of 15 to 30 (median 20) animals following published closed system respirometry protocols for this species [51].P crit was determined by breakpoint analysis of the O 2 draw down curve [52].Full protocols are provided in S1 Text.

Dynamic model
The pools and fluxes of O 2 in a generic water breather are described by a nonlinear system of 8 ordinary differential equations.For each set of model parameters, simulations are performed across the temperature range from 0˚C to 30˚C until P crit can be determined by breakpoint analysis from the rate of O 2 draw down.All simulations were carried out in the Python language using the solve ivp function in Scipy [53] for numerical integration.The full model description is provided in S1 Text.

Ventilation and circulation data
We compiled data on ventilation rates (n = 8), ventilation frequency (n = 18), ventilation stroke volumes (n = 6), circulation rates (n = 4), heart rates (n = 20), and heart stroke volumes (n = 2) of aquatic water breathers measured at 2 or more temperatures at atmospheric O 2 levels.Estimates of the sensitivity parameters E V , E C were obtained through least square fits of Arrhenius functions to the experimental data using the curve fit function and density estimates were obtained using the gaussian kde function in Scipy.A detailed description of the compiled data is provided in S1 Text, and all estimates are available in S1 Data.

State-space habitats
State-space habitats were obtained by pairing species location data downloaded from the Ocean Biodiversity Information System [37] on June 6, 2022, with monthly temperature and O 2 conditions from the World Ocean Atlas [54,55] according to the procedure described in [7].All available state-space habitats are shown in Figs G and H in S1 Text.Additional information is also provided in S1 Text. of ventilation, circulation, and metabolic rates of aquatic water breathers from published results.(XLSX) and Fig A in S1 Text).These species span 4 different phyla, have multiple modes of oxygen supply, and regularly encounter temporal or spatial variability in environmental pO 2 .Following published respirometry protocols (Materials and methods), metabolic rates underlying Fig 4 as well as Fig F in S1 Text were obtained from published results and can be found in S1 Data.All modeling and analysis code underlying the remaining figures can be found in https://doi.org/10.6084/m9.figshare.24547255.

Fig 1 .
Fig 1. Temperature-dependent critical oxygen pressures (P crit ) of 4 marine invertebrate species display thermal optima even as metabolic demands increase exponentially.(A) New closed system respirometry experiments exhibit a minimum of P crit (mean ± SE) at intermediate temperatures, indicating a thermal optimum in hypoxia tolerance.The species include an oligochaete worm (Tubifex tubifex), a sea urchin (Lytechinus pictus), an anemone (Nematostella vectensis), and a cephalopod (Doryteuthis opalescens).(B) Normalized mass-specific metabolic O 2 demand (mean ± SD) follows a simple Arrhenius relationship in all species and cannot explain the complex patterns observed in P crit .Metabolic rate data are plotted relative to the rate at the coldest measured temperature within a species.See text for details.The data underlying this figure can be found in S1 Data.Untreated experimental data can be found in https:// doi.org/10.6084/m9.figshare.24236257.v1.https://doi.org/10.1371/journal.pbio.3002443.g001

Fig 2 .
Fig 2. Overview of the model used to investigate the effects of a multistep O 2 supply chain.(A) The model tracks the concentrations of O 2 as well as unbound and bound O 2 transporting proteins ("Hx," "HxO") in 4 compartments representing external and internal volumes of water or body fluid, which are connected through a linear O 2 supply chain with external ventilation, diffusion, and internal circulation.(B) A model run at a single temperature resembles a closed system respirometry experiment.The saturation of O 2 (solid colored) and proportion of HxO (dashed colored) decline in all compartments until the O 2 level in the metabolizing tissue (dark red) reaches a critical limit near zero, at which point metabolic consumption slows down and P crit (dashed black) can be determined from the rate of environmental O 2 depletion.(C) Effects of increasing the concentration, half-saturation pressure (P 50 ), or temperature sensitivity (ΔH) of O 2 transport protein on the P crit curve.(D) Effects of increasing the rate coefficients of biophysical supply and demand processes.A higher metabolism elevates the curve, while increasing the rate of any supply process lowers it.No experimental data underlie this figure.Code underlying this figure can be found in https://doi.org/10.6084/m9.figshare.24547255.https://doi.org/10.1371/journal.pbio.3002443.g002 Fig 3A), and decreasing if instead the temperature sensitivity of supply exceeds that of metabolism (ventilation; Fig 3B).Combining ventilation and diffusion in series results in a P crit curve that is the sum of the 2 curves corresponding to the single steps and thus exhibits a minimum at an intermediate temperature (Fig 3C), similar to the observations.

Fig 3 .
Fig 3. Modeled thermal optima in P crit curves.At any given temperature, the model P crit can be found analytically as the intersection (black dots) of the demand isocline (dashed colored) and the corresponding supply isocline (solid colored), aligning with the curve diagnosed from numerical simulations (solid black).(A) In a model with only diffusive gas exchange characterized by a smaller temperature sensitivity than metabolic demand, the isocline intersections yield an increasing P crit curve.(B) Conversely, the steeper supply isoclines lead to a decreasing pattern in a single supply step model with ventilation that accelerates faster than metabolic demand.(C) Combining the 2 supply steps in series results in a P crit curve that is the sum of the single step curves (dotted black), giving rise to a thermal optimum at intermediate temperatures.No experimental data underlie this figure.Code underlying this figure can be found in https://doi.org/10.6084/m9.figshare.24547255.https://doi.org/10.1371/journal.pbio.3002443.g003

Fig 4 .
Fig 4. Temperature sensitivities of ventilation and circulation rates estimated from published experimental data.(A) The estimates (blue, n = 58 from 35species) fall between the theoretical prediction for the sensitivity of diffusion (vertical orange line) and published estimates for the sensitivity of metabolic rates (n = 186 species; data from[7]) on average but with significant overlap, indicating the widespread potential for thermal optima in P crit .Lines show kernel density estimates of the trait frequency distributions.(B) Example for the estimation of E V from published data across the full inhabited temperature range of the polychaete worm Nereis succinea (solid black;[33]) as well as under cold and warm conditions only (blue dashed and red dotted lines, respectively).(C) Estimated temperature sensitivities E V , E C of volumetric ventilation and circulation rates from published data under cold (blue) and warm (red) conditions as illustrated in panel (B).The data underlying panels (A) and (C) can be found in S1 Data.The data underlying panel (B) can be found in[33].

Fig 5 .
Fig 5. Biogeography and state-space habitats of species with thermal optima in hypoxia tolerance.Global occurrence data for 2 species, (A) the midwater shrimp Oplophorus spinosus and (B) the flounder Platichthys stellatus, from the Ocean Biodiversity Information System are mapped on the climatological temperature (lines) and pO 2 (color field), from the World Ocean Atlas, averaged from the surface to the 95th percentile depth of reported species occurrences.The location (latitude, longitude, and depth) and month of each reported occurrence are matched to the climatological T and pO 2 , whose joint histogram is shown in the state-space habitat diagrams for each species.(C) In Oplophorus spinosus, measured P crit values (red dots) as well as the curve fit based on the Metabolic Index (dashed red line) align with the lowest inhabited pO 2 across the temperature range [10].(D) In the flounder Platichthys stellatus, the P crit curve (dashed red line) predicted from the Metabolic Index framework and physiological rates also exhibits a thermal optimum consistent with the occurrence data.The data underlying this figure can be found in https://doi.org/10.6084/m9.figshare.24547270.The code underlying this figure can be found in https://doi.org/10.6084/m9.figshare.24547255.https://doi.org/10.1371/journal.pbio.3002443.g005