Impact of deoxygenation and warming on global marine species in the 21st century

. Ocean temperature and dissolved oxygen shape marine habitats in interplay with species’ physiological characteristics. Therefore, the observed and projected warming and deoxygenation in the 21 st century of the world’s oceans in the 21 st century may strongly affect species’ habitats. Here, we implement an extended version of the Aerobic Growth Index 10 (AGI), which quantifies whether a viable population of a species can be sustained in a particular location. We assess the impact of projected deoxygenation and warming on the contemporary habitat of 47 representative marine species covering the epipelagic, mesopelagic/bathypelagic, and demersal realms. AGI is calculated for these species for the historical period and into the 21 st century using bias-corrected environmental data from six comprehensive Earth System Models. While habitat viability decreases nearly everywhere with global warming, the impact of this decrease is strongly species-dependent. Most 15 species lose less than 5% of their contemporary habitat volume over the 21 st century even at 23°C of global warming relative to preindustrial, although some individual species are projected to incur losses 2-3 times greater than that. We find that the in-habitat contemporary spatiotemporal variability of O 2 and temperature (and hence AGI) provides a quantifiable measure of a species’ vulnerability to change. In the event of potential large habitat losses (over 5%), species vulnerability is the most important indicator. Vulnerability is therefore more critical than changes in habitat viability, temperature or p O 2 levels. 20 Species’ vulnerability is the most important indicator for large (>5%) potential habitat losses - not relative or absolute changes in habitat viability (i


Introduction
Ocean temperature and dissolved oxygen (O2) are strongly linked by physical and biogeochemical processes as well as through their effects on the aerobic performance of aquatic ectotherms (Pörtner, 2010;Verberk et al., 2011;Breitburg et al., 2018;Oschlies et al., 2018;Seibel et al., 2021).Indeed, temperature and O2 are both found to be central in shaping species' distributions and are important climatic stressors to marine species worldwide (Perry et al., 2005;Doney et al., 2011;Poloczanska et al., 2013;Breitburg et al., 2018;Penn et al., 2018;Deutsch et al., 2020;Clarke et al., 2021).Observed and projected warming and deoxygenation are therefore likely to impact species.
Robust observational evidence of anthropogenically-forced deoxygenation is now becoming available as long-term O2 changes emerge from their natural variability (Frölicher et al., 2009;Long et al., 2016;Stramma et al., 2020;Buchanan and Tagliabue, 2021;Sharp et al., 2022).Specifically, an increase in the temporal and spatial resolution of observational data has allowed for the discovery and quantification of a ~2% decline in the global top-1000m O2 inventory since the 1960s (Ito et al., 2017;Oschlies et al., 2017;Schmidtko et al., 2017;Breitburg et al., 2018).This negative trend is projected to continue during the 21 st century for all climate scenarios (Bopp et al., 2013;IPCC, 2019;Kwiatkowski et al., 2020).More than 10% of deep ocean O2 is likely to be lost even if CO2 emissions were stopped in the year 2020 (Oschlies, 2021).Earth System Model simulations project an O2 loss between 100-600 meter depth of 13.27±5.28mmol m -3 for a high-emission low-mitigation SSP5-8.5 scenario and 6.36±2.92mmol m -3 loss for a lowhigh-emission highlow-mitigation SSP1-2.6 scenario by the end of the 21 st century (2080-2099 mean values relative to 1870-1899 ± the inter-model standard deviation; Kwiatkowski et al., 2020).
Ocean temperatures are changing as well: Ocean warming is a direct effect of atmospheric warming, as the ocean takes up approximately 90% of the excess heat from anthropogenic activities (Von Schuckmann et al., 2020).Global mean sea surface temperatures are observed to have increased by ~0.5 °C since the 1960s (Hersbach et al., 2020).Earth System Model simulations project sea surface warming of 3.5±0.8°Cfor a high-emission low-mitigation SSP5-8.5 scenario and 1.42±0.32°Cwarming for a low-emission high-mitigation SSP1-2.6 scenario by the end of the 21 st century (2080-2099 mean values relative to 1870-1899; Kwiatkowski et al., 2020).Most marine species will thus experience further warming of their habitat, considering that chances of limiting global atmospheric warming to 1.5°C are low even if all and unconditional pledges by countries are implemented in full and on time (IPCC, 2021;Meinshausen et al., 2022).Models that include more complex representations of species biology and ecology show that every tenth of a degree of global warming increases impacts on marine biodiversity, transforming species assemblages through changes in abundance, biomass and catch potential (Cheung et al., 2016).Moreover, the warming signal penetrates the deep ocean where it has major potential to affect marine ecosystems together with deoxygenation and ocean acidification (Levin and Le Bris, 2015).expected to decrease (Deutsch et al., 2015;Clarke et al., 2021;Gruber et al., 2021;Oschlies, 2021).Our approach newly includes depth variability as well as temporal variability of temperature and O2 in calculating AGI (Sect.2.1), which we apply to 47 representative species thanks to the generalized temperature dependence of pO2 demand in AGI (Sect.2.2).For environmental data of temperature and pO2 we use bias corrected Earth System Model projections from the latest version of the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Sect.2.3).Through this approach we assess the potential loss of viable contemporary habitat volume due to warming and deoxygenation for a representative selection of species -as well as identifying the drivers of such losses (Sect.3).

Aerobic Growth Index
We apply the Aerobic Growth Index (AGI; Clarke et al., 2021) to quantify species-specific impacts on habitat viability in response to changes in temperature and pO2.WeTherefore, we interpret AGI as a measure of habitat viability.AGI integrates growth theory, metabolic theory, and biogeography to calculate a theoretical ratio of pO2 supply over pO2 demand for each species i (Eq.1; rewritten from Eq. 14 in Clarke et al., 2021).
AGI can therefore be calculated for any species for which we have distribution data as well as environmental data for temperature and pO2.Temporal variations in pO2 and temperature are considered by using monthly climatological mean data from the World Ocean Atlas 2018 (WOA18 average of all available decades; Boyer et al., 2018;Garcia et al., 2019;Locarnini et al., 2019;Zweng et al., 2019).The horizontal distribution data are extended over the full depth range of each species species (0-200m for epipelagic species, 200-1000m for mesopelagic and bathypelagic species and bottom layer of the ocean data for the demersal species, thereby covering both deep and shallow demersal habitats; see Fig. C1 and Sect.2.3) to include the vertical variability of O2 and temperature in our estimate of pO2 threshold , T pref and hence AGI.The 0-200m depth range is used for epipelagic species, 200-1000m for mesopelagic species and the bottom layer of the ocean for the demersal species, thereby covering both deep and shallow demersal habitats (see Fig. C1 and Sect.2.3).This approach facilitates the estimation of species-specific, temperature-dependent critical pO2 levels and T pref despite the lack of observational data from multi-stressor laboratory experiments that apply to field conditions (Boyd et al., 2018;Collins et al., 2022).Different iterations of the metabolic index (Deutsch et al., 2015;Penn et al., 2018;Deutsch et al., 2020), which require laboratory-based estimates of temperature dependenttemperature-dependent critical pO2 levels, agree well with AGI in their assessment of habitat loss (Clarke et al., 2021) despite the much fewer data needed to calculate AGI.For additional details on the calculation of AGI we refer to Clarke et al. (2021).
An AGI of one implies that there is sufficient O2 supply for feeding, movement, and defence, but not growth and reproduction.
To sustain a viable population, additional aerobic scope is needed until AGI is above its critical value (AGI crit ) for a particular species.Following Clarke et al. (2021), AGI crit is the 10 th percentile of all AGI values in a species' habitat including vertical and temporal variability as done similarly for pO2 threshold and T pref .In this study a species is deemed impacted by changes in temperature and pO2 whenever AGI drops below AGI crit on an annual mean basis.All species information is listed in Table A1.The coarse resolution and the imperfect harmonization between the biogeographical, temperature and O2 data may affect the accuracy of the estimated AGI crit , as indicated by some species having AGI crit , at or below 1.We discuss in Sect. 4 how such biases may affect the results and conclusion, and how future studies can build on our results to improve the accuracy of the analysis.
Relative changes in AGI (AGI rel ) between time=t1 and time=t0 can be estimated from Eq. 2: wWhere ΔAGI is AGI(t1)-AGI(t0).Relative changes are thus entirely species-independent (in contrast to the metabolic index of Deutsch et al., 2015) and are interpreted as relative changes in habitat viability.Eqs. 1 and 2 show that j2-j1 (8000-4500=3500K) modulates the influence of the temperature effect on AGI.We maintain a reference period 1995-2014 throughout this study (i.e., AGI(t0) is the mean AGI over the years 1995-2014).Individual cContributions from pO2 and (temperature) to AGI rel are calculated by keeping temperature or (pO2) constant at its 1995-2014 mean state when calculating AGI rel .

Species data
AGI can be calculated for nearly a 1000 commercially exploited marine species due to the generalized temperature dependence of the pO2 demand .This broad applicability of AGI allows us to select 47 representative species (n=47), which are chosen such that depth level, climatic zone (tropical and temperate) and body size are represented.In addition, we selected some pelagic and deep-water wide-ranging species that inhabit both tropical and temperate regions, as well as the hypoxia toleranthypoxiatolerant species Dosidicus gigas.Three depth groups are represented through our selection: Epipelagic species (n=23) in the 0-200 meter depth range, mesopelagic/bathypelagic species (n=6) in the 200-1000 meter depth range and demersal species (n=18) at the sea floor (bottom wet layer of the models; Sect 2.3).The representativeness of our species selection is assessed in Sect.3.4.Species' pO2 threshold , T pref and AGI crit are listed in Table A1.The contemporary species distributions are based on a gridded product from Close et al. (2006) and are assumed to represent the 1995-2014 mean state (Fig. C1).
To account for mean errors and model drift, both a drift correction and bias correction were performed.First, the bottom layer ('seafloor') O2 ocean data were linearly detrended for piControl drift over the piControl years corresponding to the scenario years (1850 to 2100) as these trends are more than 10% of the scenario signal for some models.Drift in bottom layer temperature and salinity as well as upper water column O2, salinity and temperature were negligible as compared to the scenario signals and therefore not accounted for.Secondly, we performed a mean bias correction (e.g., Maraun, 2016) on all model data by subtracting the 'present-day' monthly mean climatological bias (the difference between WOA18 data vertically interpolated to the respective model levels and the respective model's 1995-2014 monthly climatology) from the entire simulated timeseries.We used the available WOA18 climatological mean product for the sea floor data because WOA18 climatologies are only available at monthly resolution until 1500m depth.The extracted original spatial resolution of the WOA18 data is 1° for O2 and 0.25° for temperature and salinity but these were all regridded to 1° to match the regridded model data grid.Note that to calculate the temperature biases, the model potential temperature was converted to in-situ temperature.
Finally, pO2 was calculated following Sect.E by Bittig et al. (2018), which is based on earlier work (Benson and Krause, 1984;Garcia and Gordon, 1992;Sarmiento and Gruber, 2006), and includes pressure correction (Taylor, 1978) and the correction for water vapor pressure (Weiss and Price, 1980) in the calculation of pO2 (Appendix B).
All results are presented at global warming levels (i.e., global mean air temperature at 2 m; e.g., Hausfather et al., 2022).In order to do so we first bias-correct modeled surface air temperatures such that the 1995-2014 global mean air temperature increase since 1850-1900 is 0.87 °C as observed (HadCRUT.5.0.1.0;Morice et al., 2021) in order to be consistent with the bias correctedbias-corrected ocean temperature and oxygen data.We then find the first year where the 15-year running mean of these bias-corrected global mean surface air temperature data is greater than or equal to the warming level of interest and calculate the 15-year running mean at that year over the data for the analyses.For warming levels above 1.5ºC, we only use the results for SSP5-8.5 as not all models reach warming of more than 1.5ºC in SSP1-2.6.

Global changes in warming, deoxygenation and habitat viability
Ocean temperature is projected to increase (Fig. 1a) and pO2 to decrease (Fig. 1b) with global warming of the atmosphere.
These changes occur in all three depth layers considered here, and for all CMIP6 models.Noticeably, the response of ocean warming and deoxygenation to global atmospheric warming is approximately linear (Fig. 1).From a linear fit to the multimodel mean CMIP6 changes in Fig. 1 we find that the epipelagic realm warms the most with 0.50±0.03°C per degree of atmospheric warming (standard deviation given across the individual model fits).This signal is dampened toward depth to 0.18±0.02°C per degree of global warming in the mesopelagic/bathypelagic realm and 0.08±0.01°C per degree of global warming in the demersal realm (Fig. 1a).In addition to the warming, we find that the epipelagic realm loses 0.40±0.55mbar of pO2, the mesopelagic/bathypelagic loses 1.35±0.89mbar of pO2 and the demersal realm loses 1.17±0.97mbar of pO2 on average and per degree global warming of the atmosphere (Fig. 1b).hence Thaving the largest changes in pO2 realized are projected at depth in contrast to ocean warming.The warming (Fig. 1a) and deoxygenation (Fig. 1b) drive reduce AGI relative to its contemporary state (i.e., a negativea reduction in AGI rel ; Fig. 1c)), which we interpret as a loss of habitat viability (Sect. 2.1; Fig. 1c) that is independent of species (Eq.2).In the epipelagic, AGI rel is decreasesreduced by 2.17±0.69% per degree of global warming (Fig. 1c), while AGI decreases 2.33±1.64% per degree of global warming .The global mean reduction in habitat viability (i.e., AGI rel ) in the mesopelagic/bathypelagic realm is 2.33±1.64% per degree of global warming.Last, The demersal decrease in AGI loss of habitat viability is 0.86±0.48% per degree of global warming, making it is the least pronounced of the three studied depth intervals.The approximately linear response of marine warming, deoxygenation, and loss of habitat viability to global atmospheric surface warming (Fig. 1) highlights and confirms that all any additionally realized atmospheric warming will affect the marine environment (Cheung et al., 2016).

Local changes and drivers of habitat viability
A relative reduction in habitat viability (i.e., a reduction negativein AGI rel ; Fig. 1c) is projected to occur almost everywhere at 2 ºC of global warming (Fig. 2a-c; see Fig. C3 and C4 for 1.5 and 3 ºC of global warming, respectively), indicating that for most habitats and therefore species we expect a reduction in habitat viability.The relative reduction in habitat viability is generally largestr in the epipelagic and mesopelagic/bathypelagic realms (Figs. 1c,2a,b), but the larger spatial heterogeneity at mesopelagic/bathypelagic depths reveals that locally mesopelagic/bathypelagic AGI rel reductions may far exceed those in the epipelagic, particularly in the North Pacific (Fig. 2b).Hence the location of a species' habitat, both vertically and horizontally, is key to projected changes in habitat viability for a specific species.Note that the patterns in each of the panels of Fig. 2 remain similar for higher degrees of global warming, only the intensity of change increases (not shown), which agrees with the approximately linear response of the global average AGI rel to global warming (Fig. 1c).When considering the contribution from the two drivers of AGI change, pO2 and temperature changes, AGI rel at 2 ºC of global warming is driven mostly by temperature in the epipelagic and by pO2 in the mesopelagic/bathypelagic and demersal realms (Fig. 2d-i).The AGI rel at 2 ºC of global warming due to pO2 is -0.16±5.12%for the epipelagic, -2.52±6.96%for the mesopelagic/bathypelagic, and -0.62±2.02%for the demersal realm (Fig. 2d-f), while the respective AGI rel due to temperature are -2.32±1.36%,-0.91±1.18%,and -0.39±0.95%(Fig. 2g-i).In the mesopelagic/bathypelagic the drivers of loss in habitat viability depend more strongly on location (Fig. 2e,h).Hence, ~Globally, on average 9487% of AGI rel is on average driven by the relatively pronounced warming in the epipelagic (black lines in Fig. 1c) since changes in pO2 are minor (Fig. 1b).This is because because the epipelagic realm is generally well ventilatedwell-ventilated with O2-rich surface waters.For the mesopelagic/bathypelagic (demersal), warming accounts for only 27% (39%) of the total AGI rel is for only 274% (397%) explained by warming (black lines in Fig. 1c).In the mesopelagic, the drivers of loss in habitat viability depend more strongly on location (Fig. 2e,h).Nevertheless, theThe larger contribution from pO2 to AGI rel increases uncertainty for the mesopelagic/bathypelagic and demersal realms because model projections are uncertain for pO2 (Fig. 1b, C52).In some regions the effects of pO2 and temperature on AGI rel may compensate each other and result in negligible changes in AGI.We find examples of this in the Northern Indian Ocean at epipelagic depths, in the Gulf of Guinea at mesopelagic depths, and in the North Atlantic around Iceland at demersal depths.
AGI rel is has large model uncertaintyuncertain for species having a large part of their habitat in eastern-boundary equatorial upwelling regions or around Antarctica at epipelagic depths, the western equatorial Pacific at mesopelagic/bathypelagic depths, north of the equator in the Indian Ocean at epipelagic and mesopelagic/bathypelagic depths or regions scattered across all ocean basins for demersal depths (hatched areas in Fig. 2a-c).Most of this uncertainty is coming from pO2 (Fig. which agrees with the approximately linear response of AGI rel to global warming (Fig. 1c).
Besides considering the model uncertainty, we performed a sensitivity analysis of AGI rel to the choice of generalized temperature dependence parameters (i.e., j2-j1).If j2-j1 is adjusted to represent an arbitrary low temperature sensitivity of 1000K, global mean AGI rel is 34% of the standard case j2-j1=3500K in the epipelagic, 67% in the mesopelagic and 73% in the demersal realm.On the other hand, for an arbitrary high temperature sensitivity (j2-j1=6000K), global mean AGI rel is 165% of the standard case j2-j1=3500K in the epipelagic, 118% in the mesopelagic and 126% in the demersal realm.Projections for epipelagic species are therefore most sensitive to the choice of j2-j1, as temperature changes are largest there.Further work is needed to explore the uncertainty in j2 and j1.

Impacts of AGI changes on habitat volume of individual species
The overall decrease negativein AGI rel and hence relative loss of habitat viability with global warming (Figs.1c and 2) causes loss of contemporary habitat volume (i.e., newly exposed volume with AGI<AGI crit ) and hence local extinction for species at each of the studied depth ranges (Figs. 3 and 4).Habitat loss is expressed relative to its contemporary volume (Fig. 3) to facilitate comparison between wide-ranging and more narrowly distributed species.Loss of contemporary habitat is generally less than ~150% at 21.5°C global warming, and mostly under 5%, but increases to up to ~25% for individual species at 3°C (Fig. 3).Wide-ranging epipelagic species (e.g., Acanthocybium solandri, Coryphaena hippurus, Katsuwonus pelamis, Thunnus obesus, or Elagatis bipinnulata;: Fig. C1) experience losses of contemporary habitat volume of less than 5% for any of the analyzed warming levels, while more narrowly distributed species experience the largest losses of up to ~25% of their contemporary habitat at 3°C global warming (e.g., Micromesistius poutassou, Thunnus atlanticus, Sebastes mentella, or Anarhichas denticulatus).Notably, species that have the largest contemporary habitat loss at 1.5°C generally are those species that also lose the most at 3°C of global warming, which is in line with the earlier findings of approximately linear response of relative changes in habitat viability to warming and deoxygenation (Fig. 1).Any early (i.e., 1.5°C) response of a species to warming and deoxygenation is therefore a warning indicator for additional loss of contemporary habitat at increased levels of global warming.
We separately assess the impact of the uncertainty of AGI crit on these results by calculating habitat loss with an AGI crit of a) minimum AGI, b) 5 th percentile, c) 10 th percentile (i.e., the default case), d) 15 th percentile and e) 20 th percentile of in-habitat AGI.We find that even when including much higher thresholds (AGI crit as 20 th percentile), our results are similar with a few species having large losses but most losing less than 5% at 2°C of warming relative to the 1995-2014 state (Fig. C63).
Moreover, a sensitivity analysis for species Thunnus atlanticus and Gadus morhua shows that our median result is robust to the choice of the generalized temperature dependence parameters j2-j1 (we explored j2-j1 ± ~710%; Fig. C7).
Absolute losses in habitat volume (i.e., loss expressed in volumetric terms instead of a percentage) show that small relative losses (Fig. 3) often correspond to the largest volumetric losses (Fig. 4).As an example, median Thunnus alalunga habitat loss is less than ~2.5% at any of the analyzed warming levels (Fig. 3), while absolute losses are the largest of all 47 species at ranging from ~0.25 to -1.5e 6 km 3 -depending on the global warming level (Fig. 4).On the other hand, we find species like Sebastes mentella for which relative losses are large (median ~8-26% of the contemporary habitat depending on global warming level; Fig. 3) while absolute losses are comparably small (~0.6-1.8e 5 km 3 ) because the contemporary volume of Sebastes mentella is relatively small (Fig. 4).Note that epipelagic species lose habitat volume in the order of a million km 3 .: In comparison, the entire Black Sea has a volume of about 0.5 million km 3 .Depending on the location of viable contemporary habitat loss, for species of commercial interest such large absolute loss can be particularly impactful to local fisheries.For most species, temperature is the main driver of habitat loss (black stars in Fig. 3).Exceptions exist for example in the mesopelagic/bathypelagic, where pO2 drives about half of the habitat loss for the two species with the largest loss (i.e., Sebastes mentella and Aphanopus carbo) as well as for the demersal species Anarhichas denticulatus.Even though most of the realized loss can be explained by warming, not all species have large losses despite warming being relatively uniform although dampened toward depth (Kwiatkowski et al., 2020).We understand Tthese differences can be explained from by considering the original spatial and temporal pO2 and temperature variability in each species' habitat, which shapes their vulnerability to change.This is investigated next.(Sect.3.4).

Drivers of habitat volume loss of individual species
The differences in habitat loss between species as shown in Figs.We can quantify theis "vulnerability" of a species to changes in AGI by calculating the cumulative sum of the PDFs (i.e., the cumulative density function, CDF;, conceptual Figs.5c,d and species-specific results in C95).The slope of the CDF at a cumulative density of 0.1 (i.e., 10% of the volume where AGI crit is defined) indicates the potential loss in habitat for a certain change in AGI (Fig. 5 and C95).If the slope of the CDF is steep at the critical threshold, the species is relatively vulnerable to warming and deoxygenation as: oOnly a small reduction in habitat viability (i.e., AGI) will push a relatively large volume below the critical threshold.AIn example is given in Fig. 5, where for an identical change in mean in-habitat DAGI of 0.1 just 1% of the volume is pushed below AGI crit for a species with a small slope of 0.14 (Fig. 5b,d, 'Thunnus obesus' schematic), while the same change in AGI results in 9% volume loss for a species with a large slope of 1.67 (Fig. 5a,c, 'Thunnus atlanticus' schematic).Changes in the slope of some species' CDFs indicates that different vulnerabilities exist for different parts of that species' habitat (Fig. C9).Hence, in habitat areas which that are represented by a part of the CDF with a relatively steep slope, a relatively small change in AGI is needed to bring a relatively large volume closer to AGI crit .Nevertheless, onlyOnly the CDF slope at AGI crit relates directly to viable habitat volume loss as only AGI values below AGI crit are considered to have an impact on habitat volume.
Indeed, projected habitat volume loss increases with species' vulnerability (i.e., the CDF slope at AGI crit ; Fig. 6d), as well as to a lesser extent with warming and deoxygenation (Fig. 6b,c).Notably, the largest absolute reductions of mean in-habitat AGI do not indicate those species who lose most contemporary habitat volume (Fig. 6a).On the contrary, the environmental state of the contemporary habitat as captured in the PDFs and thus the slope of the CDFs and vulnerability is the strongest indicator for impact: 87% of the variance in volume loss at 2 °C global warming can be explained by vulnerability (R 2 of linear fit in  Habitat viability thus strongly depends on the variability of temperature and O2 in the habitat of the species as captured by the species' vulnerability (Figs.C95 and 6d).Therefore, reports of relative losses in habitat viability based on pO2 supply over pO2 demand ratios (e.g., Deutsch et al., 2015;Oschlies, 2021) should not be interpreted as leading to actual reductions of viable habitat for individual species as they do not include species-specific thresholds nor their vulnerabilities.
Considering the range captured in Fig. 6a-d we expect that our selection of species is representative of a wide range of marine ectotherms.
Interestingly, species with high vulnerability and loss (red markers in Fig. 6) all have a high pO2 threshold above ~150 mbar (Table A1, Fig. C126).Thus, even though warming explains most of the loss of contemporary habitat (black stars in Fig. 3), loss is only high for vulnerable species (Fig. 6d) -which in turn all are sensitive to pO2 as evidenced by their pO2 threshold above ~150 mbar (Fig. C126).A high sensitivity to pO2 and hence a high pO2 threshold is therefore an indicator of vulnerability, although not all species with high pO2 threshold are vulnerable (Fig. C126).The high vulnerability for species with a high pO2 threshold shows that also species in well-oxygenated regions can be vulnerable to climate change as their natural pO2 range is limited.We further note that vulnerability does not depend on the depth realm of a species.Resilient species (blue markers in Fig. 6) have strong spatiotemporal variability of AGI (broad PDF in Fig. C84) such that even large mean changes of AGI (Fig. 6a) do not expose a large volume to subcritical AGI values.Noticeable is that the two species with a PDF skewed to the right (Fig. C8; Hippoglossus hippoglossus and Theragra chalcogramma) are both in this group, while all other species tend to have a leftskewed PDF of AGI values in their habitat.These two species are both demersal-dwelling and are very pO2 tolerant (i.e., low pO2 threshold ; Table A1) and have a wide range of different AGI values in their habitat, with a relatively large volume of high-AGI values causing the right skew (Fig. C8) and resilience (Fig. C12).Whether AGI is the right metric for determining habitat viability for these two species needs further investigation that goes beyond the scope of this study.The six species with relatively high vulnerability but small habitat volume losses (yellow markers in Fig. 6) experience relatively small AGI changes in their habitats even at 3°C global warming (Fig. C84) thereby preventing large habitat losses.

Discussion
We introduce a new version of AGI that adds vertical temporal variability in the calculation of pO2 threshold , T pref , and AGI crit ,which makes it possible to assess volumetric habitat changes.The original AGI applies and assesses temporal variability in the horizontal (surface or bottom ocean layers for pelagic and demersal species, respectively) direction only (surface or bottom ocean layers for pelagic and demersal species, respectively; Clarke et al., 2021), as commonly practiced (e.g., Bryndum-Buchholz et al., 2019;Tittensor et al., 2021).In other words, e, i.e., either surface or sea floor data were applied for the calculation of pO2 threshold , T pref , AGI crit and hence AGI in earlier work.To assess the differences between our new approach and the original approach we repeated the analysis as presented in Fig. 3, (Fig. C137), now using surface ocean data only for mesopelagic/bathypelagic and epipelagic species as well as calculating pO2 threshold , T pref , AGI crit from the surface monthly mean WOA18 data only (Fig. C13).We find that the sensitivity to global warming of all species is higher for the original AGI as compared to our new approach which includes vertical and seasonal variability of temperature and pO2.This is understood from the combination of a) limited spatial variability of surface ocean pO2 as well as temperature, leading to higher T pref and pO2 threshold estimates and therefore stronger sensitivity to warming and deoxygenation as compared to our new approach and b) larger AGI rel changes closer to the surface.We expect that including temporal and vertical spatial variability in calculating AGI provides a more realistic estimate of the pO2 and temperature variability experienced by a species and therefore a better estimate of its sensitivity to warming and deoxygenation.Nevertheless, we acknowledge that further increasing spatiotemporal resolution (e.g., using daily-mean data and including interannual variability) can may increase affect variability (Deser et al., 2009;Baumann et al., 2015) which can affect estimates of T pref , pO2 threshold and AGI crit .Unfortunately, no established theory exists yet to decide what temporal variability in environmental parameters best captures species' T pref , pO2 threshold or AGI crit .By considering WOA18 monthly mean climatological data as the basis for our estimates of T pref , pO2 threshold and AGI crit , we are consistent with the time resolution of the CMIP6 model data (monthly mean).
RFurther, regarding the choice of the 10 th percentile threshold and impact of its uncertainty on our results (Fig. C63), we consider an AGI crit threshold above the 20 th percentile of in-habitat AGI values unlikely as then by definition already 20 percent of the habitat would be unsuitable to sustain a viable population of that species.Nevertheless, for species where AGI crit is very close to 1 or even below 1 (Table A1), a higher percentile may be warranted as a meaningful critical value.At the 10 th percentile, some uncertainty in the species-specific physiological parameters is considered.We find for most species that the 10 th percentile is located at an AGI above which habitat volume steeply increases suggesting it acts as an appropriate threshold (Fig. C84).
Regarding species' data, we assume that our results can be generalized to commercial fish and invertebrates worldwide, as they are based on representative species from different climatic zones (tropical, temperate, polar), vertical habitat (epipelagic, mesopelagic, bathypelagic, demersal), geographic range breadths, taxonomic groups (fish and invertebrates) and size classes.
Species distribution ranges were generated by an algorithm developed by the Sea Around Us project (see Close et al., 2006;Cheung et al., 2008).The resulting distributions, and the parameters used for their construction are available at http://www.seaaroundus.org.These distributions have been used to project climate-impacts on fishery resources in a great number of studies (Cheung et al., 2009;Cheung et al., 2010;Fernandes et al., 2013), and are assumed to represent species distributions in over the period 1995-2014 (Tai et al., 2021).Our assumption to extend the 2D distributions provided by Close et al. (2006) over the entire depth range of each species' depth realm is driven by data sparsity and reliability of 3D species distributions for our selection of species.When reliable 3D habitats, or even time-varying habitats, can be constructed from species' observations these could be included (e.g., distribution data are continuously collected in the Ocean Biodiversity Information System but are currently are too sparse to provide 3D distribution data).Some species may be limited to only part of their assigned depth range or live partly (and possibly temporarily) above or below it.Nevertheless, we expect that the assigned depth range generally provides a good estimate of in-habitat pO2 and temperature variability, which affects pO2 threshold , T pref and therefore AGI and AGI crit .
Our results for the mesopelagic include two vertical migrators (Dosidicus gigas and Aphanopus carbo).As opposed to most other species, the distribution range of vertical migrators is limited at the cold boundary of the distribution because of their low aerobic scope in cold waters (Seibel and Birk, 2022).Therefore, the temperature sensitivity of these species is likely not captured by the generalized temperature dependence in AGI, and contemporary habitat loss due to warming and deoxygenation as estimated for Aphanopus carbo is likely overestimated.We nevertheless project negligible loss of contemporary habitat for Dosidicus gigas (Fig. 3) due to its low vulnerability and low pO2 threshold , which is in good agreement with the findings of Seibel and Birk (2022) despite the generalized temperature dependence of AGI.In addition,S species-specific thresholds pO2 threshold and AGI crit and preference T pref are calculated based on the in-habitat spatiotemporal variability of pO2, temperature and AGI respectively.This is done in lack of observation-based thresholds and preferences that translate to field conditions (Boyd et al., 2018;Collins et al., 2022).Detrimental effects from deoxygenation such as reduced visionas visual hypoxia actually become relevant at much higher pO2 than (near) lethal pO2 levels (Mccormick and Levin, 2017), while only the latter is often what is assessed in the lab.As an effect, the exact threshold of impact remains unknown and probably depends on many factors including the impact itself, and the abruptness, magnitude, intensity, duration, heterogeneity, and recurrence of exposure to subcritical values (Gruber et al., 2021), as well as timing of and adaptability to unfavorable temperatures, subcritical pO2 and hence subcritical AGI.
Through the bias correction of the CMIP6 model data all monthly mean biases relative to WOA18 are removed from our analysis.We acknowledge the influence of observational uncertainties as well as resolution mismatch between our model and the observational WOA18 data used in our bias correction (Casanueva et al., 2020).More complex bias adjustment such as correction for variance biases is prevented by the spatial and temporal resolution of the model and observation data at the global scale.The ongoing effort to collect, compile, and quality-control O2 data in open-access repositories (e.g., Grégoire et al., 2021) will hopefully make it possible to do more advanced bias correction in the future.Until that time the strong temporal variability and spatial heterogeneity of O2 trends complicate the comparison between model and observational data.
Nevertheless, the remaining forced response of the models likely underestimates deoxygenation (Andrews et al., 2013;Oschlies et al., 2017;Oschlies et al., 2018;Buchanan and Tagliabue, 2021) and overestimates atmospheric warming (Tokarska et al., 2020) and therefore ocean warming for some CMIP6 models.Part of these warming biases are due to the relatively high climate sensitivities in the CMIP6 models (Meehl et al., 2020).As a further measure to limit model uncertainty, we therefore present results at different global warming levels such that they are insensitive to the differences in model climate sensitivity (Hausfather et al., 2022).We last acknowledge the relatively coarse resolution of the CMIP6 data (typically ca.100km in the ocean) which for species with a highly local distribution (Fig. C1) may lead to higher model uncertainties, especially along the coasts where model disagreement is larger (Fig. C52).
Our approach may give a conservative estimate of contemporary habitat loss since a) crossings of the critical thresholds on timescales shorter than a year are excluded from our analysis, b) CMIP6 projections likely underestimate deoxygenation (Andrews et al., 2013;Oschlies et al., 2017;Oschlies et al., 2018;Buchanan and Tagliabue, 2021), but considering the importance of temperature in driving habitat loss (Fig. 3), especially in the epipelagic realm, the uncertainty of pO2 projections likely has a relatively small effect on our results, and c) we do not include other potential stressors on species' habitats in our analysis such as acidification, changes in ecosystem structure, overfishing, marine phenology, disease pressure, food resources, predation pressure, pollution or eutrophication (e.g.; Poloczanska et al., 2016;Bindoff et al., 2019;Whalen et al., 2020).
Examples of crossings of the critical thresholds on timescales shorter than a year would be short hypoxic events and marine heatwaves (Frölicher and Laufkötter, 2018;Jacox et al., 2020;Cheung et al., 2021).Projected deoxygenation and particularly hypoxic or anoxic events have the potential to worsen and even surpass the effects of warming, marine heatwaves, and acidification (Gruber et al., 2021;Sampaio et al., 2021).
On the other hand, for some species the impact will be overestimated if they are able to adapt to future warming and deoxygenation (Cheung et al., 2009;Pinsky et al., 2013;García-Molinos et al., 2016;Palumbi et al., 2019;Collins et al., 2021;Liao et al., 2021).Further note that we considered potential loss of contemporary habitat only: Mobile species have been observed to redistribute based on the rate and direction of climate change (Pinsky et al., 2013) which can preserve the species range area if they are able to expand into newly suitable areas -however this can alter the original ecosystem structure and function.
For most species we find a loss of habitat volume of less than 10%.It is found for example by Gotelli et al. (2021) that only a small percentage of species drives the observed changes in marine species assemblages, showing that even when only a few species experience large losses, impacts can be profound for the ecosystem as a whole.For the individual species however, the loss of only a small fraction of their contemporary habitat likely provides adaptation opportunities.Our results imply that species that are deemed vulnerable due to their limited range of in-habitat pO2 and temperature are likely to be the most impacted by global warming (i.e., 'vulnerable species' in Fig. 6 and species with steep CDF slopes in Fig. C9).Our study can therefore inform e.g., fisheries management by identifying species particularly vulnerable to ocean warming and deoxygenation.Such identification provides species-specific information complementing earlier studies that found reduced impact on fisheries at lower levels of global warming (Cheung et al., 2016).Indeed, for every tenth of a degree ofany additional global warming, our study shows increased marine deoxygenation and warming as well as increased loss of contemporary habitat across all species albeit with a strongly species-specific magnitude.These results confirm the need to limit global warming levels to the minimum to prevent loss of contemporary habitat and support the identification of the species that would be most vulnerable to marine deoxygenation and warming.

Conclusions
• Marine warming and deoxygenation are projected to intensify with global warming and drive a relative decrease in global habitat viability penetrating to all depths (Fig. 1 and 2).
• The generally negative relative decreasechanges in habitat viability (i.e., AGI rel ) areis dominated by warming at the surface while deoxygenation becomes increasingly important with depth (Fig. 2).
• Species' Lloss of contemporary habitat is driven mostly by warming in the epipelagic realm, while in the mesopelagic/bathypelagic and demersal realms reduced pO2 is also contributing for some species (Fig. 3).
• Deoxygenation and warming cause most species to lose less than 5% of their contemporary habitat volume over the 21 st century relative to preindustrial (Fig. 3).Some individual species are however projected to incur losses 2-3 times greater than that at 1.5 and 2 ºC of global warming and 4-5 times greater at 3ºC of global warming.At 2 ºC of global warming, epipelagic losses are generally in the order of 0.1-0.5 million km 3 , while mesopelagic losses are 0.01-0.15 million km 3 and demersal losses are in the order of about 0.00025 million km 3 .
• Species' vulnerability is shown to be the most important indicator for potential large (>5%) potential habitat lossesnot relative or absolute changes in AGI, pO2 or temperature (Fig. 6).A species' pO2 threshold above ~150 mbar is an indicator for high species vulnerability to warming (Fig. C126).Our approach of quantifying vulnerability can help identify those species most vulnerable to marine warming and deoxygenation.
Where The term exp ( ) (unitless) is the pressure correction term for O2 sat .We then first calculate the saturation concentration of O2 in seawater (Garcia and Gordon, 1992) in mol kg -1 in using Eq.B2.

Figure 1
Figure 1 Simulated gGlobal mean changes in ocean insituin-situ temperature in ºC (a), pO2 in mbar (b) and AGI rel in % (c) for different depth layers and global warming levels, where global warming is calculated as global surface air temperature increase relative to the 1850-1900 mean.The multi-model running mean is given in opaque blue (SSP1-2.6)and red (SSP5-8.5)and has for several decades the corresponding 20-year multi-model mean year labeledlabelled.Individual models are shown in transparent light blue and red without taking a running mean.AGI rel is given relative to the 1995-2014 mean as in the remainder of the manuscript, and the mean contribution from temperature only (excluding the small effect of temperature on pO2) is indicated by the black line in panel (c) and calculated by keeping pO2 constant at its 1995-2014 mean state when calculating AGI rel .AGI rel (%) is entirely species-independent (Eq.2) and values that exceed 1000% or are below -1000% were excluded during the calculation of the global mean AGI changes to omit several grid-cells with extreme outliers caused by very small absolute changes in O2 causing very large changes in AGI rel .

Figure 2
Figure 2 Multi-model mean AGI rel relative to the 1995-2014 mean at 2 °C global warming (using SSP5-8.5 simulations), vertically averaged over the epipelagic and mesopelagic realms, and shown for the demersal realm (a-c).AGI rel is split up into the contributiondivided into contributions from (d-f) pO2 and (g-i) temperature.Data are hatched where more than 2 out of the 6 models disagree about the sign of change.Note that sea floor depth and thus demersal depth dependsdepend strongly on location.Contributions from pO2 (temperature) are calculated by keeping temperature (pO2) constant at its 1995-2014 mean state when calculating AGI rel .Further note that since [O2] depends on temperature too, the contribution to AGI rel from pO2 also contains a minor temperature component.
2 d-f, Fig.C52), with temperature contributing to uncertainty in the North-Atlantic south of Greenland and in the western equatorial Pacific at mesopelagic/bathypelagic depths.Exceptions to the decrease in AGI rel are limited to some smaller parts of the world's oceans including equatorial regions and the North-Atlantic south of Greenland in the epi-and mesopelagic/bathypelagic, and around the Antarctic continent in the epipelagic.Model disagreement is generally large in these regions of positive AGI rel increase and is mostly attributable to projected increases in pO2 which have large uncertainties (hatching in Fig.2a-f and model range in Fig.C52).Notably, in some regions the effects of pO2 and temperature on AGI rel may compensate each other and result in negligible changes in AGI: We find examples of this in the Northern Indian Ocean at epipelagic depths, in the Gulf of Guinea at mesopelagic/bathypelagic depths, and in the North Atlantic around Iceland at demersal depths.Note that the patterns in each of the panels of Fig.2remain similar for higher degrees of global warming, only the intensity of change increases (not shown),

Figure 3 Figure 4
Figure 3 Percentage of remaining Habitat change (%) of contemporary (1995-2014) habitat volume for different levels of global warming, with negative values indicating habitat loss and positive values indicating habitat gain.Note the different y-axis scale for 3 or 4 are better understood from the probability density of contemporary(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) in-habitat AGI for each species (conceptual Fig.5and species results in Fig.C84).The spatial variability of the contemporary pO2 and temperature in each species' habitat results in a species-specific probability density function (PDF) for AGI (black lines in Fig.5a,b).Depending on this shape, a given reduction in AGI (ΔAGI) exposes a relatively large or small part of the species' habitat to subcritical AGI values (red lines and stippling in Fig.5a,b), thereby causing volume loss.We can quantify this "vulnerability" of a species to changes in AGI by calculating the cumulative sum of the PDFs (i.e., the cumulative density function, CDF, Figs.5c,d and C5).

Figure 5
Figure 5 Conceptual figure based on Thunnus atlanticus (a,c) and Thunnus obesus (b,d) showing the difference in impact (change in volume ΔV) of an example mean AGI reduction of ΔAGI of 0.1 (i.e., habitat-mean ΔAGI=0.1)below the 1995-2014 contemporary mean (black lines).This difference is shown to be related to the shape of the PDF and the slope of the CDF at 0.1 (i.e., at AGI crit ), which we refer to as the species' "vulnerability".

Fig. 6d )
Fig. 6d).This result holds across different levels of global warming.:At 1.5°C of global warming, 85% of the variance in volume loss can be explained by vulnerability, and at 3°C of global warming this is 88% (see Figs. C10 and C11).The correspondent linear equation taken across all depth realms is   (%) = 7.31 * vulnerability − 0.10.

Figure 6
Figure 6 Multi-model mean in-habitat changes at 2°C of global warming of (a) AGI, (b) temperature, (c) pO2 (in SSP5-8.5) and (d) vulnerability (CDF slope at a cumulative density of 0.1 based on 1995-2014 mean data, Fig. C9) plotted against loss of contemporary habitat volume for each species (model range indicated by error bars).Species with > 5% loss marked in red, more than -0.1 ΔAGI in blue, and volume loss < 5% as well as vulnerability > 0.3 in yellow.There is no uncertainty in the vulnerability calculation because all models have the same 1995-2014 CDF slope due to the WOA18 bias-correction.From a linear regression to the data which is plotted in dashed grey we find an R 2 of 0.0% for (a) which line is therefore not plotted, 18% for (b), 21% for (c) and 87% for (d).

Figure
Figure C2 Global mean changes in ocean in-situ temperature in ºC (a), pO2 in mbar (b) and AGI rel in % (c) for years 1850-2100.The 575

Figure
Figure C4 Same as Figure C3, but for 3°C global warming (and therefore using SSP5-8.5 simulations only).

Figure
Figure C52 Multi-model range of AGI rel at 2°C global warming for the three depth intervals studied.

Figure
Figure C63 Habitat change (%) Percentage of remaining contemporary (1995-2014) habitat volume for 2°C global warming shown similar to Fig. 3, including 5 levels of AGI crit in every species' boxplot (number of datapoints n=5 AGI crit levels * 6 models = 30): AGI crit is taken as the minimum in-habitat AGI value, the 5 th percentile, the 10 th percentile, the 15 th percentile and the 20 th percentile, respectively.Note the different y-axis when comparing to Fig. 3.Each boxplot indicates the median in orange and a box bounded by

Figure
Figure C84 Probability density of AGI for each species for the contemporary reference period 1995-2014 (in black) and for 3°C 615

Figure
Figure C95 Cumulative density function of AGI at 1995-2014 for each species with "vulnerability" annotated in upper right corner 620

Figure
Figure C10 Multi-model mean in-habitat changes at 1.5°C of global warming of (a) AGI, (b) temperature, (c) pO2 (mean over SSP1-2.6 and SSP5-8.5)and (d) vulnerability (CDF slope at a cumulative density of 0.1 based on 1995-2014 mean data, Fig. C9) plotted against loss of contemporary habitat volume for each species (model range indicated by error bars).There is no uncertainty in the vulnerability calculation because all models have the same 1995-2014 CDF slope due to the WOA18 bias-correction.From a linear regression to the data which is plotted in dashed grey we find an R 2 of 0.0% for (a) which line is therefore not plotted, 25% for (b), 27% for (c) and 85% for (d).

Figure
Figure C11 Same as Fig. C10, but for 3°C global warming (and therefore using SSP5-8.5 simulations only).

Figure C12 For
Figure C12For each species T pref , AGI crit and pO2 threshold (see also TableA1) plotted against their vulnerability (=slope at cumulative density of 0.1; i.e., at AGI=AGI crit ; see also Fig.C9, Fig.6, TableA1).Colors as in Fig.6:For 2 °C global warming, species with > 5% loss marked in red, more than -0.1 ΔAGI in blue, and volume loss < 5% as well as vulnerability > 0.3 in yellow.

Figure
Figure C137 Habitat change (%) of contemporary (1995-2014) habitat volume for different levels of global warming, with negative values indicating habitat loss and positive values indicating habitat gain.Percentage of remaining contemporary (1995-2014) habitat volume of the surface layer (or bottom layer/sea floor for demersal species) for different levels of global warming.This figure is the same as Fig. 3, but here applying surface values only in calculating pO2 threshold , T pref and AGI crit and hence AGI and volume<AGI crit .