Long-term impacts of global temperature stabilization and overshoot on exploited marine species

Global warming alters ocean conditions, which can have dramatic consequences for marine species. Yet, the centennial-scale effects and reversibility of habitat viability for marine species, particularly those that are important to fisheries, remain uncertain. Using the Aerobic Growth Index, we quantify the impacts of warming and deoxygenation on the contemporary habitat volume of 46 exploited marine species in novel temperature stabilization and overshoot simulations until 2500. We demonstrate that only around half of the simulated loss of contemporary (1995–2014) habitat volume is realized when warming levels are first reached. Moreover, in an overshoot scenario peaking at 2 °C global warming before stabilizing at 1.5 °C, the maximum decrease in contemporary habitat volume occurs more than 150 years post-peak warming. Species' adaptation may strongly mitigate impacts depending on adaptation rate and pressure. According to our study, marine species will be affected for centuries after temperature stabilization and overshoot, with impacts surpassing those during the transient warming phase.

1 Introduction

Marine fisheries play a pivotal role in the global food system, generating economic benefits (Sumaila et al., 2012), supplying essential nutrients to humans (Hicks et al., 2019), and supporting people's livelihoods (Teh and Sumaila, 2013). Recent international assessments emphasize the predominately adverse impacts of climate change on marine ecosystems and fisheries (Bongaarts, 2019; IPCC, 2019, 2021). Projections indicate that ocean warming (Cheng et al., 2022) and deoxygenation (Frölicher et al., 2009; Kwiatkowski et al., 2020; Schmidtko et al., 2017) will likely result in a global decline in total animal biomass (Lotze et al. 2019) and potential fishery catch (Cheung et al. 2016), as well as alterations in species distributions (Deutsch et al., 2020; Howard et al., 2020; Hodapp et al., 2023; Morée et al., 2023; Mongwe et al., 2024). Despite considerable research on the transient impacts of global warming and deoxygenation throughout the 21st century on marine fisheries, there remains a considerable gap in understanding the long-term repercussions of centennial-scale warming and deoxygenation. Moreover, the increasing probability of a temporary temperature overshoot, despite efforts to limit global warming to the Paris Agreement's 1.5 °C long-term goal (Rogelj et al., 2018), raises critical concerns – particularly regarding our limited understanding of the potential reversible or irreversible impacts of such a temporary overshoot on marine species (Meyer et al., 2022; Meyer and Trisos, 2023).

To comprehensively assess the multi-centennial impacts on marine species and potential reversibility after a temporary temperature overshoot, policy-relevant multi-centennial temperature stabilization and overshoot simulations are needed but are currently lacking (King et al., 2021; Nature Geoscience Editorial, 2023). Recent progress has been made in understanding the physical and biogeochemical reversibility and hysteresis under idealized scenarios (Frölicher and Joos, 2010; Jeltsch-Thömmes et al., 2020; Schwinger et al., 2022). These studies are, however, not extended to species impacts and are not designed to achieve specific global warming levels (Keller et al., 2018). Modeling studies that consider impacts on global marine species are often limited to transient warming (i.e., a non-stabilized and therefore temporary warming and oxygen level, Cheung et al., 2016; Deutsch et al., 2015; IPCC, 2019; Morée et al., 2023) or idealized overshoot scenarios (Meyer et al., 2022; Santana-Falcón et al., 2023) and do not apply to specific species (Deutsch et al., 2020; Oschlies, 2021; Santana-Falcón et al., 2023). Moreover, modeling of species' distributions generally relies on the environmental state of the sea surface or sea floor only (Garciá Molinos et al., 2016; Heneghan et al., 2021; Hodapp et al., 2023; Meyer et al., 2022), thereby neglecting the crucial vertical dimension when estimating changes in marine species distributions (Duffy and Chown, 2017) and hence largely neglecting changes in the ocean interior for disturbances in marine species distributions.

Recent advances in climate scenarios and impact assessments on marine species under ocean warming and deoxygenation, such as the adaptive emission reduction approach (AERA; Terhaar et al., 2022) and the Aerobic Growth Index (AGI); Clarke et al., 2021; Morée et al., 2023), enable temperature-targeted modeling of the viability of contemporary distributions of marine species. The AGI indicates the potential habitat conditions that theoretically support the aerobic scope required for the growth of marine-water-breathing ectotherms, represented by the ratio of pO₂ supply to metabolic demand. AGI, when above the species-specific critical level (AGI^crit), signifies the potential for a viable population of the considered species to be sustained with respect to temperature and the partial pressure of oceanic oxygen, pO₂. This makes AGI a metric for habitat viability (Clarke et al., 2021; Morée et al., 2023). AERA achieves policy-relevant temperature stabilization at any desired level and temporary temperature overshoot at any desired magnitude and duration by iteratively adjusting CO₂-forcing-equivalent (CO₂-fe) emissions in Earth system models.

Here, we apply AERA to the Earth system model GFDL-ESM2M to simulate novel global warming scenarios spanning the period 1861 to 2500, encompassing stable temperatures of 1.2, 1.5, 2, and 3 °C as well as 2 and 3 °C overshoot scenarios followed by a return to stabilized 1.5 °C warming relative to the 1861–1900 baseline (Fig. 1a, b). We quantify changes in the contemporary habitat volume of 46 representative marine species that are important to fisheries across the epipelagic, mesopelagic, and demersal realms by utilizing AGI. Additionally, we examine the potential effects of species' adaptation to ocean warming and deoxygenation on their habitat volume loss, as adaptation is likely an important factor at the centennial timescale considered here (Pinsky et al., 2020).

2 Methods

2.1 Earth system model

The temperature stabilization and overshoot simulations were carried out with the fully coupled Earth system model GFDL-ESM2M (Dunne et al., 2012, 2013). The GFDL-ESM2M couples an ocean, atmosphere, land, and sea-ice model. Ocean physics, including sea ice, is simulated with the Modular Ocean Model version 4p1 at a 1 ° horizontal resolution, which increases to up to $\mathrm{1}/\mathrm{3}\mathit{°}$ meridionally at the Equator and is tripolar above 65 ° N, with 50 vertical levels (Griffies, 2012). Mesoscale eddies are parameterized. Marine biogeochemistry is represented by Tracers of Ocean Productivity with Allometric Zooplankton version 2 (TOPAZv2) and consists of the cycles of carbon, nitrogen, phosphorus, silicon, iron, O₂, alkalinity, and particulate matter as well as surface sediment calcite. Three phytoplankton groups and one zooplankton group are explicitly represented in TOPAZv2. Denitrification is included under suboxic conditions, while in the absence of both O₂ and nitrate, additional respiration is accumulated as a negative O₂ concentration. The atmospheric module (AM2) consists of 24 vertical levels and a horizontal resolution of 2 ° latitude ×2.5 ° longitude (Anderson et al., 2004). The land model (LM3.0) represents the hydrological, energy, and carbon cycles on land.

2.2 Model simulations with the adaptive emissions reduction approach

The model employed the AERA (Terhaar et al., 2022, 2023) to conduct novel simulations spanning 1861 to 2500 that stabilize global mean 2 m air temperature anomalies at 1.2, 1.5, 2, and 3 °C, including temporary overshoots to 2 and 3 °C and then returning to stable 1.5 °C warming relative to 1861–1900. The temperature changes were selected to align with the global warming levels commonly used in the IPCC reports – 1.5, 2.0, and 3.0 °C (IPCC, 2022) – as benchmarks for assessing impacts. Additionally, we included a stabilization scenario at the current level of warming (i.e., 1.2 °C) to investigate committed impacts. The warming rates leading to these temperature levels closely follow (at least initially) those observed during the historical period. The ocean physical and carbon cycle responses in those scenarios are discussed in Lacroix et al. (2024) and Silvy et al. (2024).

All simulations branch off from an emissions-driven simulation over the historical period from 1861 to 2005. Post-2005, fossil fuel CO₂ emissions follow observed emissions until 2020 and projected emissions from the nationally determined contributions (Climate Action Tracker: https://climateactiontracker.org/global/temperatures/, last access: 1 December 2021) from 2021 to 2025. From 2025 onward, prescribed fossil fuel CO₂ emissions are obtained every 5 years from AERA by subtracting CO₂-forcing-equivalent (CO₂-fe) emissions from prescribed non-CO₂ radiative agents and land use change from the AERA-derived total CO₂-fe emission curve. Non-CO₂ radiative forcing and land use change are prescribed to follow the RCP2.6 scenario until 2100, with no changes thereafter. The AERA approach adapts CO₂-fe emissions successively to converge to pre-defined warming levels. AERA calculates the realized simulated anthropogenic warming and the resulting remaining CO₂-fe emission budget every 5 years. This remaining emission budget is then distributed across future years using a cubic polynomial function. The temporary overshoot simulations are branched off from the 2 and 3 °C stabilization simulations, gradually reducing the warming level to 1.5 °C by circa 2300 following the approach described in Terhaar et al. (2022) and conducted in Lacroix et al. (2024) in GFDL-ESM2M. The resulting CO₂-fe emission pathways that forced the model to the different warming levels show a maximum change in annual CO₂-fe emissions of −2.1, −1.5, −0.7, and −0.8 Pg C yr⁻² for the 1.2, 1.5, 2, and 3 °C stabilization scenarios, respectively, and a maximum change in annual CO₂-fe emissions of −0.7 and −1.6 Pg C yr⁻² for the 2 and 3 °C overshoot scenarios (Fig. A1).

2.3 Model evaluation

We briefly summarize previous evaluations of GFDL-ESM2M's physical and biogeochemical state with a focus on the variables most relevant for this study. In general, the model tends to better simulate physical variables than biogeochemical variables. The simulated ocean overturning structure and rate, ocean heat transport, ocean temperature and salinity, sea surface height, water mass age distributions, and mixed-layer depth in general show fair agreement with observational data (Dunne et al., 2012). Some potentially relevant mean biases include a thermocline that is overall too deep, Antarctic Bottom Water formation through open ocean convection, weak deep Pacific ventilation, and North Pacific thermocline ventilation that is too strong (Dunne et al., 2012). Furthermore, GFDL-ESM2M performs within the CMIP5 range with respect to O₂ projections, despite the global mean O₂ concentration being 5 % lower than observed due to comparably large low-O₂ water volumes (Bopp et al., 2013). The spatial present-day distribution of O₂ agrees well with observations, and most of the bias in O₂ is limited to the oxygen minimum zones in the eastern equatorial Pacific at depth (Frölicher et al., 2020), similar to other Earth system models (Cabré et al., 2015). The globally integrated net primary production (NPP) is 74 PgC yr⁻¹ in GFDL-ESM2M compared to 53 Pg C yr⁻¹ on average in the observation-based estimates (Le Grix et al., 2022). GFDL-ESM2M overestimates NPP, especially in the low latitudes, but despite these differences, GFDL-ESM2M succeeds in simulating the NPP mean spatial pattern of higher values in the low latitudes and lower values in the Southern Ocean and the subtropical gyres. As a result of the mean temperature and oxygen biases in the model, all mean biases have been corrected for by subtracting the mean bias between World Ocean Atlas 2018 data and the model from all model data (see below).

A comparison of the relative changes in the Aerobic Growth Index (AGI_rel; see Sect. 2.5) between CMIP6 model simulations and GFDL-ESM2M at transient 1.5 °C warming allows for additional evaluation of GFDL-ESM2M's performance. In CMIP6 models (Morée et al., 2023), the multi-model mean AGI_rel changes are slightly smaller than in GFDL-ESM2M, while patterns are generally similar. We attribute part of the relatively weaker signal in the CMIP6 ensemble compared to GFDL-ESM2M to the later timing of reaching (transient) 1.5 ° warming in GFDL-ESM2M. Nevertheless, pelagic AGI_rel largely falls within the CMIP6 range, while demersal AGI_rel may be of larger magnitude than CMIP6.

2.4 Model data and analysis

The annual mean GFDL-ESM2M data used in this study consist of dissolved O₂ concentration, potential temperature, salinity, ideal age, and the particulate organic carbon flux at 100 m depth. The model data were first drift corrected by subtracting a linear trend obtained from the corresponding preindustrial control simulation over 1861–2500. The simulated O₂, temperature, and salinity data were bias corrected using World Ocean Atlas 2018 (WOA18) observational data (Garcia et al., 2019; Locarnini et al., 2019; Zweng et al., 2019) as follows: climatological mean WOA18 data were regridded horizontally and vertically to the model grid, after which the bias between the WOA18 data and the time-mean model data over the years 1995–2014 was subtracted from the full time series. In situ temperature is calculated from the potential temperature, and pO₂ is calculated from in situ temperature and the salinity and O₂ data (Bittig et al., 2018; Morée et al., 2023) .

The “first year of a stable warming level” of its respective global warming level refers to the year in which the 31-year running mean of the global mean 2 m air temperature is within 1 standard deviation of the aimed-for warming level (standard deviation of 0.074 °C obtained from the preindustrial control simulation). The 31-year mean was chosen to exclude interannual-to-decadal natural variability to a large degree. To make a fair comparison between the temperature stabilization scenarios with respect to the maximum amount of equilibration time available after reaching stable warming, committed habitat change is the change in the 296 years after hitting the respective stable warming level. Since the 3 °C stabilization simulation reaches its warming level the latest – specifically in the year 2162 – there are 296 years remaining until the end of the simulation in 2500. Consequently, we adopt this 296-year timescale for assessing committed changes in all stabilization scenarios. Habitat change in a certain year is expressed as the forward 31-year mean in that year.

For the species-specific results, we estimate the combined uncertainty (e.g., whiskers in Fig. 3) from the timing of first reaching a stable warming level (varied by ±10 years) and the species-specific control simulation variations in subcritical habitat volume, which are combined using the root sum square of the respective standard deviations. For the overshoot data, only the control standard deviation is considered as there is no stable warming level hit year.

2.5 Aerobic Growth Index

We applied two alternative formulations of AGI as indicators of species-level and ecosystem-level vulnerabilities to the impacts of ocean warming and deoxygenation (see Eqs. 1 and 2, respectively).

For each species i, the AGI is calculated as the ratio between the pO₂ supply and demand (Eq. 1) (Clarke et al., 2021; Morée et al., 2023):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{AGI}}_i& ={\displaystyle \frac{p{\mathrm{O}}_{\mathrm{2}}^{\mathrm{supply}}}{p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{demand}}}}\\ & ={\displaystyle \frac{p{\mathrm{O}}_{\mathrm{2}}^{\mathrm{supply}}}{p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}\cdot (\frac{\mathrm{1}}{\mathrm{3}}{)}^{\mathrm{1}-d}\cdot \exp \left(\frac{j_{\mathrm{2}}-j_{\mathrm{1}}}{T_i^{\mathrm{pref}}}-\frac{j_{\mathrm{2}}-j_{\mathrm{1}}}{T}\right)}},\end{array}\end{array}$$

with partial pressure pO₂ (units in Pa) and in situ temperature T (K). The generalized temperature dependence is represented by the variables j₁ (the anabolism activation energy divided by the Boltzmann constant, 4500 K) and j₂ (the catabolism activation energy divided by the Boltzmann constant, 8000 K), scaled by the metabolic scaling coefficient d (0.7). Aerobic scope limits the species' distribution. Therefore, by using biogeographical data, we inferred the oxygen threshold necessary to support a viable population of the species. The species-specific critical threshold of pO₂ ($p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}$) is calculated as the volume-weighted 10th percentile of 1995–2014 time-mean pO₂ within the three-dimensional distribution of the species. The estimated pO₂ threshold for growth is proportional to critical pO₂ measured from physiological experiments (Clarke et al., 2021). The preferred in situ temperature ($T_i^{\mathrm{pref}}$) is calculated as the volume-weighted 50th percentile of the 1995–2014 time-mean temperature within the three-dimensional distribution of the species and ${\mathrm{AGI}}_i^{\mathrm{crit}}$ as the volume-weighted 10th percentile of 1995–2014 time-mean AGI (Clarke et al., 2021; Morée et al., 2023).

Using AGI_i and ${\mathrm{AGI}}_i^{\mathrm{crit}}$, we quantify the percentage of contemporary habitat available to sustain a viable population of the species (i.e., where ${\mathrm{AGI}}_i>{\mathrm{AGI}}_i^{\mathrm{crit}}$). Our use of this extended version of AGI (Morée et al., 2023) includes vertical variability in O₂ and temperature in the estimate of species' $p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}$, $T_i^{\mathrm{pref}}$, and hence AGI_i and ${\mathrm{AGI}}_i^{\mathrm{crit}}$, as well as considering fully three-dimensional changes in O₂ and temperature (as opposed to surface or sea floor values only). The generalized temperature dependence of $p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{demand}}$ may cause over- or underestimation of an individual species' response to temperature. The sensitivity of loss of contemporary habitat to this parameterization (j₂–j₁ in Eq. 1) is minor (Morée et al., 2023).

Changes in AGI between a time t₁ and t₀ relative to t₀ are species independent and provide a sense of direction and magnitude of change in habitat viability (Morée et al., 2023):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathrm{AGI}}^{\mathrm{rel}}& ={\displaystyle \frac{\mathrm{\Delta }\mathrm{AGI}}{\mathrm{AGI}\left(t_{\mathrm{0}}\right)}}={\displaystyle \frac{p{\mathrm{O}}_{\mathrm{2}}\left(t_{\mathrm{1}}\right)}{p{\mathrm{O}}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)}}\\ & \cdot \exp \left(\left(j_{\mathrm{2}}-j_{\mathrm{1}}\right)\cdot \left({\displaystyle \frac{\mathrm{1}}{T\left(t_{\mathrm{1}}\right)}}-{\displaystyle \frac{\mathrm{1}}{T\left(t_{\mathrm{0}}\right)}}\right)\right)-\mathrm{1}.\end{array}\end{array}$$

Here, t₀ refers to the mean over the period 1861–1900 in this study. Note that small changes in pO₂ can lead to large changes in AGI_rel if reference pO₂ (i.e., pO₂ at t=0) is low (Eq. 2).

We recognize the potential habitat volume, and the reversibility of its changes may not be realized by the species because of other biogeographic constraints, such as dispersal potential, availability of suitable prey, or other environmental limitations beyond temperature and oxygen.

2.6 Species data

Spatial distribution data for the 46 representative exploited species (species names are indicated in Fig. 5; Morée et al., 2023; Palomares et al., 2004) are used to calculate $T_i^{\mathrm{pref}}$, $p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}$, and ${\mathrm{AGI}}_i^{\mathrm{crit}}$. These data also form the reference habitat for assessing changes in contemporary habitat volume. The species were selected such that they provide a representative range in body size, climatic zone (tropical, temperate), habitat size, and depth range. The 46 species also cover a broad range of vulnerabilities to warming and deoxygenation, with the most vulnerable species having a ∼30 times larger change in volume per unit change in AGI than the least vulnerable species (Morée et al., 2023). We include 23 species with their predominant occurrence in the epipelagic (0–200 m depth), 5 species that mostly inhabit the mesopelagic (200–1000 m depth), and 18 demersal species which live on or just above the sea floor (for which we use the deepest ocean model layer). Some pelagic and deep-water wide-ranging species were selected that inhabit both tropical and temperate regions. We assess contemporary habitat volume change by extending the two-dimensional habitats over their depth range assuming that the distribution is the same throughout the water column. We acknowledge that some species may occupy only part of their assigned depth range or may temporarily reside outside it, either above or below. Nevertheless, we believe that the assigned depth ranges generally provide a reasonable estimate of in-habitat pO₂ and temperature variability. The reason for not just using the three-dimensional model output across the entire depth range of each species' depth realm is primarily due to the current lack of reliable three-dimensional species distributions for our selected species.

All habitat changes are changes in contemporary habitat volume since the 1861–1900 mean, to be consistent with the atmospheric temperature anomalies, as well as being presented as a percentage of viable contemporary habitat volume to facilitate comparison between species.

2.7 Adaptation

Species will likely acclimatize and/or adapt to global warming by improving their temperature and O₂ stress resistance on the timescales considered in this study (Pinsky et al., 2020). Such acclimatization and adaptation to environmental changes have already been observed over the last few decades, for example for warm-water corals (Lachs et al., 2023; Logan et al., 2021). We base our adaptability assessment on a running climatology approach, which allows for the evaluation of both adaptation pressure to a species' habitat and the timescale of a species' adaptation to such pressure. We thereby follow the approach applied to corals by Logan et al. (2014), as outlined in the paragraphs below.

Adaptation pressure is evaluated by considering for each species a time-dependent ${\mathrm{AGI}}_i^{\mathrm{crit}}$ (${\mathrm{AGI}}_{i,t}^{\mathrm{crit}})$ that changes when in-habitat AGI (and hence pO₂ and/or temperature) changes. As AGI_i and thus ${\mathrm{AGI}}_{i,t}^{\mathrm{crit}}$ depend on $T_i^{\mathrm{pref}}$ and $p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}$, we also calculate a time-dependent $p{\mathrm{O}}_{\mathrm{2},i}^{\mathrm{threshold}}$ ($p{\mathrm{O}}_{\mathrm{2},i,t}^{\mathrm{threshold}}$) and $T_i^{\mathrm{pref}}(T_{i,t}^{\mathrm{pref}}$). We do so by following the standard approach (described in Sect. 2.5) but instead using a moving 20-year time-mean time window of pO₂, temperature, and AGI before each year t. In this way, ${\mathrm{AGI}}_{i,t}^{\mathrm{crit}}$ is calculated for each species and each year over the period 1861 to 2500. It also follows that ${\mathrm{AGI}}_{i,t=\mathrm{2014}}^{\mathrm{crit}}$ is equal to the standard ${\mathrm{AGI}}_i^{\mathrm{crit}}$. For the first 20 years of our dataset, we use the time mean of all years available before and including that particular year.

We then consider the timescale of adaptation and its uncertainty by using mean ${\mathrm{AGI}}_{i,t}^{\mathrm{crit}}$ over 40, 60, 80, and 100 years (the range in Logan et al., 2014), as well as 140, 180, and 220 years (to extend to species with longer adaptation times due to, for example, slower generational overturn) before year t in our quantification of contemporary habitat volume change in year t. A 40-year adaptation timescale thus means that a species has adjusted (i.e., adapted) its ${\mathrm{AGI}}_{i,t}^{\mathrm{crit}}$ to the mean ${\mathrm{AGI}}_{i,t}^{\mathrm{crit}}$ over the previous 40 years.

3 Results

3.1 Committed loss of contemporary habitat and overshoot impacts

In the stabilization simulations, the global atmospheric surface warming levels of 1.2, 1.5, 2, and 3 °C are reached in the years 2019, 2054, 2091, and 2162, respectively (Fig. 1a). Afterwards, atmospheric temperatures are stabilized. Concurrently, the ocean is projected to warm and lose oxygen (Fig. 1c, e). By 2400–2500, the habitats of the 46 representative marine species have warmed by a median 0.8, 1.1, 1.5, and 2.2 °C (Fig. 1c) and lost 1.2, 1.3, 1.8, and 3.1 mbar of pO₂ in 2400–2500 for the 1.2, 1.5, 2, and 3 °C stabilization scenarios, respectively (median across all species; Fig. 1e).

Figure 1 Simulated temperature stabilization and overshoot scenarios and associated impacts on median in-habitat temperature and pO₂ changes and on contemporary habitat volume across 46 representative marine species. (a, b) Global mean 2 m air temperature change since 1861–1900 for (a) temperature stabilization scenarios at 1.2, 1.5, 2, and 3 °C warming and (b) overshoot scenarios that peak at 2 and 3 °C warming before returning to 1.5 °C warming. (c–h) Median change in (c, d) in situ temperature and (e, f) pO₂ in the respective habitats of the 46 species and in (g, h) contemporary habitat volume (% of volume) since 1861–1900 for the (c, e, g) stabilization and (d, f, h) overshoot scenarios. The whiskers in (c)–(h) indicate the 25th to 75th percentile range across the species habitat at the respective year of the starting temperature stabilization and in year 2500. The thick lines show the median across the individual species' 31-year running mean data and thin lines the annual mean time series.

With habitat warming and deoxygenation, species' habitat volume, defined here as the volume of contemporary habitat with the Aerobic Growth Index (AGI) above a species-specific threshold (see Methods), shrinks across the 46 species under all stabilization scenarios (Fig. 1g). Contemporary habitat volume continues to decrease even after atmospheric temperatures have stabilized. As a result, only about half of the total loss of contemporary habitat volume is realized when first reaching the respective stable warming levels: the median fraction of realized habitat loss upon hitting stable warming levels, compared to the loss 296 years thereafter, amounts to 0.45 (with 25th–75th percentiles across the species at 0.19–0.76) for the 1.2 °C scenario, 0.41 (0.25–0.67) for the 1.5 °C scenario, 0.67 (0.36–0.93) for the 2 °C scenario, and 0.60 (0.40–0.81) for the 3 °C stabilization scenario (Table 1). Adverse impacts on habitat viability of marine species thus continue to worsen for centuries despite stabilization of the global mean atmospheric temperature.

Table 1 Simulated median habitat change (ΔHabitat) across the 46 representative species, presented in percentage of contemporary habitat volume at (a) the time of warming level hit and (b) 296 years after the time of warming level hit. (c) The difference between (a) and (b) and (d) the ratio between (a) and (b). In all columns the 25th and 75th percentiles across the species are indicated in brackets, and all values are 31-year forward running means.

After 296 years of temperature stabilization, the median decline in contemporary habitat among the species ranges from 1.06 % (with 25th–75th percentiles across the species at 2.32 % to 0.89 %) for the 1.2 °C stabilization scenario to 4.59 % (7.13 % to 3.14 %) for the 3 °C stabilization scenario (Fig. 1g; Table 1). While these changes may appear small, the corresponding volumetric losses are substantial, ranging from thousands (order for the demersal species) to millions of cubic kilometers (order for the epipelagic species), depending on the species and scenario (Morée et al., 2023). Notably, the 25 % of species with the largest losses experience 2–5 times larger losses depending on the species and warming level (whiskers in Fig. 1g).

Temporarily overshooting the 1.5 °C warming level to 2 °C or 3 °C (Fig. 1b) drives bigger loss of contemporary habitat volume than direct stabilization at 1.5 °C (Fig. 1h) due to associated overshoots of in-habitat warming (Fig. 1d) and deoxygenation (Fig. 1f). Peak loss of contemporary habitat volume during the 2 and 3 °C overshoot scenario is −2.1 % and −3.2 % in the years 2270 and 2246, respectively (medians across the species, Table 2). This corresponds to −0.5 % and −1.7 % more loss of contemporary habitat than under the 1.5 °C stabilization scenario at that time. Notably, peak loss of contemporary habitat is realized 166 years (for the 2 °C overshoot scenario) and 66 years (3 °C overshoot) after maximum atmospheric warming in the years 2104 and 2180 (Fig. 1b, h). By the time the maximum median loss of contemporary habitat is realized, atmospheric warming has therefore already returned to 1.7 and 2.5 °C, respectively. When atmospheric warming has returned to 1.5 °C in 2400–2500, median contemporary habitats are 0.06 % larger (0.31 % smaller) in the 2 °C (3 °C) overshoot scenario than in the 1.5 °C stabilization scenario. This indicates that for most species, overshoot impacts are largely reversible by this time (Table 2). Nevertheless, the strong adverse impacts during the overshoot underscore the advantages of stabilizing at 1.5 °C warming over overshooting this level.

Table 2 Simulated median habitat change (ΔHabitat) across the 46 representative species, presented in percentage of contemporary habitat volume (a) in the year of maximum change for the respective scenario and (b) in 2400–2500. (c) The difference with the 1.5 °C stabilization scenario in 2400–2500 and (d) the ratio relative to the 1.5 °C stabilization scenario in 2400–2500. In (d), the overshoot legacy is the difference between the overshoot scenario and the 1.5 °C warming scenario relative to the change in the 1.5 °C warming scenario. In all columns the 25th and 75th percentiles across the species are indicated in brackets, and all values are 31-year forward running means.

3.2 Spatial patterns and drivers of habitat viability changes

To generalize these median impacts beyond the 46 analyzed species and to investigate the drivers of changes, we calculate the local relative changes in habitat viability expressed as the changes in the Aerobic Growth Index relative to the period 1861–1900 (AGI^rel). These relative changes in habitat viability are fully species independent and indicate the direction and extent of change in habitat viability across the world's oceans. We focus on the 1.5 °C stabilization and the 2 °C overshoot scenario, but results are qualitatively similar for the other scenarios.

Figure 2 Simulated committed change in the 1.5 °C stabilization scenario (a, b, c) and overshoot legacy in the 2 °C overshoot scenarios (d, e, f) in relative changes in habitat viability for three depth realms. (a–c) Committed change in habitat viability is expressed as the difference between AGI^rel 296 years after temperatures have stabilized and when the 1.5 °C warming level is hit (i.e., the committed period). (d–f) Overshoot legacy in habitat viability is expressed as the difference between AGI^rel in the 2 °C overshoot and the 1.5 °C stabilization scenarios in 2400–2500. Results are shown as the vertical mean over the epipelagic (0–200 m), mesopelagic (200–1000 m), and demersal (sea floor) realms.

In the 1.5 °C stabilization scenario, committed decreases in AGI^rel (i.e., after temperatures have stabilized) intensify with depth and are comparable in magnitude to the changes during the period until 1.5 °C warming is first reached (Fig. 2a–c). The most substantial committed declines in the epipelagic realm occur in the North Atlantic (Fig. 2a), driven by a sustained rise in ocean temperature (Fig. 3a) linked to the recovery of the Atlantic Meridional Overturning Circulation during the stabilization phase (Frölicher et al., 2020; Lacroix et al., 2024). Elsewhere in the pelagic and demersal realms, committed changes in habitat viability are predominantly driven by O₂ changes through changes in apparent oxygen utilization (i.e., changes in ocean ventilation and organic matter remineralization) (Figs. 3 and A2). In the Pacific demersal realm (Fig. 2c), there is a committed reduction in habitat viability due to O₂ decrease tied to reduced deep ocean ventilation, as indicated by older water ages (Fig. 3). Conversely, there is a committed increase in habitat viability in the North Pacific mesopelagic realm and Southern Ocean mesopelagic and demersal realms. The North Pacific increase is linked to increased ventilation and reduced biological oxygen consumption (Figs. 2 and A2). The Southern Ocean increase is mainly driven by enhanced ventilation, as indicated by younger water ages (Fig. 2).

Figure 3 Simulated committed changes in the environmental state for the 1.5 °C stabilization scenario for three depth realms. (a–c) Committed change in in situ temperature, (d–f) pO₂, (g–i) AGI^rel due to changes in O₂ only (calculated by keeping temperature constant at its 1861–1900 mean), (j–l) apparent oxygen utilization (AOU $={\mathrm{O}}_{\mathrm{2}}^{\mathrm{sat}}-{\mathrm{O}}_{\mathrm{2}}$), and (m–o) ideal age of the water. As a proxy for change in ventilation time, we used the ideal age tracer, which is set to zero in the surface ocean and ages at a rate of 1 yr yr⁻¹ below that. Committed change is the difference between the environmental state 296 years after temperatures have stabilized and when the 1.5 °C warming level is hit. Results are shown as the vertical mean over the epipelagic (0–200 m), mesopelagic (200–1000 m), and demersal (sea floor) realms.

The overshoot legacy for the 2 °C overshoot is determined by comparing the changes in the overshoot scenario with the 1.5 °C stabilization scenario without overshoot in the period 2400–2500, when global warming levels are similar (Fig. 2d–f). The most pronounced overshoot legacies are simulated in the northern flank of the eastern equatorial Pacific epipelagic (Fig. 2d) and in the mesopelagic of the North Pacific (Fig. 2e). These regions of habitat viability increase are caused by higher O₂ levels resulting from reduced O₂ consumption and rejuvenating mesopelagic waters (Fig. 4j, k), which is consistent with previous idealized model simulations (Li et al., 2020; Santana-Falcón et al., 2023). In the demersal realm, relative changes in habitat viability are generally more negative under the overshoot scenario than in the stabilization scenario by 2400–2500 (Fig. 2f). This negative legacy is mainly driven by reduced deep ocean ventilation, as indicated by increased water ages (Fig. 4).

Figure 4 Simulated overshoot legacy for the 2 °C overshoot in the environmental state for three depth realms in 2400–2500. (a–c) Overshoot legacy in in situ temperature, (d–f) pO₂, (g–i) AGI^rel due to changes in O₂ only (calculated by keeping temperature constant at its 1861–1900 mean), (j–l) apparent oxygen utilization (AOU $={\mathrm{O}}_{\mathrm{2}}^{\mathrm{sat}}-{\mathrm{O}}_{\mathrm{2}}$), and (m–o) ideal age of the water. As a proxy for change in ventilation time, we used the ideal age tracer, which is set to zero in the surface ocean and ages at a rate of 1 yr yr⁻¹ below that. Overshoot legacy is the difference between the state in the overshoot scenario and 1.5 °C stabilization scenario in 2400–2500. Results are shown as the vertical mean over the epipelagic (0–200 m), mesopelagic (200–1000 m), and demersal (sea floor) realms.

3.3 Impacts on contemporary habitat volume of individual species

The committed loss of contemporary habitat volume is generally the largest for demersal species (median of 1.7 % in the 1.5 °C stabilization scenario and up to 6.4 % for Spectrunculus grandis) and typically the smallest in the epipelagic realm (median of 0.7 %) (gray bars in Fig. 5a). When reaching stable 1.5 °C warming in 2054, only a median 47 % of the total loss is realized in the epipelagic realm, 56 % in the mesopelagic realm, and 28 % in the demersal realm (Fig. A3a). Both the committed loss of contemporary habitat volume and the contrast in the impact between the different depth realms generally strengthen with increased levels of global warming (Fig. A4). For most species, the committed loss of contemporary habitat volume is driven by ongoing decreases in dissolved O₂ (black stars in Fig. 5a). However, the magnitude of loss of contemporary habitat for individual species also depends on species-specific vulnerability and not only on the magnitude of changes in AGI^rel, temperature, or O₂ (Morée et al., 2023).

Figure 5 Simulated committed changes in contemporary habitat volume for the 1.5 °C stabilization scenario (top panels) and overshoot legacy in the 2 °C overshoot scenario (bottom) for 46 representative marine species and three different depth realms. (a) Committed change (gray bars) expressed as the difference between contemporary habitat volume changes 296 years after temperatures have stabilized and when the 1.5 °C warming level is hit (i.e., the committed period). (b) Overshoot legacy (red bars) expressed as the difference between the 2 °C overshoot scenario and the 1.5 °C stabilization scenario in 2400–2500. The black stars in (a) and (b) indicate contemporary habitat changes that are driven by temperature changes only (i.e., keeping O₂ values at 1861–1900 conditions). Maximum changes during the overshoot are indicated with the horizontal red bars in (b). The whiskers in (a) indicate the uncertainty as the combined uncertainty coming from the uncertainty in the timing of warming level hit and the species-specific control simulation variability (in b, only the species-specific control simulation variability is considered in the whiskers as no stable warming level hit year is considered; see Methods).

When temporarily overshooting stable 1.5 °C warming to 2 °C, median peak overshoot habitat loss exceeds the impact of 1.5 °C warming by 1.1 % and at most 11.4 % for Thunnus atlanticus (bars in Fig. 5b; see Fig. A4 for the 3 °C overshoot). After overshoot reverses in 2400–2500 from 2 °C to stable 1.5 °C warming, the species contemporary habitat volume is between −2.6 % and +2.3 % of stable 1.5 °C warming (Fig. 5b), with a median loss of contemporary habitat volume of 0.10 % in the epipelagic, 0.11 % in the mesopelagic, and 0.36 % in the demersal realms (medians across depth realms in Fig. 5b). The overshoot legacy is predominantly explained by enhanced (de)oxygenation in the overshoot simulation relative to the stabilization simulation in 2400–2500 (stars in Fig. 5b). Some species like Thunnus atlanticus and Hippoglossus hippoglossus even benefit from overshooting 1.5 °C warming, mainly due to local oxygenation relative to the 1.5 °C stabilization scenario.

3.4 Potential adaptation

We consider both the adaptation pressure to a species which would incite adaptation (as the in-habitat change in habitat viability quantified by AGI) and the timescale of adaptation (the number of years needed to adapt a species' thresholds, evaluated for 40, 60, 80, 100, 140, 180, and 220 years). Consequently, after being under high (low) adaptation pressure, a species will be relatively insensitive (sensitive) for a time that increases with the timescale of adaptation.

Adaptation substantially reduces the projected loss of contemporary habitat compared to non-adaptation (Fig. 6). In the 2 °C stabilization scenario that includes adaptation, loss of contemporary habitat temporarily peaks at 0.63 %–1.53 % around 2100 depending on the timescale of adaptation, before gradually recovering to 0.26 %–0.40 % loss by 2500 (Fig. 6a, qualitatively similar results are obtained for the other scenarios). In comparison, the loss of contemporary habitat without adaptation is 2.54 % by 2500. The maximum loss of contemporary habitat is larger and later for longer adaptation timescales. However, after the 21st century, loss of contemporary habitat in the adaptation scenarios is reduced again because critical thresholds in AGI have had time to decrease enough (i.e., species could adapt to tolerate less favorable conditions).

In the overshoot scenarios, adaptation also causes reduced impact when compared to the non-adaptive scenarios (Fig. 6b). By 2400–2500, adaptation reduces the 3 °C (2 °C) overshoot legacy in terms of loss of contemporary habitat to 0.04 % (0.28 %), while this is 2.11 % (1.75 %) without adaptation. Following our approach, we find that the stronger the overshoot a species is exposed to, the lower the species' critical threshold becomes through adaptation (i.e., adapting to lower AGI^crit and hence becoming more tolerant to low-oxygen or high-temperature environments). This results in a “double overshoot” in the 3 °C overshoot simulation, where habitat volume impacts after ∼2250 are smaller than in an adaptive 1.5 °C stabilization scenario because species are still adapted to the (much) less favorable conditions in this recovery phase (Fig. 6b). In all scenarios with adaptation, median impact on contemporary habitat volume does not reduce to zero because some species are exposed to a continued decrease in habitat viability even at the end of the 25th century.

4 Discussion and conclusions

Through novel policy-relevant stabilization scenarios using an Earth system model, we show that only about half of the fish habitat changes have been realized when the temperature stabilization level is first hit. In overshoot scenarios, it may take over 150 years after the peak temperature for the most substantial impacts on marine species to manifest. Additionally, our findings suggest that rapid adaptation to changing conditions by species could mitigate some of the effects, contingent upon the rate of adaptation and external pressure.

An important assumption in our analysis is that the employed single Earth system model GFDL-ESM2M simulates changes in temperature and dissolved O₂ in a sufficiently realistic manner. We consider our global analysis for the warming and deoxygenation to be robust, especially given the good agreement of the model with observations and with other Earth system models (see Methods). However, projected local changes in oxygen may be less robust (Kwiatkowski et al., 2020), especially in low-O₂ zones, and are likely to be underestimated (Buchanan and Tagliabue, 2021). Additional ESMs should be run under policy-relevant temperature stabilization and overshoot scenarios to quantify the uncertainty in pO₂ changes. Nevertheless, multiple modeling studies indicate a centennial to multimillennial response timescale of O₂ to global warming (Battaglia and Joos, 2018; Bertini and Tjiputra, 2022; Frölicher et al., 2020), making it likely that multi-centennial O₂ changes are indeed important for committed impacts and overshoot legacies. Therefore, our foremost conclusion – that a major part of the total loss of contemporary habitat volume occurs after global temperatures have stabilized – likely does not depend on the model used, but the quantitative results may vary when using other models.

We highlight a multi-decadal to centennial delay in peak impact in contrast to peak global warming in the overshoot scenarios. This is much longer than the 8-year delay found in Meyer et al. (2022) for marine species in a 2 °C overshoot. Their use of sea surface data likely causes an underestimation of the peak delay due to faster reversibility of temperature and O₂ at the surface as compared to depth (Santana-Falcón et al., 2023; Schwinger et al., 2022), providing an additional argument for using three-dimensional species distributions (Duffy and Chown, 2017). The depth of a species' distribution is also central to the quantification of the legacy impact of an overshoot after return to stable 1.5 °C warming, since changes in ocean conditions in greater depths take substantially longer to reverse following overshoot scenarios (Lacroix et al., 2024; Santana-Falcón et al., 2023; Schwinger et al., 2022). Furthermore, the steep vertical gradients of pO₂ and temperature as well as the dampened climate change signal at depth provide additional arguments to include the third, vertical dimension when estimating (changes in) contemporary species distributions (Duffy and Chown, 2017).

The Aerobic Growth Index faces limitations in its application due to several factors, notably the absence of species-specific fully three-dimensional distribution data, its neglect of critical stressors beyond warming and deoxygenation, and its confinement to assessing solely the loss of contemporary habitat. While we aim to capture habitat variability in O₂ and temperature relevant for species-specific critical thresholds by extrapolating two-dimensional habitats across different depth ranges, these approximations may not fully capture the complete picture. Additionally, the response of marine organisms to global warming may encompass sensitivities to other effects, such as changes in low temperatures, acidity, nutrient availability, phenology, disease, predation pressure, invasive species, and (over)fishing (Gissi et al., 2021; IPCC, 2023). Multi-stressor research aims to assess and predict the cumulative effects of these stressors, including synergies and antagonistic relationships between them, but is still bound by many challenges before full assessments can be made (Gissi et al., 2021; Orr et al., 2020). The impacts of changing temperature and oxygen levels on species' habitat volume may be underestimated under the overshoot scenarios. Specifically, species may not be able to shift to viable habitats because of biogeographic constraints not represented by the AGI, e.g., dispersal potential and trophic interactions. Such uncertainties are particularly notable for species with large transient changes in habitat volume, such as blackfin tuna (Thunnus atlanticus), and those with relatively lower dispersal potential, such as coral hind (Cephalopholis miniata). Furthermore, the AGI, in its current form, only represents the limitation of temperature and oxygen on the warm-temperature edge of fish distributions and does not represent the observed reduction in aerobic scope at the low-temperature edge (Pörtner, 2010; Clarke et al., 2021; Deutsch et al., 2020). Thus, our results do not include habitat expansion or contraction due to ocean warming or cooling at the cold edge of species distributions. Development and application of AGI that incorporate low-temperature limitation of species' metabolism and growth will help quantify such uncertainties. Last, we note the importance of knowledge of species-specific critical pO₂ thresholds and preferred temperature environments. Relative changes in pO₂ supply to pO₂ demand ratios have been implied to assess ecosystem-level impact (Battaglia and Joos, 2018; Deutsch et al., 2015; Oschlies, 2021; Santana-Falcón et al., 2023), but species-specific thresholds and preference windows are needed for such estimates (Morée et al., 2023). In the future, the application of the AGI to more species will allow for the assessment of a wider range of interspecies responses.

Figure 6 The impact of potential adaptation on simulated change in contemporary habitat volume across 46 representative marine species in the 2 °C stabilization (a) and 2 and 3 °C overshoot (b) scenario. (a) Median habitat change since 1861–1900 for different adaptation timescales and the non-adaptive 2 °C stabilization scenario. (b) Median habitat change since 1861–1900 for the 2 and 3 °C overshoot scenarios and the 1.5 °C stabilization scenario assuming a 100-year adaptation timescale, as well as their non-adaptive counterparts. In both panels, whiskers indicate the 25th and 75th percentile limits across the species for the 40-year adaptation timescale in the year 2091 when first reaching the 2 °C warming level and in the year 2500. Time series are smoothed with a 31-year running mean.

Our sensitivity analysis shows that adaptation has large potential to mitigate the impacts of global warming on the habitat viability of marine fishes, although observations of a reduced equatorward extent of species' habitats suggest that adaptation may be too slow, particularly for larger species such as marine fishes (Hastings et al., 2020). Furthermore, we strengthen the notion that both the pathway and the magnitude of temperature overshoot are important for future ecological impact (Meyer and Trisos, 2023) and additionally show that warming pathways can further modulate impact through their effect on adaptation pressure. However, there is still a large gap in knowledge about marine fishes' adaptation to warming and deoxygenation. Particularly, the rate of adaptation to warming and deoxygenation would depend on many biological and ecological factors that are not considered in this study, e.g., the existing genetic diversity, species' life history traits, and population/meta-population structure and connectivity. In addition, one of many critical questions regarding how adaptation shapes biodiversity and biogeography under global change is the role of extreme events and tipping points in species adaptation. However, these dynamics remain challenging to study due to the scarcity of data on tipping events and evolutionary responses (Grant et al., 2017). Our analysis represents a first step in exploring the potential implications of adaptation for marine biogeography under centennial-scale climate change. Our findings could inspire future research to delve deeper into these complex and pressing issues.

Commonly used CMIP6 overshoot scenarios projected a rapid rise to around 2 °C of global warming, followed by reversal within 20–50 years (Pfleiderer et al., 2024). However, such a rapid reversal, requiring substantial atmospheric CO₂ removal, is unlikely given current limitations in carbon dioxide removal technologies (Schleussner et al., 2024). In our 3 °C overshoot scenario, where we assume that CO₂ removal occurs over more than 150 to 200 years, negative CO₂-fe emissions of up to 9 Pg C yr⁻¹ would still be required to bring temperatures back to 1.5 °C (Fig. A1). This scale of removal far exceeds what is achievable in the foreseeable future (Fuss et al., 2018). Even in the 2 °C overshoot scenario, the peak negative CO₂ emissions would need to reach approximately 3 Pg C yr⁻¹. We conclude that losses of marine species' contemporary habitats continue for centuries beyond reaching stable global warming levels and after peaking global warming in a temporary overshoot. Any impact assessed at transient warming levels (Hausfather et al., 2022; Morée et al., 2023) thus largely underestimates the total impact on marine species.

Appendix A: Additional figures

Figure A1 Prescribed annual CO₂-forcing-equivalent emissions in the GFDL-ESM2M simulations. (a) CO₂-fe emissions in the temperature stabilization scenarios. (b) CO₂-fe emissions in the temperature overshoot scenarios and the reference 1.5 °C stabilization scenario.

Figure A2 Simulated committed changes in the 1.5 °C stabilization (a) and overshoot legacy in the 2 °C overshoot scenario (b) of particulate organic carbon export flux at 100 m depth. (a) Committed changes are expressed as the ratio between the change during the committed period (the 296 years after stable warming level hit) and at the time when the 1.5 °C warming level is first hit. (b) Overshoot legacy of the 2 °C overshoot scenario as the difference between the 2 °C overshoot and the 1.5 °C stabilization scenario divided by the 1.5 °C stabilization scenario export flux in 2400–2500.

Figure A3 Simulated committed changes in the 1.5 °C stabilization scenario (a) and overshoot legacy in the 2 °C overshoot scenario (b) expressed in relative changes in contemporary habitat volume for 46 representative marine species and for three different depth realms. (a) Relative committed change as the ratio of habitat change at the time of stable warming level hit to the absolute total change in habitat after stabilization 296 years later. Values of −100 % or +100 % (−100 % is highlighted by a horizontal gray line) indicate that all change is realized at the time of reaching the warming level. (b) Relative overshoot legacy (red bars) is expressed as the ratio between the 2 °C overshoot scenario and the 1.5 °C stabilization scenario in 2400–2500. Hippoglossus hippoglossus and Theragra chalcogramma have a relative overshoot legacy of 165 % and 199 %, respectively. The black stars in (a) and (b) indicate contemporary habitat changes that are driven by temperature changes only (i.e., keeping O₂ values at 1861–1900 conditions; see Methods). The whiskers in (a) indicate the uncertainty as the combined uncertainty coming from the uncertainty in warming level hit timing and the species-specific control simulation variability (in b, only the species-specific control simulation variability is considered in the whiskers as no warming level hit year is considered; see Methods).

Figure A4 Simulated changes in contemporary habitat volume (%) for 46 representative marine species for different stabilization (a) and overshoot (b) scenarios. (a) Committed change expressed as the difference between contemporary habitat volume changes 296 years after temperatures have stabilized and when the 1.5 °C warming level is hit (i.e., the committed period). (b) Overshoot legacy as the difference between the 2 or 3 °C overshoot and 1.5 °C stabilization scenario in 2400–2500. Maximum changes during the overshoot (10-year running mean filter) are indicated with horizontal bars. Whiskers in (a) indicate the uncertainty as the combined uncertainty coming from the uncertainty in warming level hit timing and the species-specific control simulation variability (in b, only the species-specific control simulation variability is considered in the whiskers as no warming level hit year is considered; see Methods).