Different drivers, common mechanism; the distribution of a reef fish is restricted by local-scale oxygen and temperature constraints on aerobic metabolism

The metabolic index (ϕ) captures how temperature and oxygen interact to limit marine organism distributions at local scales. These findings indicate that ϕ has utility as a conservation tool for local management.


Introduction
Climate change is driving the redistribution of species across the globe (Poloczanska et al., 2013;Pecl et al., 2017). Marine species are particularly responsive to climate-driven distribution shifts because they are exposed to environmental conditions closer to their limits and have fewer colonization barriers than their terrestrial counterparts (Sunday et al., 2012;Pinsky et al., 2019). The general warming trend of the ocean has, on a broad scale, already driven the distribution of fishes to higher latitudes (Perry et al., 2005) or deeper depths (Dulvy et al., 2008), but this directional pattern is far from ubiquitous (Pinsky et al., 2020). Predicting the direction and magnitude of climate-driven distribution shifts of marine organisms is important to guide management and conservation strategies (Brander, 2007;Link et al., 2011) but is complicated by variable, species-specific distribution responses to prevailing climatic conditions (Pinsky et al., 2020). Identifying the fundamental physiological mechanisms governing how organisms respond to environmental variability may overcome these complications and improve predictive accuracy (Evans et al., 2015).
Ectothermic organism performance is regulated by the environment through its influence on internal physiological rates and processes (Bozinovic and Pörtner, 2015;Horodysky et al., 2015). All physiological and biochemical processes where energy is consumed or stored are expressed as an organism's metabolism (Nelson, 2016). These energy conversions are fuelled by ATP (adenosine triphosphate) which is generated up to 30 times more efficiently through aerobic versus anaerobic metabolic pathways (Richards, 2009;Clarke, 2019). Any constraints to an organism's aerobic metabolism thus correspond to reductions in overall energy available to provision processes related to performance. When localized environmental conditions change, organisms can seek areas where rates of aerobic metabolism can be maintained, ultimately resulting in a net distribution shift (Walther et al., 2002;Parmesan, 2006;Cheung et al., 2009). Determining the environmental conditions that constrain aerobic metabolism is therefore an important physiological consideration that can improve interpretability and forecasting accuracy of distribution shifts (Kearney and Porter, 2009;Teal et al., 2018;Rodríguez et al., 2019).
Temperature and oxygen availability are the primary controlling and limiting factors on the aerobic metabolism of ectotherms (Claireaux and Chabot, 2016) and are also the two environmental variables most pervasively influenced by anthropogenic climate change (Gruber, 2011). Temperature regulates the rate of aerobic metabolism as more energy is required for physiological processes driven faster by increasing kinetic energy at high temperatures (Clarke and Fraser, 2004). Oxygen availability poses an upper limit on aerobic metabolism, as the rate at which an organism can consume the required oxygen to fuel physiological processes is set by rates of diffusion via a pressure gradient from the environment across gill epithelium (Ern, 2019). The temperature effect on aquatic ectotherm aerobic metabolism is often depicted as absolute or factorial aerobic scope [the difference or ratio between the maximum (MMR) and standard metabolic rates (SMR) respectively] and is thought to represent the excess metabolic energy above that required for maintenance at a given temperature (Fry, 1971;Farrell, 2016). The oxygenand capacity-limited thermal tolerance (OCLTT) hypothesis posits that the often-bell-shaped aerobic scope-temperature relationships are driven by differences in oxygen supply and oxygen demand and define the thermal envelope within which an organism can persist (Portner and Farrell, 2008). There is, however, mounting evidence against the universality of the OCLTT hypothesis, including mismatches between thermal optima for aerobic scope and various important ecological and physiological processes (Gräns et al., 2014;Norin et al., 2014;Verberk et al., 2016). Whilst the OCLTT concept suggests that thermal limits are set by failures in the oxygen delivery system, it does not explicitly incorporate how this system is influenced by ambient oxygen availability, which together with temperature are the primary variables governing metabolism. A quantitative model that can evaluate how both these environmental variables interact to drive an aquatic organisms' aerobic metabolism has only relatively recently been developed, in the form of the metabolic index (φ) (Deutsch et al., 2015).
The metabolic index (φ) represents the ratio of oxygen supply to resting oxygen demand, incorporating the effect of both temperature and oxygen availability. Oxygen demand is quantified as the critical oxygen partial pressure (pO 2crit ) which represents the minimum level of oxygen availability to sustain a standard metabolic rate (Ultsch and Regan, 2019) and oxygen supply is taken as the prevailing oxygen partial pressure (pO 2 ). The critical pO 2 for maximum metabolic rates of aquatic organisms has evolved to match prevailing pO 2 and thus at normoxia φ is akin to factorial aerobic scope (MMR/SMR) (Seibel and Deutsch, 2020). Unlike factorial aerobic scope, φ quantitatively incorporates oxygen availability together with temperature effects on aerobic metabolism and offers a powerful framework to explore how variability in these two environmental variables can control the distribution of viable habitat for a species (Somero et al., 2016). Indeed, Deutsch et al. (2015) showed that a mean φ threshold (φ crit ) ranging between 2 and 5 limits the equatorward distribution of four marine species from diverse habitats. Penn et al. (2018) argue that historical temperature and oxygen levels at the end of the Permian Period exceeded φ crit , accounting for the spatial variation of mass marine extinctions at the time. These studies explored spatial trends in φ across broad north to south extents and from coarse-resolution ocean models, with the patterns observed largely attributed to increasing temperature driving metabolic demand above supply.
At local scales, organisms respond to abrupt environmental heterogeneity rather than broader climate signals, which occur on scales of metres to kilometres, particularly around  . For example, localized upwelling can pull deeper oxygen-deprived water towards the surface resulting in anomalously low coastal oxygen partial pressures and sea temperatures compared to surrounding coastal waters (Lutjeharms et al., 2000;Grantham et al., 2004). Here, we use the considerable spatial variability in temperature and oxygen availability along South Africa's coastal zone (Roberts, 2005) as the model system to explore if φ can explain the contemporary distribution patterns of an endemic marine fish (Chrysoblephus laticeps) and predict future distribution responses. We hypothesize that φ will be limited towards (1) the warm-temperate edge of C. laticep's distribution through a temperature-induced increase in oxygen demand and (2) the cool-temperate edge through the low oxygen availability in this region.

Study area and species profile
The South African coastal zone is characterized by a significant local-scale spatial variability in temperature and oxygen availability signals making it a model system to explore how these two important environmental variables may interact to limit species distributions. The coastal zone is subdivided into three main bioregions (cool-temperate, warm-temperate and subtropical) with these characterized by spatial and temporal heterogeneous ocean temperatures and oxygen availability (Potts et al., 2015). The east coast is subtropical, with coastal temperatures driven by the warm Agulhas current, which has intensified and warmed since the 1980s (Rouault et al., 2009). Oxygen concentration is high (>4 ml.l −1 ) relative to temperatures throughout most parts of the east coast (Pretorius et al., 2016). By contrast, the south coast is warm-temperate and is characterized by temperature variability patterns driven by intermittent upwelling (Goschen and Schumann, 2011). Although there is a negative trend in mean annual sea surface temperature along areas of the south coast associated with increases in upwelling (Rouault et al., 2010), some localized areas of warming are also reported (Lima and Wethey, 2012). Oxygen availability is relatively high throughout the south coast (>4 ml.l −1 ), but there are localized areas (towards the west of the south coast) where oxygen levels can get to critically low (0-2 ml.l −1 ) levels (Roberts, 2005). The southwest coast is cool-temperate, with environmental patterns dominated by the more permanent Benguela upwelling cell that is situated north of Cape Town. Sea temperatures are cooler than the south coast, and oxygen concentrations are low (0-2 ml.l −1 ) (Roberts, 2005;Jarre et al., 2015). Here, the Benguela upwelling has strengthened since the 1970s resulting in a negative mean annual sea surface temperature trend (Santos et al., 2012).
We used the roman seabream (Chrysoblephus laticeps) as a model species for this study. Typical of most species in the sparid genus, C. laticeps is considered slow growing, has a maximum age of 19 years (Götz et al., 2008), attains 50% maturity as a female between ages 2.5 and 4.27 (Buxton 1987, 1990, Götz et al. 2008, and undergoes a protogynous sex change (Buxton, 1990) between ages 8 and 10.25 (Götz et al., 2008). Importantly, this species is highly resident, with the probability of being recaptured within the Tsitsikamma National Park Marine Protected Area estimated at 0.94 (Kerwath et al., 2007a) and a home range size estimated using acoustic telemetry of between 1 and 3 km 2 (Kerwath et al., 2007b). Because Chrysoblephus laticeps forms an important component of the commercial fisheries in South Africa, a high-resolution spatial data of where the species is caught available (Kerwath et al., 2013).

Chrysoblephus laticeps
The metabolic index (φ) (Equation 1) represents the ratio of oxygen supply to oxygen demand where supply is the prevailing oxygen partial pressure (pO 2 ) and demand is the minimum pO 2 required to sustain a standard metabolic rate (pO 2crit ) taking into account temperature and mass (Deutsch et al., 2015).
The parameters of φ (A o and −E o ) are derived from the slope and intercept of the linear relationship between critical oxygen partial pressure (pO 2crit ) which has been mass standardized [divided by mass (B) and its scaling exponent (n)] and the inverse of k B T [the product of temperature (T) in kelvin and the Boltzmann constant (k B )]. These parameters are species-specific and thus need to be derived from laboratory experiments that quantify a species' pO 2crit across a temperature range.
Fifty C. laticeps specimens were caught from either the Tsitsikamma Marine Protected Area (n = 25, date = 23 September 2016) or offshore of the Noordhoek Ski-boat Club outside of Port Elizabeth (n = 25, date = 20 September 2016) and brought back alive to the Aquatic Ecophysiology Research Platform laboratories of the South African Institute for Aquatic Biodiversity. Fish were left to adjust for ∼6 weeks to holding tank conditions kept at the average bottom sea temperature of sampling locations (16 • C) before standard and maximum metabolic rates were quantified. These metabolic rate measurements formed part of a complimentary study quantifying the metabolic response of C. laticeps to abrupt temperature swings which are a common occurrence through its distribution (see Duncan et al. 2019 for detailed methods). The protocol involved placing fasted individuals inside respirometers at holding temperatures (16 • C), allowing ∼12 h to adjust to respirometer conditions and then subjecting them to an acute (1 • C per hour) temperature change to either 8, 12, 16 (no change), 20 or 24 • C test temperatures. At these test temperatures, intermittent flow respirometry (15-min flush, 5-min measure) was run for ∼20 h whereafter individuals were chased to exhaustion, exposed to air and placed back inside the respirometer to elicit maximum metabolic rate. The standard metabolic rate was determined as the quantile that assigned the bottom 20% of all metabolic rate measurements prior to the chase protocol . Following the elicitation of maximum metabolic rate each individual was left in its respirometer and given time to recover (∼4 to 5 h) to metabolic rates approximating its standard metabolic rate. At this point, the protocol to determine the critical oxygen saturation (O 2crit ) was started. Respirometer flush pumps were turned off, and oxygen was depleted within a respirometer by the organism until distinct signs of stress or impaired performance were observed. Metabolic rates were measured in 5-min intervals during this progressive hypoxia and combined with all earlier metabolic rate measurements to quantify the critical oxygen level.
O 2crit was defined as the oxygen saturation level where standard metabolic rates could no longer be maintained and below which metabolic rates decreased in proportion with prevailing oxygen saturation levels. To quantify O 2crit , we used the 'calcO2crit' function, developed by Claireaux & Chabot (2016), where O 2crit is taken as the point of intersection between standard metabolic rate and the linear declines in metabolic rates below the standard metabolic rate in hypoxic water ( Figure S1.1 in the Supporting Information). First, all metabolic rate measurements (mg O 2 .min −1 .kg −1 ) from Duncan et al. (2019) were combined with metabolic rate measurements during progressive hypoxia (this study) and paired with the average oxygen saturation of respirometer water during each measurement period. The 'calcO2 crit ' function was run for each individual dataset with the maximum number of points for the regression below the standard metabolic rate set to 15 to limit the potential for including too many points above 'true O2 crit ' and overestimating. Subsequent O2 crit plots were visually inspected, and O2 crit was considered overestimated if a large number of points were included in the regression below SMR, but metabolic rates were not decreasing in proportion to oxygen saturation levels. This can happen if metabolic rate measurements greater than standard metabolic rate are lacking at low oxygen levels and, in such cases, the maximum number of points was limited to four (one more than the minimum required for a regression). All O2 crit thresholds (% saturation) were converted to corresponding critical oxygen partial pressures (pO 2crit , kPa).
To standardize mass effects on pO 2crit , we estimated the mass scaling exponent (n) for C. laticeps by first temperature standardizing data (done by dividing with the Arrhenius function) and finding the best fit mass scaling exponent from the mass-pO 2crit power-law relationship (Deutsch et al., 2015). To test for differences in pO 2crit between the Tsitsikamma (TNP) and Port Elizabeth (PE) sampling areas, the linear relationship between pO 2crit and temperature was modelled including area as an additive effect and its interaction with temperature, but only for temperatures between 12 and 24 • C to ensure the residuals of the linear model were normally distributed. We also followed the same approach to determine if sex had a significant effect on pO 2crit . The sex and stage of each study organism were determined by inspecting the gonads of euthanized individuals following the O 2crit trial. Without microscopic gonad sections, it was difficult to determine the sexual stage (male, female or intersex) with 100% accuracy. Sex was therefore classified as either female (F), male (M), predominately male but possibly intersex stage (M/I) or predominately female but possibly intersex stage (F/I), and included as an additive and interactive effect with temperature in a linear model with pO 2crit as the dependent variable . If no mass and sex stage effects were found, all data were pooled to estimate the parameters of φ.

Environmental and occurrence data
Spatial environmental data were obtained from a highresolution (0.25 • horizontal resolution) coupled oceanbiogeochemistry model (Yool et al., 2015;Popova et al., 2016) for the periods between 2005 and 2009 (contemporary) and projections towards 2095-2099 (future). The model framework is built on the physical Nucleus for European Modelling of the Ocean (NEMO) model (Madec, 2008) coupled with a Model of Ecosystem Dynamics, Nutrient Utilisation, Sequestration and Acidification (MEDUSA 2.0) (Yool et al., 2013). The forward projection is based on the IPCC representative concentration pathway 8.5, which is the worst case, business-as-usual scenario of future greenhouse gas concentrations (IPCC, 2014).
Mean monthly contemporary and future oxygen concentration (mmol.m −3 ), temperature ( • C), salinity (psu) and depth data were extracted from the layer of cells closest to the sea floor because C. laticeps is a benthic reef-associated species. Data were extracted for depths between 0 and 100 m below sea level and from around the southern African coast (26.5-37.5 • S and 13.5-37.5 • E). Oxygen concentrations were converted to partial pressure (kPa) based on grid cell temperature, oxygen concentration, salinity and pressure at depth using the 'o2_unit_conv' function in the 'presens' package (Birk, 2016). Pressure at depth was independently determined based on grid cell depth and latitude using the 'swPressure' function in the 'oce' package (Kelly and Richards, 2018). The mean of surrounding cells was used to estimate values for those cells with no data near the coast and the resolution of environmental data was increased to 0.125 • horizontal resolution (13.31-11.27 km) by locally resampling using bilinear interpolation with the 'raster' package (Hijmans, 2016). Monthly φ grid cells were generated by projecting the calibrated φ equation across the spatial and temporal extent of temperature and oxygen partial pressure data from the ocean model. contemporary (2005)(2006)(2007)(2008)(2009)) and future (2095-2099) model predictions. Gridded bathymetry data were obtained from the General Bathymetric Chart of the Oceans (GEBCO) website at a 30 arc-second resolution and resampled to a 0.125 • horizontal resolution (13.31-11.27 km) by taking the mean depth of every cell within the 0.125 • grid and clipped between 0 and 100 m below sea level.
Occurrence points for C. laticeps were obtained from the National Marine Linefish System (NMLS) for the period between 2000 and 2010, made available by The Department of Environment, Forestry and Fisheries [DEFF, formally the Department of Agriculture Fisheries and Forestry (DAFF)]. The NMLS contains mandatory species level catch and effort data from the boat-based commercial fisheries operating in South African waters since 1985 and is one of the largest geo-referenced marine datasets in the world (DAFF, 2016). The boat-based fishery operates throughout South Africa's coastal zone, and we are confident that the occurrence points from the NMLS accurately represent where C. laticeps occurs.

Distribution modelling
We used the random forest algorithm (Breiman, 2001) to model the current distribution, quantify critical φ levels (φ crit ) and predict future distributions of C. laticeps up to 2099 using the 'randomForest' package (Liaw and Weiner, 2002). Random forests are an extension of traditional classification trees, where many trees are grown using bootstrapped samples of data, and a random selection of predictor variables at each tree node is used for classification (Cutler et al., 2007). We used random forests because they are not especially sensitive to predictor collinearly or the distribution of data and are commonly used and considered accurate for species distribution modelling (Cutler et al., 2007;Mi et al., 2017).
We followed the guidelines of Barbet-Massin et al. (2012) for applying random forests to model the distribution of C. laticeps. A dataset consisting of the known occurrence points of C. laticeps was combined with an equal number of pseudo-absence data points. Pseudo-absence data points were generated by randomly selecting cells throughout the extent of the predictor variables but excluding cells adjacent to known occurrence points (Barbet-Massin et al., 2012). This presencepseudo-absence dataset was then randomly subdivided into a model training (80%) and a model testing (20%) dataset. The predictive performance of each model was assessed based on the accuracy of correctly classifying the test dataset and variable importance of the predictors ranked according to the mean decrease in accuracy based on random permutations of each predictor variable.
To assess what value of φ may limit the distribution of C. laticeps (φ crit ), 10 full random forest models were trained and assessed with maximum, minimum and mean monthly φ, and depth as predictor variables, for the contemporary period. Each of the 10 models was trained with a unique generation of pseudo-absence data and subsequent random training/testing data split. Ten more reduced random forest models were then re-trained and tested using only the top two predictor variables from the full models. The predictive accuracy (measured as a proportion of correctly identified test data cells) of the full versus the reduced models were compared and the simplest combination of predictor variables with the highest accuracy subsequently used to project current and future distributions of C. laticeps. Current and future distributions of C. laticeps were estimated by combining all 10 model projections, trained on unique combinations of presence-pseudo-absence data (described above), with most confidence given to areas where all 10 models agreed. To quantify the critical threshold of φ (φ crit ), the feature contribution (the contribution of a data point to the final model prediction) of each training data point was first extracted from the random forest models using the 'forestFloor' package (Welling et al., 2016) for all 10 current model runs and combined. A feature can contribute either positively (towards an occurrence prediction) or negatively (towards an absence prediction) to the random forest model prediction. In order to fit a logistic regression through the φ feature contribution data, they were converted to binary data by reclassifying negative feature class contributions to zero and positive contributions to one. φ crit was estimated as the φ value where the logistic regression model = 0.5. Other predictor variable thresholds were estimated as the points where the majority of feature class contributions transitioned from positive to negative.

Comparing φ as a predictor variable to other common metrics
To compare where and how φ may differ from distribution predictions over other commonly used predictor variables, we followed the distribution modelling steps outlined above but for distribution models including minimum, maximum and mean values of either absolute aerobic scope (AAS), factorial aerobic scope (FAS) or prevailing temperature as predictor variables. Aerobic scope environmental layers were generated by fitting quadratic polynomials through the relationship between mass standardized standard (SMR) or maximum metabolic rates (MMR) and temperature, from data corresponding to the same individuals in this study but measured in Duncan et al. (2019), and projecting these models through the temperature layers of the ocean model. Absolute aerobic scope was taken as MMR minus SMR for a given temperature grid cell and FAS taken as MMR divided by SMR. For each predictor variable (AAS, FAS or temperature), we ran 10 random forest models with minimum, mean and maximum values from the transformed ocean model as predictor variables and added model outputs together resulting in spatial distribution predictions ranging from 100 (all 10 models agree for the area to be suitable) to 0 (no models indicate area to be suitable).

pO 2crit experiments
Of the 50 specimens obtained from the wild, we ran 39 (20 from PE and 19 from TNP) O 2crit trials across the temperature range tested (Table S1.1 in the Supporting Information). Specimen mass ranged from 0.32 to 1.55 kg and was evenly distributed across temperature treatments. O 2crit appeared to be overestimated in only 5 of the 39 trials and was adjusted to four regression points accordingly (Table S1.1 in the Supporting Information). Of the 39 specimens, we identified 10 as males, 8 as female, 16 as females but possibly intersex and 3 as males but possibly intersex, and 2 were not identified (Table S1.1 in the Supporting Information). Sex classifications were evenly distributed among temperature treatments ( Figure S1.4 in the Supporting Information).

Calibration of the metabolic index
We identified a weak positive relationship between mass and pO 2crit , with the mass scaling exponent estimated at 0.17 ( Figure S1.2 in the Supporting Information). pO 2crit data revealed a positive linear relationship with temperature (12-24 • C), which was not significantly different between sampling areas and their interaction with temperature (Tables S1.3, P-value > 0.05, Figure S1.3 in the Supporting Information) or sexual stage and its interaction with temperature (Tables S1.3, P-value > 0.05, Figure S1.4 in the Supporting Information). All data were subsequently pooled to determine the parameters of φ for C. laticeps. The differing slopes of the pO 2crit -temperature relationship (8-12 and 12-24 • C) permitted fitting a piecewise relationship to estimate φ parameters on either sides of 12 • C [ Fig. 2a, 12 • C = 40.69 (1/k B .T)]. The estimated φ model parameters A o, E o , taken from the intercept and slope respectively, could thus be estimated for each linear relationship. At temperatures greater than 12 • C, we estimated A o to be 1.55951E+9 and E o to be −0.4885, and below 12 • C A o was estimated to be 5.109275E−18 and E o to be 1.01. A graphical representation of this calibrated piecewise φ for C. laticeps is presented across a matrix of temperatures ( • C) and oxygen levels (% saturation), indicating how temperature or oxygen availability can limit φ (Fig. 2b).

Environmental and occurrence data
A total of 14 934 occurrence points were obtained and reduced to 206 points once duplicates were removed from the 0.125 • grid cell resolution. Occurrence points were distributed throughout the South African south coast, in accordance with the published distribution of C. laticeps (Fig. 1)  Mean contemporary (2005)(2006)(2007)(2008)(2009)) spatial bottom temperature (Fig. 3a) and oxygen partial pressure (kPa) (Fig. 3b) were variable throughout the extent of the modelling domain. Mean current bottom temperatures ranged from 11.9 to 27.7 • C and were the coolest along South Africa's south coast but increasing in the east and west towards higher latitudes (Fig. 3a). The lowest monthly temperature of 10.33 • C occurred in a cell off the south west coast in the Benguela upwelling zone, and the highest monthly temperature of 31.2 • C occurred in a cell towards the north eastern extent of the environmental data. Mean contemporary oxygen availability ranged from 13.6 to 21.5 kPa, with low oxygen zones occurring offshore of South Africa's west and central south coasts (Fig. 3b). The lowest monthly oxygen availability was 6.4 kPa and occurred along the south-west coast in the Benguela upwelling zone. The spatial variability in temperature and oxygen resulted in a heterogeneous spatial distribution φ for C. laticeps throughout the modelling domain (Fig. 3c). The mean φ for the period between 2005 and 2009 ranged from 2.17 to 4.78 and was the greatest along the south coast as warm temperatures along the east coast and low oxygen availability along the west coast limited φ.

Species distribution modelling with φ
The full random forest models based on all magnitudes (minimum, maximum and mean) of contemporary φ and depth had a mean predictive accuracy of 0.93 across all 10 runs. Depth and minimum φ were the two most important predictor variables in all full-model runs. The mean predictive accuracy of the reduced model with just depth and minimum φ as predictor variables was 0.92 and was statistically indistinguishable from the full models (student T test, P-value > 0.05, Table S1.4 in the Supporting Information). The reduced models were used to predict distribution changes and φ crit for C. laticeps. Extracting contributions from parameters of the reduced models indicated that suitable depths for C. laticeps occur between around 13 and 75 m below sea level ( Fig. 4a), and the threshold of minimum φ (φ crit ) below which C. laticeps can no longer persist is 2.78 (Fig. 4b).
The contemporary modelled distribution indicated a core range for C. laticeps from Cape Town in the west towards Port St. Johns in the east, with a break and a small patch of suitable habitat in the east (Fig. 5a). Projections indicate that this core distribution will persist up until 2100, but the western and eastern edges of the species distribution will contract slightly (Fig. 5b).

Species distribution models with AAS, FAS and temperature
The distribution models with different sets of predictor variables [minimum, maximum and mean of either factorial aerobic scope (FAS), absolute aerobic scope (AAS) or temperature] predicted a similar core distribution for C. laticeps to the φ model ( Figure S1.5 in the Supporting Information). The absolute aerobic scope (AAS) models however were unable to adequately delimit the warm edge of C. laticeps' distribution with more than 50% of random forest models predicting suitable habitat to occur in Mozambican tropical waters that does not support the robust occurrence data [ellipses, Figure S1.5 (d) in the Supporting Information]. Around the cooler, low oxygen western edge of C. laticeps' range, the factorial aerobic scope (FAS), absolute aerobic scope (AAS) and temperature only models predicted distributions to stretch past Cape Town into areas where catches of C. laticeps are not common (Fig. 6a-d). The φ distribution model most accurately delimited the spatial area where C. laticeps abundantly occurs across both edges of its known distribution, which was further supported by the φ distribution model iterations having a higher mean predictive accuracy over models with other predictor variables ( Figure S1.6 in the Supporting Information). This difference in predictive accuracy among sets of distribution models was significant [one-way ANOVA, φ or temperature predictor models (Tukey's HSD, adjusted P-value < 0.05).

Discussion
The geographic distribution of C. laticeps is constrained by a minimum critical metabolic index (φ crit ) threshold of 2.78, whereby φ is limited in the west by supply (low oxygen) and in the east by demand (high temperatures). Importantly, we show that the distribution of C. laticeps is limited by different environmental variables across its contemporary distribution but through their influences on a common mechanismaerobic metabolism.
Identifying consistent patterns in physiological limits on species distributions across taxa is rare (Bozinovic et al., 2011;Evans et al., 2015). The monthly φ crit value for C. laticeps (2.78) found here; however, is similar to the summer φ crit (3.2) found for the related sharpsnout seabream (Diplodus puntazzo) in the north Atlantic (Deutsch et al., 2015). Deutsch et al. (2015) do show consistent patterns between taxa regarding the warm equatorward limit of φ, with summer φ crit ranging from 2.2 for common eelpout (Zoarces viviparus) to 3.3 for Atlantic cod (Gadus morhua), due to the antagonistic effect of rising temperature on φ through exponentially increasing oxygen demand. These results suggest that, on average, marine fishes require environmental conditions that permit a factorial aerobic scope of around three, although the number of fishes tested is small. Estimating in situ metabolic rates of wild fishes (field metabolic rates) remains a technological challenge (Treberg et al., 2016). Recent advances in otolith isotope analysis have revealed that the average field metabolic rate (wild fishes) to standard metabolic rate (measured in the laboratory) ratio for Atlantic cod, G. morhua, ranged between 1.9 and 3.5 (Chung et al., 2019). This isotope field metabolic rate work suggests that, on average, the wild metabolic rates of G. morhua will begin to be compromised if environmental conditions limit factorial aerobic scope below 3.5 which is comparable to annual φ crit of 3.7 for G. morhua distribution (Deutsch et al., 2015). The φ crit will vary per species by lifestyle and morphology, but for low-energy teleost fishes that are not obligate swimmers, the aerobic supply to demand ratio of ∼ 2-3.5 appears to be an important constraint (Deutsch et al., 2015).
A bidirectional pattern between hypoxia tolerance and temperature has largely gone undetected, perhaps due to a disproportionate focus on warm-temperature experiments in an era of climate warming (Szekeres et al., 2016). This pattern is new to fishes but has been reported for an intertidal Anthozoan cnidarian, Diadumene lineata (Boag et al., 2018). A bidirectional relationship between hypoxia tolerance and temperature (e.g. Fig. 2a, 12 • C) indicates that O 2crit is not primarily driven by temperature effects on increased oxygen demand (standard metabolic rate) but rather by the oxygen supply capacity and how it scales with temperature compared to the maximum metabolic rate (Seibel and Deutsch, 2020). The decrease in oxygen supply capacity at cold temperatures may be driven by the increased viscosity and thicker boundary layers of cold water that can impede oxygen diffusion rates (Verberk et al., 2011). When such a bidirectional pattern is present, the utility of φ is extended beyond predicting only warm edge habitat limits and towards including both warm and cool edge limits across the entire distribution of an organism. In this case, however, there were few temperature points below 12 • C associated with minimum annual φ for locations where C. laticeps occurs indicating that low oxygen associated with lower temperatures is primarily limiting φ (Fig. 7). When the temperature and oxygen availability associated with minimum φ layers where C. laticeps occurs are mapped, it mirrors the pO 2crit -temperature bidirectional relationship. This correspondence highlights the interactive effects of both variables in setting the suitability of a habitat patch (Fig. 7). Given the increased utility of φ when a bidirectional O 2crit -temperature relationship is present, we advocate for generating O 2crit data across the full range of temperatures experienced by an organism and at rates of exposure that are ecologically relevant, for future experimental physiology.
Biodiversity is affected by multiple dimensions of climate change across space and time that may not be captured in metrics of magnitude such as means (Garcia et al., 2014;Waldock et al., 2018). We find that the minimum monthly φ best explained the distribution limits of C. laticeps rather than means or maximums in all 10 random forest model runs. Whilst we rely on magnitudes (mean, maximum, minimum) of φ for this analysis, the monthly temporal resolution gives an indication of φ temporal availability (the total duration of a specific environmental event [e.g. Waldock et al., 2018)] that C. laticeps can tolerate. Whilst it is possible that monthly temporal resolution is too coarse to quantify absolute lower thresholds of φ, we did observe good recoveries of C. laticeps following acute (over 48 h) temperature and oxygen stressors in the lab, suggesting a capacity to withstand sub-optimal φ on a daily scale. The dimension of 'time' is an important consideration for future lab studies to calibrate φ as longer-term acclimation can reduce the temperature sensitivity of metabolic measurements (Semsar-kazerouni and Verberk, 2018;Bates and Morley, 2020). It is important that the rates of change in environmental stressors in the lab and the temporal resolution of spatial layers are ecologically relevant to accurately quantify organism responses to prevailing environmental change. At best, we can conclude that C. laticeps requires temperature and oxygen conditions that result in an average factorial aerobic scope of ∼three, over monthly periods.
The effect of temperature on species distributions is well studied in the marine environment (e.g. Lasram et al., 2010;Albouy et al., 2013), but since oxygen is more difficult to measure in the field and manipulate in the laboratory, it has been neglected, despite its importance in limiting aerobic energy production in organisms (Bianucci et al., 2016;Stortini et al., 2017). For example, Verberk et al., (2016) consider both variables and show that in situ oxygen availability influences the thermal suitability of a habitat patch for two Mayfly species, which corresponds to laboratory studies where temperature tolerance is reduced under hypoxia. The limited use of oxygen as a variable may be attributed to our inability to measure in situ oxygen levels from satellites and the lack of precision in oxygen data from the ocean models (Robinson et al., 2011). Even so, recent advances in numerical ocean models have vastly enhanced information regarding ocean oxygen availability (Rose et al., 2017), but improved predictions at finer temporal and spatial scales are still needed (Breitburg et al., 2018). Even when accurate oxygen and temperature data are available, a complex non-linear interaction between the two variables (e.g. O 2crit ∼ temperature , Fig. 2a) on habitat suitability is often poorly represented by regression or envelope models (Golding and Purse, 2016). Advances in machine learning classification models can incorporate these non-linear complexities and overcome these difficulties (Elith et al., 2006).
We advocate the use of physiological models to transform environmental data into indices of habitat suitability and then to combine these with correlative species distribution modelling (SDM) approaches to bridge the gap between correlation and process and filter out interaction complexities (Austin, 2007;Kearney et al., 2008;Kearney and Porter, 2009;Dormann et al., 2012;Rodríguez et al., 2019). Physiological indicators, like aerobic scope, can project beyond the range of data and have been incorporated into predictive models for fisheries species (Peck et al. 2016;Marras et al. 2015). For example, Marras et al. (2015) used experimentally derived aerobic scopes of two fish species, the seabream (Sarpa salpa) and the marbled spinefoot (Siganus rivulatus), in the Mediterranean to develop a spatial thermal habitat suitability index and predict the distributions of both species. Cucco et al. (2012) developed an aerobic scopebased physiological model and combined it with temporally resolved oceanographic data to demonstrate that the migration of flathead mullet (Mugil cephalus) in the Mediterranean generally tracks the environmental conditions that optimize aerobic scope. Although these physiological-based mechanistic models overcome many of the limitations associated with purely correlative SDMs, they often assume that temperature is the only environmental variable that modulates how a fish population relates to its immediate environment (Sinclair et al., 2016).
By including oxygen availability together with temperature, φ overcomes the temperature-dependence limitations of aerobic scope. Whilst the aforementioned aerobic scope models have predictive accuracy, an increasing number of studies have found mismatches between temperatures where aerobic scope is optimized versus other important ecological/physiological processes such as development, growth or thermal preference, highlighting that aerobic scope is not always a reliable indicator of potential climate responses (Clark et al., 2013;Gräns et al., 2014;Norin et al., 2014). These aerobic scope-temperature mismatches usually occur at the warm edges of an organism's thermal range (Jutfelt et al., 2018). Our distribution modelling results support this notion as most of the absolute aerobic scope models over-predicted C. laticeps' distribution into warm tropical waters, far from where it occurs naturally ( Figure S1.5d in the Supporting Information). When aerobic scope models fail to explain ecological patterns, it may be because oxygen availability is not incorporated in a quantifiable way. Hypoxia has been shown to reduce the preferred temperature of teleost fish (Steffensen and Lomholt, 1991;Enders et al., 2019), but whether this is mediated through hypoxia effects on thermal optima for aerobic scope remains to be investigated. It is worth noting that where aerobic scope is found not be ecologically relevant, the study species sometimes occur in hypoxic waters like the Atlantic halibut (Hippoglossus hippoglossus) (Gräns et al., 2014) which occupies hypoxic areas in the north Atlantic Ocean (Breitburg et al., 2018) or barramundi (Lates calcarifer) (Norin et al., 2014), which are found in hypoxic mangroves in Australia, for example. However, this pattern is far from ubiquitous, and several studies based on normoxic species, such as the pink salmon (Oncorhynchus gor-buscha) (Clark et al., 2011) or coho salmon (Oncorhynchus kisutch) (Raby et al., 2016), have found absolute aerobic scope increases throughout temperatures that are not experienced by the species.
Distribution predictions for C. laticeps up to the end of the century indicate that the western and eastern edges will contract slightly due to the general warming trend along warm limits in the east and greater hypoxia in the west. Overall, the core distribution of C. laticeps is, however, predicted to persist up until at least 2100. This predicted contraction of C. laticeps' distribution is consistent with climate-mediated distribution shifts of fish species already recorded off the temperate South African coastal zone (Whitfield et al., 2016). The Cape anchovy (Engraulis encrasicolus) shifted its distribution from the southern Benguela to the Agulhas Bank during the mid-1990s, which was attributed to an increase in upwelling off the Agulhas Bank (Roy et al. 2007). Along the east coast, warming has resulted in a southward extension in distribution and abundance of tropical species (Mbande et al., 1995;James et al., 2008;Lloyd et al., 2012;Whitfield et al., 2016), which is consistent with the general poleward shift observed for fish species globally (Hastings et al., 2020). The extension of cool-temperate habitat and species in the west and subtropical habitat and species in the east is predicted to put a 'squeeze' on potential distribution shifts of warm-and cool-temperate species around South Africa (James et al., 2013;Potts et al., 2015;Whitfield et al., 2016) which is corroborated by the findings from this study.

12
predicted distribution contraction is a concern for any species, the most productive fishing grounds for C. laticeps in terms of catch per unit effort, west of Port Alfred and east of the Cape Town (Kerwath et al., 2013) are predicted will persist up to 2100 according to our data.
In conclusion, we highlight that the metabolic index (φ) has the potential to explain the whole distribution of marine fish, if calibrated across a wide temperature range. We also show that φ can predict the distribution of an endemic marine organism across a dynamic heterogenous marine environment. Importantly, whilst we highlight how different environmental variables (oxygen and temperature) limit different edges of the distribution for C. laticeps, this is through a common mechanism. Studies investigating the effects of climate change on marine organism distributions should therefore consider the combined, interacting effects of oxygen and temperature to be ecologically relevant (Townhill et al., 2017), and physiological indices such as φ offer a framework to achieve this.

Ethics
Ethics for this research was approved under the regulations of Rhodes University Animal Ethics (DIFS152025) and the Animal Use and Care Committee of the South African National Parks (004/16).

Supplementary material
Supplementary material is available at Conservation Physiology online.

Acknowledgements
Dr Sven Kerwath and The Department of Environment, Forestry and Fisheries are thanked for making C. laticeps catch data available from the National Marine Linefish System. Dr. Ekaterina Popova from the National Oceanographic Centre, Southampton, is thanked for providing access to the high-resolution ocean model used in this study. We thank the Aquatic Ecophysiology Research Platform of the South African Institute for Aquatic Biodiversity for use of laboratory space and equipment.

Funding
This research was funded by the Rhodes University Sandisa Imbewu Fund, the African Coelacanth Ecosystem Programme (Grant number 110762) and the SAIAB-NRF Aquatic Ecophysiology Research Platform. M.I.D. was funded by Rhodes University, the National Research Foundation and the Commonwealth scholarship commission. A.E.B. was supported by a Canada Research Chairs program.

Author contributions
MID, NCJ and WMP designed the study, MID collected the data and MID and AEB analysed the data. All authors reviewed and contributed to manuscript improvements.

Data availability
Occurrence data for Chrysoblephus laticeps is available from the National Marine Linefish System at The Department of Environment, Forestry and Fisheries in South Africa. Ocean model data is housed at the National Oceanographic Centre, Southampton. Experimental data is available in the supplementary material and analysis scripts are available from the lead author upon request.