Abstract
The ability to anticipate marine habitat shifts responding to climate variability has high scientific and socioeconomic value. Here we quantify interannual-to-decadal predictability of habitat shifts by combining trait-based aerobic habitat constraints with a suite of initialized retrospective Earth System Model forecasts, for diverse marine ecotypes in the North American Large Marine Ecosystems. We find that aerobic habitat viability, defined by joint constraints of temperature and oxygen on organismal energy balance, is potentially predictable in the upper-600 m ocean, showing a substantial improvement over a simple persistence forecast. The skillful multiyear predictability is dominated by the oxygen component in most ecosystems, yielding higher predictability than previously estimated based on temperature alone. Notable predictability differences exist among ecotypes differing in temperature sensitivity of hypoxia vulnerability, especially along the northeast coast with predictability timescale ranging from 2 to 10 years. This tool will be critical in predicting marine habitat shifts in face of a changing climate.
Similar content being viewed by others
Introduction
Habitat shifts of marine species are closely associated with changes in ocean temperature (\(T\)) and dissolved oxygen (\({{{{{{\rm{O}}}}}}}_{2}\)) concentrations, as they are two fundamental constraints on aerobic metabolism of marine species and consequently determine the viability of marine habitats1,2,3,4,5,6. Physically, warming of the oceans decreases \({{{{{{\rm{O}}}}}}}_{2}\) levels directly as its solubility decreases in warmer waters, and indirectly as it intensifies upper ocean stratification, reducing ocean ventilation - thus suppressing the supply of oxygen into the ocean interior5,7,8. On the other hand, metabolic rates of organisms are \(T\)-dependent and increase exponentially with rising temperatures - thus leading to increased aerobic demand for \({{{{{{\rm{O}}}}}}}_{2}\) (e.g., ref. 9). For marine habitats to be metabolically viable, the environmental supply of \({{{{{{\rm{O}}}}}}}_{2}\) must meet not only the minimum physiological requirements for survival but also additional requirements for essential ecological activities of marine species such as feeding, defense, and reproduction3. This implies that \({{{{{{\rm{O}}}}}}}_{2}\) concentrations may be inadequate to meet the demand of species for oxygen under warmer conditions, probably resulting in geographical or vertical viable habitats shrinking or shifting. Because \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) exhibit strong decadal variability, predictions of the corresponding marine habitat shifts are potentially possible on multi-annual time scales.
A physiologically mechanistic framework has been developed recently that enables quantifying habitat constraints for different species arising from the metabolic dependence on \(T\) and requirements for \({{{{{{\rm{O}}}}}}}_{2}\)3,10. Under this framework, an index termed the Metabolic Index (\(\varPhi\)) is defined as the ratio of \({{{{{{\rm{O}}}}}}}_{2}\) supply to an organism’s resting metabolic demand, integrating consideration of both \({{{{{{\rm{O}}}}}}}_{2}\) availability (in the form of partial \({{{{{{\rm{O}}}}}}}_{2}\) pressure: \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\)) and \(T\) effects (“Methods” section)3,10. The Metabolic Index depends on two key metabolic traits of marine species: the hypoxic tolerance (\({A}_{{{{{{\rm{o}}}}}}}\)) and its temperature sensitivity (\({E}_{{{{{{\rm{o}}}}}}}\)). These traits have been synchronously determined for a diversity of species using published laboratory and field data3. Observed marine species distributions better coincide with the spatial distribution of this index than models based on either \(T\) or \({{{{{{\rm{O}}}}}}}_{2}\) alone3,10. Occurrences of marine species are only found where \(\varPhi\) is above a critical value (\({\varPhi }_{{{{{{\rm{crit}}}}}}}\)), which equals to 1 for resting metabolism and is larger than 1 for active metabolism in support of growth and essential ecological activities3,10. The capacity of the \(\varPhi\) to identify viable habitats of marine species has been successfully employed to explore how variability in the environmental variables can affect the distribution of viable habitats for the past, contemporary, and projected future oceans3,5,6,11. In regions of strong climatic variability of \({{{{{{\rm{O}}}}}}}_{2}\), its contribution to \(\varPhi\) is critical to explain observed variations in population abundance and geographic range12,13, suggesting that \({{{{{{\rm{O}}}}}}}_{2}\) may greatly enhance the predictability of marine habitat shifts.
Global Earth System Models (ESMs) have demonstrated skill in forecasting physical and biogeochemical variables important to marine species on seasonal to decadal time scales14,15. The predictive capabilities of ESMs in combination with the \(\varPhi\) enables the development of forecasts that explicitly consider impacts of \(T\) and \({{{{{{\rm{O}}}}}}}_{2}\) availability on habitat viability. This index in the upper-400 m ocean is projected to reduce by ~20% globally and by ~50% in northern high-latitude regions by the end of this century, due to the combined effects of warming and deoxygenation under a high emissions scenario3. The projected reductions would generate spatial shifts and contractions of habitats and habitable seasons, posing a significant risk for marine ecosystems16,17,18,19. Global and local extinction risks for various marine biota are assessed based on this index and compared across a range of emission scenarios over the next few hundreds of years20. While the long-term trends of projections are worrisome, there is likely important variability on the interannual-to-decadal timescale, and predictability of this index and its mechanistic controls on this time scale from the physical and biogeochemical driver variables are still unknown. This time horizon, however, is critical for marine resource management to reduce impacts, promote resilience, and maximize the value of living marine resources in the face of changing ocean conditions14,21.
We focus on the normalized Metabolic Index (\(\phi=\varPhi /{\varPhi }_{{{{{{\rm{crit}}}}}}}\)) in this work for a consistent critical value of the \(\phi\) (\({\phi }_{{{{{{\rm{crit}}}}}}}\) = 1) for resting and active metabolisms. Like \(\varPhi,\) \(\phi\) is characterized for a species by two measurable traits. In this case, \({A}_{{{{{{\rm{c}}}}}}}\) equals to \({A}_{{{{{{\rm{o}}}}}}}/{\varPhi }_{{{{{{\rm{crit}}}}}}}\), which is the hypoxic tolerance normalized to its minimum value needed to sustain ecological activity. Thus, ecological habitability in a natural environment is defined by \(\phi\) > 1, in the same way that physiological habitability is defined by \(\varPhi\) > 1. The two metabolic traits of marine species exhibit approximately normal distributions with the \({A}_{{{{{{\rm{c}}}}}}}\) and \({E}_{{{{{{\rm{o}}}}}}}\) mainly within the ranges of 0 to 20 atm−1 and −0.2 to 1.0 eV respectively, based on the available trait database (Supplementary Fig. 1).
We examine the predictability of the \(\phi\) on the interannual-to-decadal timescale in two depth habitats (0-200 m, the surface layer or epipelagic zone where most of the visible lights exist, and 200-600 m, the thermocline layer within the mesopelagic zone or twilight zone) of the upper ocean focusing on the 11 North American Large Marine Ecosystems (LMEs; Fig. 1). A decadal prediction system with embedded ocean biogeochemistry entitled the Community Earth System Model - Decadal Prediction Large Ensemble (CESM-DPLE)22 is employed, which includes a Forced Ocean-Sea Ice (FOSI) reconstruction for prediction initialization (“Methods” section). The reconstruction is skillful in representing observed spatial and temporal variability of \(T\) and \({{{{{{\rm{O}}}}}}}_{2}\) in the upper-200 m ocean, although biases exist (Supplementary Fig. 2). Given the temporally and spatially sparse \({{{{{{\rm{O}}}}}}}_{2}\) observations, we evaluate the initialized CESM-DPLE forecasts against the FOSI reconstruction instead of real observations, so the predictability in this work actually stands for potential predictability, an upper limit of predictability we could actually obtain21,23,24. This decadal prediction system has shown significantly enhanced predictive capacity of physical, biogeochemical, and ecological variables in the ocean21,22,24. Our results suggest that the \(\phi\) is potentially predictable in the upper-600 m ocean, showing a substantial improvement in average predictability timescale from 2 to 6 years against a simple persistence forecast (Fig. 1). Except over some high-latitude coastal regions, predictability of the \(\phi\) is mainly attributed to its \({{{{{{\rm{O}}}}}}}_{2}\) component rather than the \(T\) or salinity (\(S\)) component for the ecotype considered. Distinct differences in predictability exist among ecotypes with different \({E}_{{{{{{\rm{o}}}}}}}\) traits in the LMEs, especially along the northeast coast with the average predictability timescales ranging from 2 to 10 years for low-\({E}_{{{{{{\rm{o}}}}}}}\) species.
Results
Interannual viable habitat shifts indicated by ϕ
We focus on viable habitats of marine species with different \({E}_{{{{{{\rm{o}}}}}}}\) traits (across the range of −0.2 to 1.0 eV) and the same \({A}_{{{{{{\rm{c}}}}}}}\) trait (medium value of 10 atm−1), as the \({A}_{{{{{{\rm{c}}}}}}}\) trait works as a linear scaling coefficient for \(\phi\) (its impact on \(\phi\) is homogeneous in temporal and spatial dimensions). We refer to each trait combination (\({A}_{{{{{{\rm{c}}}}}}}\) and \({E}_{{{{{{\rm{o}}}}}}}\)) as an ‘ecotype’, each of which can characterize multiple named and recognized species. We first select three representative ecotypes with the \({E}_{{{{{{\rm{o}}}}}}}\) traits equal to −0.2, 0.4, and 1.0 eV, termed low (e.g., sea squirt, a cosmopolitan tunicate), medium (e.g., common littoral crab), and high (e.g., northern/deep-water shrimp) \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, respectively. These traits correspond to different viable habitat distributions in space, as well as their interannual shifts, integrated with spatial differences and interannual variabilities of the environmental variables (e.g., T and O2; Fig. 2a–f).
The low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes mainly inhabit (\(\phi\)>1) the subtropical regions covering the LMEs of Insular Pacific Hawaiian (IPH), Gulf of Mexico, and Southeast U.S. Shelf (SEUS), instead of the southwest coast off Mexico (which has extremely low \({{{{{{\rm{O}}}}}}}_{2}\)) and northern high-latitude regions (which have generally high \({{{{{{\rm{O}}}}}}}_{2}\) and low \(T\)), e.g., the LMEs of Eastern Bering Sea (EBS), Aleutian Islands, and Labrador-Newfoundland (LN; Fig. 2a). On the interannual timescale, the viable habitats of low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes have large meridional contractions or expansions in the northern high-latitudes and longitudinal shifts along the southwest coast off the North America (Fig. 2b).
In contrast, the medium- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes mainly inhabit the northern high-latitude regions as well as the North Pacific Current region (around 45°N; Fig. 2c–f), as they usually require lower critical \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\) to meet metabolic demand in these relatively colder regions, a pattern of which generally reflects low \(T\) and abundant \({{{{{{\rm{O}}}}}}}_{2}\). Southern regions that have either low \({{{{{{\rm{O}}}}}}}_{2}\) or high \(T\), have limited habitability for these ecotypes, whereas the low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotype has relatively more viable habitats. Thus, the interannual shifts of viable habitats for medium- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes are generally meridional, suggesting the ecotypes would have northward contractions of viable habitats with the southern boundaries retreated northward in response to warming or deoxygenation and opposite southward expansions of viable habitats in response to cooling or oxygenation.
Within all these LMEs, aerobic habitat viability varies with depth and ecotypes. Shoaling or deepening habitats occurs on the interannual timescale, in response to the interannual variations of environmental \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) (Fig. 2g, h). For the northwest coast LMEs, habitat viability of these ecotypes generally decrease as depth increases in the upper-600 m ocean, with the viable habitats extending from surface to deeper depths for higher-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. The interannual shoaling or deepening of viable habitats is due to the vertical lower boundaries, with higher-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes having larger vertical shifts. The southwest coast LMEs show similar vertical profiles of \(\phi\) as the northwest coast LMEs, but with shallower lower boundaries and less difference in interannual viable habitat shifts across different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. The southeast coast LMEs and the IPH LME suggest distinct layer preference for different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, with the upper-200 m layer more habitable for lower-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes and the lower 200-600 m layer for higher-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. Thus, viable habitats in these LMEs would generally shoal for lower-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, deepen for higher-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, and even with both upper and lower boundaries shrink inward for some medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (e.g., in the IPH LME), in response to warming or deoxygenation. Although lower-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes have less temperature sensitivity, their viable habitats still vary strongly with \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\) on the interannual timescale. The northeast coast LMEs are habitable for almost all the \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes considered in the upper 600-m ocean, except those with low \({E}_{{{{{{\rm{o}}}}}}}\) traits in the subsurface (below ~50-m depth). However, their interannual viable habitat shifts can span hundreds of meters at depth, due to nearly vertical profiles of \(\phi\) which are around the critical value of habitat viability (\({\phi }_{{{{{{\rm{crit}}}}}}}\) = 1).
Interannual-to-decadal predictability of ϕ
To assess the interannual-to-decadal predictability of habitat viability indicated by \(\phi\), we present two metrics of prediction skill - the Anomaly Correlation Coefficient (ACC; a measure of degree of associations and the most commonly used metric for forecast accuracy) and Normalized Mean Absolute Error (NMAE; a scale-normalized measure of distance or difference corroborating the ACC), between yearly anomalies of the CESM-DPLE ensemble-mean forecast and the FOSI reconstruction (formally, the predictand; “Methods” section). The resulting skills are then compared against those of a simple persistence forecast using the FOSI reconstruction, which measures the lowest predictability we can obtain (“Methods” section).
We first evaluate the interannual-to-decadal predictability of habitat viability for the medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotype. The retrospective forecasts of \(\phi\) have significantly (at the 95% confidence level; “Methods” section) higher predictability than the simple persistence forecasts in the upper-600 m ocean (Fig. 1), with improvement varying between regions, depths, and lead years (Figs. 3–4 and Supplementary Figs. 3–5). The two prediction skills - ACC and NMAE, suggest consistent spatial and vertical patterns for the persistence and DPLE forecasts, with higher ACCs corresponding to lower NMAEs (Fig. 3). The persistence forecasts only have predictability timescale up to ~2 years in both depth layers of the eleven LMEs, with only a few spots that have longer timescales of predictability (~3 years), e.g., the northern Gulf of Alaska and the Labrador shelf and slope (Fig. 3a, b and Supplementary Fig. 3a, b). The DPLE forecasts provide significantly higher (at the 95% confidence level) prediction skills in almost all the LMEs, depth habitats, and lead years, with the average ACC improvement (\(\Delta {{{{{\rm{ACC}}}}}}\)) ranging from 0.1 in the IPH and SEUS to 0.4 in the Aleutian Islands (Supplementary Fig. 4). The average predictability timescales also increase up to 1 or 2 years in the IPH and SEUS and up to ~6 years in the Aleutian Islands using the DPLE forecasts. Although spatial difference in predictability exists within some LMEs, e.g., the mid-shelf of EBS, northern Gulf of Alaska, southern and northern California Current LME, and Labrador shelf have higher ACCs and longer predictability timescales against the rest regions of the corresponding LME depth habitats (Fig. 3c, d and Supplementary Figs. 3c, d and 4).
Improvement in prediction skill stems mainly from subsurface habitat, consistent with the preponderance of low-frequency variability in the ocean interior. For example, the DPLE forecast has achieved predictability timescale up to 9 years in the Gulf of Alaska at depths of 300-500 m and up to 7 years in the Aleutian Islands at depths of 100-300 m (Fig. 4). In the neighboring EBS LME, the DPLE forecast also has achieved higher prediction skills at depth than the persistence forecast, with the average predictability timescale of ~3 years against 1 year. The northeast coast LMEs not only have improved predictability from the DPLE at depth (below 100-m depth; up to ~4 years of the predictability timescale), but also higher prediction skills and longer timescales of predictability near the surface, e.g., up to ~7 years of predictability timescale in the upper-100 m of the Northeast U.S. Shelf (NEUS), Scotian Shelf, and the LN (mainly the Labrador shelf). The southeast coast LMEs and the IPH LME have substantial layer difference in prediction skills of the DPLE, with low prediction skills in ACC (~0.1) in the upper-200 m layer and much higher ACCs (~0.8) in the deeper 200-600 m layer. The significantly improved prediction skills (e.g., ACC increase by ~0.3, NMAE decrease by ~0.4, up to ~6 years against 1 or 2 years of predictability timescale in the SEUS) are mainly below the 200-m depth (Fig. 4 and Supplementary Fig. 5).
The LMEs of Gulf of California and California Current along the southwest coast suggest significant prediction skill improvement at depth as well (average ACC increase by ~0.4 and ~0.2, respectively), but the California Current LME exhibits a distinct prediction skill pattern from the DPLE forecast. Instead of steadily decreasing with longer lead years, prediction skill initially declines but rebounds at longer (4-10) lead years (Fig. 4e). The reemergent ACCs are located at ~100 m depth at lead year 4 and gradually extend deeper to ~300 m as the lead year increases to 10. The near-surface (<50 m) layer of the IPH LME also suggest similar (but weaker) reemergent ACCs at longer (4-10) lead years. These two LMEs’ reemergence in prediction skills at 4 to 10 lead years is probably associated with the mean gyre-scale circulation in the North Pacific, i.e., mean advection of subsurface water mass anomalies25. Outside of these LMEs, the Labrador Sea (up to 6 years) and some open ocean regions (e.g., up to 9 years in the region south of Aleutian Islands) also exhibit long predictability timescales of the \(\phi\), indicating the advantage of using the DPLE forecasts in these regions as well.
Mechanistic controls on interannual-to-decadal predictability
We investigate which driver variable lends the predictability to \(\phi\) for the medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotype by linearly decomposing \(\phi\) into all possible driver variable components - \({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\) (see “Taylor Linear Decomposition” in “Method” section) and assessing these components’ prediction skills relative to those of the \(\phi\). We find that the \({{{{{{\rm{O}}}}}}}_{2}\) component is the major contributor to the predictability of \(\phi\) in most regions, depths, and lead years (Fig. 5 and Supplementary Fig. 6-7). However, for some LMEs, predictability of the \({{{{{{\rm{O}}}}}}}_{2}\) component is very limited compared to that of \(\phi\), for example, the inner shelves of EBS and the southeast coast LMEs (Gulf of Mexico and SEUS) at shorter (1-2) lead years, and the Labrador shelf at all lead years (Fig. 5a, d and Supplementary Figs. 6a and 7a). Outside of the LMEs, limited predictability of the \({{{{{{\rm{O}}}}}}}_{2}\) component is also found in the upper-200 m layer of the Atlantic subtropical coastal regions (e.g., the Caribbean Sea) and both depth layers of the Labrador Sea (Supplementary Figs. 6a, b and 7a, b). The contribution to the predictability of \(\phi\) from its \(T\) component cannot be neglected in these regions, especially in the northern high-latitude regions along the inner shelf of the EBS and the Labrador shelf where the \(T\) component even becomes the dominant driver variable (Fig. 5b and Supplementary Figs. 6c and 7c). The \(T\) component contributes comparative predictability to \(\phi\) with the \({{{{{{\rm{O}}}}}}}_{2}\) component in some other regions as well, though mainly within the near-surface (<100 m) layer, for example, the northeast (NEUS and Scotian Shelf) and southwest (California Current and Gulf of California) coast LMEs and the IPH LME (Fig. 5b, e). The impact from the \(S\) component, which arises from impact of \(S\) on \({{{{{{\rm{O}}}}}}}_{2}\) solubility (\({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\)) in calculation of \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\), is more limited compared to those of \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) and mainly confined to the Labrador shelf and sea (Fig. 5c, f and Supplementary Figs. 6d and 7d).
A variance budget analysis of \(\phi\) and its components (“Methods” section) further confirms the dominant role of the \({{{{{{\rm{O}}}}}}}_{2}\) component in most regions of the study domain and the nonnegligible impacts from the \(T\) component along the inner shelf of EBS and the northeast coast (Scotian Shelf and LN), where interannual variance of the \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) components is partly contradicted by their negative covariance at both depth habitats (Fig. 6a–e and Supplementary Fig. 8a–e). A further decomposition of the \({{{{{{\rm{O}}}}}}}_{2}\) component into contributions associated with \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\) and apparent oxygen utilization (AOU, “Methods” section) suggests that the AOU is the dominant component in driving the interannual variability of the \({{{{{{\rm{O}}}}}}}_{2}\) component at both depth habitats in the Pacific regions (Fig. 6f–h and Supplementary Fig. 8f–h). While for the northeast coast regions in the Atlantic, interannual variability of the \({{{{{{\rm{O}}}}}}}_{2}\) component is jointly contributed by the \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\) and AOU components as well as their covariance at both depth habitats. Variance of the \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\) component has relatively larger magnitude than that of the AOU component at the upper-200 m layer of the northeast coast regions, which is similar for the further decomposition of the covariance between the \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) components at these two layers (Fig. 6f–j and Supplementary Fig. 8f–j). Thus, the \({{{{{{\rm{O}}}}}}}_{2}\) variability, especially the AOU variability, are critical factors for skillful predictions of \(\phi\) for the medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes in most regions on the interannual-to-decadal timescale; while in those high-latitude coastal regions - the EBS inner shelf and the northeast coast regions (especially Scotian Shelf and LN), the physical drivers (variability of \(T\) and \(T\)- and \(S\)- dependent \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\)) can be equally important and cannot be neglected, especially in the near-surface (<100 m) layer.
Predictability differences across different E o ecotypes
The interannual-to-decadal predictability of aerobic habitat viability may vary largely among different marine ecotypes, due to their difference in the temperature sensitivity of hypoxia tolerance (i.e., \({E}_{{{{{{\rm{o}}}}}}}\)). We calculate prediction skill for ecotypes across a broad range of the \({E}_{{{{{{\rm{o}}}}}}}\) trait (from −0.2 to 1.0 eV with an interval of 0.1 eV) and find that the predictability differences among different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes are most notable in the northeast coast LMEs (Fig. 7 and Supplementary Fig. 9). For example, the upper depth layers of the northeast coast LMEs have predictability timescale up to only 1 or 2 years for ecotypes with the \({E}_{{{{{{\rm{o}}}}}}}\) traits equal to 0.1-0.3, 0.2-0.3, and 0.6-1.0 eV for the NEUS, Scotian Shelf, and LN, respectively. But they have much longer predictability timescales for some higher- or lower- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, e.g., up to 10 years in the NEUS LME for ecotypes with the \({E}_{{{{{{\rm{o}}}}}}}\) equal to 0 eV (Fig. 7a). The deeper depth layer of the northeast coast LMEs maintains the predictability differences similar to those of the upper layer. For example, the deeper layer of the NEUS has up to 4 years of predictability timescale for medium-to-high \({E}_{{{{{{\rm{o}}}}}}}\) (above 0.2 eV) ecotypes but has up to 8-10 years for low-\({E}_{{{{{{\rm{o}}}}}}}\) (below 0.1 eV) ecotypes (Supplementary Fig. 9a).
The upper depth layer of the northwest coast LMEs and the deeper depth layer of the SEUS LME also have notable predictability differences across different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. The former has longer predictability timescale up to ~6 years for medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (0.4, 0-0.4, and 0.1-0.2 eV for the LMEs of EBS, Aleutian Islands and Gulf of Alaska, respectively) and very limited predictability timescale up to only 1-2 years for low- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (e.g., \({E}_{{{{{{\rm{o}}}}}}}\) below 0 and −0.1 eV in the LMEs of EBS and Gulf of Alaska, respectively, and above 0.5, 0.7, 0.8 eV in the LMEs of EBS, Aleutian Islands and Gulf of Alaska, respectively; Fig. 7a). Similarly, the latter has longer predictability timescale up to 6-7 years for medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (0.1–0.6 eV) and limited predictability timescale up to only 3-4 years for low- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (below −0.1 eV and above 0.8 eV; Fig. 7e and Supplementary Fig. 9e). The rest depth habitats of LMEs have relatively little differences (less than 1 or 2 years) in predictability (especially prediction skills) across different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, although significance tests of prediction skills may result in some difference in the predictability timescales (e.g., in the southwest coast LMEs).
The mechanistic controls to the predictability differences across different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes in these regions can be explained by relative contributions of the driver variable (\({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\)) components. Generally, the \({{{{{{\rm{O}}}}}}}_{2}\) component is the major contributor to the predictability of \(\phi\) in most LMEs (especially for the lower depth layer of the Pacific LMEs), but not for all the ecotypes across this \({E}_{{{{{{\rm{o}}}}}}}\) range (Fig. 7 and Supplementary Fig. 9). For example, in the upper-200 m of the northwest coast LMEs, the trait-space predictability differences are mainly contributed by the \({{{{{{\rm{O}}}}}}}_{2}\) component, especially for the medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. The \(T\) component also contributes to the predictability but only for the low- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (e.g., below 0.1 eV and above 0.8 eV in the LME of EBS with limited predictability timescale of 1 or 2 years; Fig. 7a–d). Similarly, at the deeper depth layer of the SEUS, the higher prediction skills (up to 6 years of predictability timescale) for medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (0.2-0.4 eV) are mainly contributed by the \({{{{{{\rm{O}}}}}}}_{2}\) component. The \(T\) component does contribute to the predictability but only for relatively low- (below −0.1 eV) and high- (above 0.8 eV) \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes at short-time (1 or 2) lead years (Supplementary Fig. 9e–h).
In contrast to other regions, in the northeast coast LMEs, the \(T\) component plays a strong role in trait-space predictability differences (Fig. 7a–d and Supplementary Fig. 9a–d). The \({{{{{{\rm{O}}}}}}}_{2}\) component only has comparative predictability (to that of the \(\phi\)) for habitat predictions of some medium-to-high \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, e.g., above 0.2, 0.4, and 0.6 eV in the upper layer of NEUS, Scotian Shelf, and LN, respectively. The \(S\) component also presents as an important contributor to the predictability in these LMEs, but only for some relatively low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, e.g., below 0.2, 0.3, and 0.4 eV in the upper-200 m of NEUS, Scotian Shelf, and LN, respectively. It suggests a substantial impact from the \(T\)- and \(S\)- dependent \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\) on the predictability of habitat viability for these relatively low \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. The dominance of the \({{{{{{\rm{O}}}}}}}_{2}\) or \(T\) component in predictability differences of some southwest and southeast coast LMEs is generally less clear than the northern high-latitude LMEs, due in part to the overall limited predictability of \(\phi\) (e.g., the southeast coast LMEs and the IPH LME at the upper-200 m layer). The \(T\) component, however, does contribute to the predictability in the LMEs of California Current and Gulf of California at the upper-200 m layer, especially for relatively high \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (e.g., above 0.6 eV in the Gulf of California; Fig. 7e–h).
Discussion
Recent advances in modeling and forecasting the earth system (e.g., refs. 14,26,27) enable rapid expansion of marine ecological and habitat forecasts. Such forecasts have typically focused on \({T}\) as the primary if not sole determinant of ecological niche (e.g., ref. 28). This study expands the scope and skill of these tools by highlighting the potential that \({{{{{{\rm{O}}}}}}}_{2}\) can play in deriving forecasts of marine habitats in the upper-600 m of the ocean. A key normalized Metabolic Index (\(\phi\)) enables quantifying habitat constraints for different species arising from their metabolic dependence on \({T}\) and requirements for \({{{{{{\rm{O}}}}}}}_{2}\)3,10. We combine this eco-physiologically mechanistic framework with a full suite of initialized retrospective decadal prediction system embedded with ocean biogeochemistry, to investigate the interannual-to-decadal predictability of aerobic habitat viability (\(\phi\)) for diverse marine ecotypes, its mechanistic linkages with the environmental driver variables, and the predictability differences across different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes.
We select three representative marine ecotypes with low (−0.2 eV; e.g., sea squirt), medium (0.4 eV; e.g., common littoral crab), and high (1.0 eV; e.g., northern/deep-water shrimp) \({E}_{{{{{{\rm{o}}}}}}}\) traits to demonstrate the model performance in capturing the interannual-to-decadal habitat shifts, which exhibit large spatial differences in the habitat viability (\(\phi\)) among different ecotypes, as well as their habitat shifts on this timescale. Due to spatial differences in the environmental driver variables (e.g., \({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\)), the high-latitude coastal regions have high habitat viability with deep vertical extension for the medium- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes but limited viability for the low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes (Fig. 2), while the subtropical regions have the opposite (e.g., the Gulf of California) or layer-specified (e.g., the Gulf of Mexico and SEUS) habitat viability for these ecotypes. In response to interannual variations of environmental variables, the viable habitats of medium- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes mainly have their southern boundaries shifted meridionally (e.g., retreat northward in response to warming or deoxygenation), while those of the low-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes have both northern and southern or southeastern boundaries shifted not only meridionally but zonally (e.g., along the west coast of California). The habitat contraction or expansion trends are consistent with the observed relocation of various economically important species along both west and east coasts of North America. For example, the medium-\({E}_{{{{{{\rm{o}}}}}}}\) species of American lobster and summer flounder are both shifting northward along the NEUS during recent years (e.g., refs. 29,30,31).
Combined with this mechanistic framework, the decadal prediction system is used to evaluate the interannual-to-decadal predictability of habitat viability at both depth habitats and vertically within the LMEs, providing significantly higher (at the 95% confidence level) prediction skills and longer predictability timescales than the simple persistence forecast. For marine ecotypes with the medium-\({E}_{{{{{{\rm{o}}}}}}}\) trait, prediction skills of the DPLE are spatially inconsistent, with higher ACCs over the mid-shelf of EBS, northern Gulf of Alaska, southern and northern California Current LME, and the Labrador shelf (Fig. 3). Vertically, these improved prediction skills are mainly in the subsurface, e.g., below 100-m depth for the west (northwest and southwest) coast LMEs and 200-m depth for the southeast coast LMEs (Fig. 4). For habitat viability in the northeast coast LMEs, the DPLE has achieved both higher ACCs at depth and in the near-surface (<100 m) layer. For the LMEs of California Current and IPH, the decadal prediction system not only has achieved higher prediction skills at shorter (1-3) lead years but realized reemergent skillful ACCs at lead years from 4 to 10. These reemergent ACCs may be associated with the mean advection of subsurface water mass anomalies in the North Pacific at comparable timescales and warrants further dynamical investigation.
The decadal prediction system combined with the mechanistic framework also enables investigating the relative mechanistic controls to the predictability of habitat viability from the environmental driver variables (\({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\)). The Taylor linear decomposition and variance budget analysis suggest that the \({{{{{{\rm{O}}}}}}}_{2}\) component of \(\phi\) is the primary control on the predictability in most regions (especially the Pacific) for the medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, which demonstrates the importance of reliable \({{{{{{\rm{O}}}}}}}_{2}\) forecasts and measurements for skillful forecasts of habitat viability at a range of depths (e.g., ref. 32). However, contribution of the \(T\) component to the predictability can also be important or even determinant over the high-latitude regions, e.g., the EBS inner shelf and the Labrador LME. This is probably because of relaxation of the \({{{{{{\rm{O}}}}}}}_{2}\) constraint on habitat viability due to strong ventilation processes in these high-latitude regions. While the physical variables such as \(T\) and \(T\)- and \(S\)- dependent \({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}\) exert larger impacts on habitat viability in these regions, these physical impacts are mainly in the near-surface (<100 m) layer, i.e., in the northeast coast LMEs (Fig. 5). This work also identifies the importance of AOU for interannual-to-decadal habitat predictions. Improved simulation and representation of the processes responsible for AOU variability - ventilation and organismal \({{{{{{\rm{O}}}}}}}_{2}\) production and consumption, would also further improve our habitat forecasts. The former could be enhanced with higher resolutions, while the latter remains unresolved and requires more physiological, ecological, and modeling efforts.
This study also provides a trait-space view of predictability differences across marine ecotypes with different temperature sensitivities of hypoxia vulnerability - the \({E}_{{{{{{\rm{o}}}}}}}\), which can be explained by the relative dominance of the environmental driver variables. The trait-space predictability differences are prominent in the upper-200 m of the northwest coast LMEs, both depth layers of the northeast coast LMEs, and lower 200-600 m of the southeast coast LMEs (Fig. 7 and Supplementary Fig. 9). It suggests predictability of habitat viability in these regions may need case-by-case investigation for certain species and the DPLE system may have extremely high performance in decadal habitat predictions of some species, e.g., low \({E}_{{{{{{\rm{o}}}}}}}\) (0 eV) ecotypes in the NEUS LME. For most LMEs that have notable predictability differences (the upper layer of the northwest coast LMEs and the deeper layer of the SEUS LME), the \({{{{{{\rm{O}}}}}}}_{2}\) component is the major contributor for medium-\({E}_{{{{{{\rm{o}}}}}}}\) ecotypes, which suggests much longer predictability timescales than those of the \(T\) component for low- and high- \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes. While predictability differences of the northeast coast LMEs are mainly contributed by the \(T\) component (especially for low-to-medium \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes with up to 10 years of predictability timescale) and impacts from the \({{{{{{\rm{O}}}}}}}_{2}\) component exist only for certain \({E}_{{{{{{\rm{o}}}}}}}\) range (medium-to-high) of ecotypes. We focus on habitat shifts and viability predictions for different \({E}_{{{{{{\rm{o}}}}}}}\) ecotypes with the same \({A}_{{{{{{\rm{c}}}}}}}\) (10 atm-1) trait. Due to the linear scaling effect of the \({A}_{{{{{{\rm{c}}}}}}}\) trait, these results are applicable to any other ecotypes with the same \({E}_{{{{{{\rm{o}}}}}}}\) but different \({A}_{{{{{{\rm{c}}}}}}}\) traits. For example, the interannual habitat shifts are similar for ecotypes with the same \({E}_{{{{{{\rm{o}}}}}}}\) but smaller (larger) \({A}_{{{{{{\rm{c}}}}}}}\) traits, only having the habitat boundaries moving northward (southward) or towards high (low) \(\phi\) regions. The prediction skill assessments in ACC are also the same between the same-\({E}_{{{{{{\rm{o}}}}}}}\) different-\({A}_{{{{{{\rm{c}}}}}}}\) ecotypes, as the linear scaling coefficient \({A}_{{{{{{\rm{c}}}}}}}\) is self-canceled.
Our results demonstrate the potential of using an initialized ESM to retrospectively predict habitat viability for diverse marine ecotypes on the interannual-to-decadal timescale. With skillful interannual-to-decadal predictions of the \(\phi\), spatial extents of viable habitats and the corresponding areas and volumes for diverse marine ecotypes can be estimated. The generated habitat forecast products can be fully employed for living marine resource management and decision making in response to the changing ocean conditions (e.g., ref. 33). For example, northward habitat contractions can make it infeasible to fish from certain ports, while southward habitat expansion opens ports that would need to be prepared for processing the landings. Changes in depth distribution can also necessitate changing the type of fishing gear used. Additionally, unequal spatial and vertical habitat shifts of species due to their different metabolic traits (e.g., temperature sensitivity) may lead to substantial changes in prey-predator dynamics and ecosystem structure (e.g., overlap of viable habitats of southern and northern silver hake with juveniles as prey and adults as predator, ref. 34), which may affect resource availability to the fishery and require adaption by stakeholders33.
Although this study presents a promising result, it is important to note that the coarse spatial resolution (~1°) may not resolve well the fine-scale oceanic processes that are important for physical and biogeochemical variations. Higher-resolution global forecast models (e.g., the developing high-resolution CESM decadal forecasts at ~0.1°) or regional models dynamically downscaled from the global initialized ESM forecasts can be used to further improve the habitat viability forecasts, based on the framework of \(\phi\). However, reliable observations or observation-based products of \({{{{{{\rm{O}}}}}}}_{2}\) with spatio-temporal variability are urgently needed to be able to truly verify the realized prediction skills and move from potential predictability to verified prediction skill. These observational products exist for \(T\), \(S\), and chlorophyll, but are just beginning to be developed for \({{{{{{\rm{O}}}}}}}_{2}\) (e.g., the GOBAI-O2 product, ref. 35) and other biogeochemical variables.
Methods
The Metabolic Index
The Metabolic Index (\(\varPhi\)) is defined as the ratio of dissolved oxygen (\({{{{{{\rm{O}}}}}}}_{2}\)) supply to an organism’s resting metabolic demand for evaluation of the aerobic energy balance of an organism, as expressed by this equation:
where \({A}_{{{{{{\rm{o}}}}}}}\) is the ratio of efficacy for \({{{{{{\rm{O}}}}}}}_{2}\) supply (\({\alpha }_{{{{{{\rm{S}}}}}}}\)) and minimum metabolic rate (\({{{{{{\rm{O}}}}}}}_{2}\) demand; \({\alpha }_{{{{{{\rm{D}}}}}}}\)) at reference temperature (\({{{{T}}}}_{{{{{{\rm{ref}}}}}}}\); \({A}_{{{{{{\rm{o}}}}}}}=\frac{{\alpha }_{{{{{{\rm{S}}}}}}}}{{\alpha }_{{{{{{\rm{D}}}}}}}}\)), \({E}_{{{{{{\rm{o}}}}}}}\) is the difference between temperature variation in the metabolic rate (\({E}_{{{{{{\rm{d}}}}}}}\)) and that in the oxygen supply (\({E}_{{{{{{\rm{s}}}}}}}\); \({E}_{{{{{{\rm{o}}}}}}}\) = \({E}_{{{{{{\rm{d}}}}}}}\)-\({E}_{{{{{{\rm{s}}}}}}}\))10, \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\) is the partial pressure of \({{{{{{\rm{O}}}}}}}_{2}\) calculated using \({{{{{{\rm{O}}}}}}}_{2}\) and temperature (\(T\)) and salinity (\(S\)) dependent oxygen solubility36, and \({k}_{{{{{{\rm{B}}}}}}}\) is the Boltzmann constant37. To account for the variation of \({E}_{{{{{{\rm{o}}}}}}}\) with \({{{T}}}\), the \({E}_{{{{{{\rm{o}}}}}}}\) value is modified by a linear term: \({E}_{{{{{{\rm{o}}}}}}}\to {E}_{{{{{{\rm{o}}}}}}}+\frac{d{E}_{{{{{{\rm{o}}}}}}}}{{dT}}(T-{T}_{{{{{{\rm{ref}}}}}}})\), with the \(\frac{d{E}_{{{{{{\rm{o}}}}}}}}{{dT}}\) equal to an average value of 0.022 eV/K for ecotypes considered in this study and the \({T}_{{{{{{\rm{ref}}}}}}}\) equal to 288.15 K10. Spatial distributions of this index compared with observed occurrence of various species with different \({E}_{{{{{{\rm{o}}}}}}}\) traits can be found in the Supplementary Information of ref. 10.
We focus on the normalized Metabolic Index (\(\phi\); \(\phi=\varPhi /{\varPhi }_{{{{{{\rm{crit}}}}}}}\)), with the \({A}_{{{{{{\rm{o}}}}}}}\) normalized to the \({\varPhi }_{{{{{{\rm{crit}}}}}}}\) at \({T}_{{{{{{\rm{ref}}}}}}}\). If \(\phi\) falls below 1, the environment can no longer support the aerobic demands of species’ energetic requirements. Conversely, values of \(\phi\) above 1 permit critical activities such as feeding, defense, growth, and reproduction.
To calculate depth and LME averages of the \(\phi\), we first calculate the \(\phi\) at each individual grid cell, and then perform depth-averaging at each layer and area-averaging within each LME, with depth- and area-weighted functions applied. To evaluate habitat viability (indicated by the \(\phi\)) on the interannual timescale, we first calculate the \(\phi\) from monthly model outputs of \({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\), and then calculate yearly averages of the \(\phi\) by averaging from January to December. To investigate interannual viable habitat shifts, we first calculate interannual standard deviation of the \(\phi\) (\({\sigma }_{\phi }\)) over the period of 1954–2017, and then identify the viable habitats (\(\phi\) > 1) in space and vertically within the LMEs, with the \(\phi\) varying between \(\phi\)±3\({\sigma }_{\phi }\). The corresponding spatial (vertical) differences in the viable habitats suggest habitat contraction or expansion (shoaling or deepening) on the interannual timescale.
The Decadal Prediction System
The CESM-DPLE is a suite of initialized retrospective Earth System Model simulations with embedded ocean biogeochemistry22. Horizontal resolution of the CESM-DPLE ocean component is nominally 1°\(\times\)1° with 60 vertical levels. The ocean and sea ice model components in the CESM-DPLE were re-initialized from a forced ocean-sea ice reconstruction (referred to as the FOSI reconstruction) each year on November 1st from 1954 to 201722. For each initialization date, an ensemble of 40 forecast members was created by creating Gaussian perturbations to the initial atmospheric temperature field (order 10−14 K) at each grid cell integrating forward for ~10 years, and model components developed as a result of the spread in the atmospheric state22.
The Biogeochemical Elemental Cycling (BEC) model is used to simulate biogeochemistry components in CESM-DPLE38,39. Note that the ocean biogeochemistry is diagnostic that has no feedback onto the simulated physical climate22. Historical radiative forcing with volcanic aerosols is included in CESM-DPLE through 2005, and projected radiative forcing (including greenhouse and short-lived gases and aerosols) is added from 2006 onward.
Drift adjustment is performed with CESM-DPLE, to correct for model drift caused by full-field initialization. The correction procedures are the same as described in ref. 22: the lead-time dependent model climatology is computed as the mean drift across all 40 ensemble members and 64 initialization dates between 1954 and 2017; and this drift was then subtracted at each grid cell from all forecasts to generate anomalies.
The FOSI reconstruction is a hindcast simulation from 1948 to 2017, with active ocean (physics and biogeochemistry) and sea ice model components from CESM that has identical grids as the fully coupled CESM-DPLE. It is forced at the surface by a modified version of the Coordinated Ocean-Ice Reference Experiment (CORE) with interannual forcing40,41. The FOSI reconstruction reproduces some key aspects of observed ocean (e.g., ocean temperature and dissolved oxygen as shown in Supplementary Fig. 2) and sea ice variability quite well, although there is no direct assimilation of either ocean or sea ice observations42,43,44. The FOSI reconstruction is used to evaluate the prediction skill of the simple persistence forecast, as a reference compared with the CESM-DPLE’s prediction skill.
Prediction skill assessments
Predictability is assessed by calculating the Anomaly Correlation Coefficient (ACC) and the Normalized Mean Absolute Error (NMAE) between yearly anomalies from the CESM-DPLE ensemble-mean forecast and the FOSI reconstruction, both of which are considered as functions of lead time in years:
where \({F}^{{\prime} }\) and \({R}^{{\prime} }\) represent yearly anomalies of the DPLE ensemble-mean forecast and the FOSI reconstruction respectively, \(\tau\) is the lead time in years ranging from 1 to 10, \(N\) is the length of the time series, and \({\sigma }_{{R}^{{\prime} }}\) is the interannual standard deviation of the \({R}^{{\prime} }\) time series. The NMAE is considered as an accurate assessment of bias when quantifying the accuracy of the forecasts21,45. Statistical significance of ACC is tested nonzero at the 95% confidence level via a Student’s t-test, considering the effective degree of freedom to account for autocorrelation of the two timeseries being correlated46. Higher values of ACC and lower values of NMAE suggest better prediction skills from the DPLE forecast.
To show the utility of the initialized forecasting system (DPLE) over a simple, low-cost forecasting method, we also evaluate these prediction skills of a persistence forecast using the FOSI reconstruction (hereafter named reconstruction persistence), which predicts the yearly anomalies of the FOSI reconstruction using predictors of itself but with \(\tau\) lead years. The ACC and NMAE equations are the same as Eqs. (2) and (3), except replacing the \({{F}^{{\prime} }}_{i}(\tau )\) with the \({{R}^{{\prime} }}_{i}\). To assess if the prediction skill of the DPLE is better over that of the reconstruction persistence, we calculate statistical significance of the difference between these two ACCs (DPLE minus persistence), which is also tested larger than 0 at the 95% confidence level following ref. 47,48, taking into account the fact that the two correlations are based on a common FOSI reconstruction time series (\({R}^{{\prime} }\))49.
The predictability timescale is quantified as the maximum lead time when the corresponding ACCs are continuously significant (at the 95% confidence level) starting from lead year 1. For example, if ACC is statistically significant at lead year 1 through 5 and at lead year 8 through 9, the predictability timescale is 5 years instead of 9 years.
Taylor linear decomposition and variance budget
To understand the mechanistic controls towards the predictability, a first-order Taylor-series decomposition method is applied to decompose the normalized Metabolic Index (\(\phi\)) into contributions from the physical (\(T\), \(S\)) and biogeochemical (\({{{{{{\rm{O}}}}}}}_{2}\)) driver variables. For a certain marine species with known metabolic traits (\({A}_{{{{{{\rm{c}}}}}}}\) and \({E}_{{{{{{\rm{o}}}}}}}\)), \(\phi\) is a function of \({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\) (Eq. 1; \({p}_{{{{{{{\rm{O}}}}}}}_{2}}\) is a function of \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\)- and \(S\)- dependent oxygen solubility), and its first-order Taylor-series decomposition is:
where \({x}_{i}\) represents one of the three driver variables (\({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\); \(n\) = 3), \(\bar{{x}_{i}}\) is the time-mean value of the variable \({x}_{i}\), and \(\bar{\phi }\) is the mean value of \(\phi\) calculated using all the time-mean driver variables. \(\frac{\partial \phi }{\partial {x}_{i}}\) is the first-order partial derivative of \(\phi\) with respect to the driver variable \({x}_{i}\). The residual term suggests higher-order terms that could be omitted. In other words, we decompose variations of \(\phi\) into a constant time-mean term (\(\bar{\phi }\)) and three driver variable components:
where \(t\) represents time. The \({\phi }_{{{{{\rm{O}}}}}_{2}}\left(t\right)\) could be further decomposed into an oxygen solubility component (\({\phi }_{{{{{\rm{O}}}}}_{2}^{{{{{\rm{sol}}}}}}}\)) and an apparent oxygen utilization (\({\phi }_{{{{{\rm{AOU}}}}}}\)) component, depending on their linear relationship: \({{{{{{\rm{O}}}}}}}_{2}\) is the difference between \({{{{{\rm{O}}}}}}_{2}^{{{{{\rm{sol}}}}}}\) and \({{{{{\rm{AOU}}}}}}\) (e.g., ref. 50; \({{{{{\rm{O}}}}}}_{2}={{{{\rm{O}}}}}_{2}^{{{{{\rm{sol}}}}}}-{{{{{\rm{AOU}}}}}}\), where \({{{{{\rm{O}}}}}}_{2}^{{{{{\rm{sol}}}}}}\) nonlinearly depends on \(T\) and \(S\), and \({{{{{\rm{AOU}}}}}}\) reflects the balance between respiration-driven oxygen consumption and ocean ventilation32,36,51):
where \({\phi }_{{{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}}=\frac{\partial \phi }{\partial {{{{{{\rm{O}}}}}}}_{2}}\cdot({{{{{{\rm{O}}}}}}}_{2}^{{{{{{\rm{sol}}}}}}}-{\bar{{{{{\rm{O}}}}}}}_{2}^{{{{{\rm{sol}}}}}})\), \({\phi }_{{{{{{\rm{AOU}}}}}}}=-\frac{\partial \phi }{\partial {{{{{{\rm{O}}}}}}}_{2}}\cdot ({{{{{\rm{AOU}}}}}}-\overline{{{{{\rm{AOU}}}}}})\), as \(\frac{\partial \phi }{\partial {{{{{\rm{O}}}}}}_{2}}=\frac{\partial \phi }{\partial {{{{{\rm{O}}}}}}_{2}^{{{{{\rm{sol}}}}}}}=-\frac{\partial \phi }{\partial {{{{{\rm{AOU}}}}}}}\). We perform the same prediction skill assessments to evaluate the predictability of each physical and biogeochemical driver component at each depth habitat and LME. These driver components are calculated using the DPLE and considered as forecasts of the original, unchanged FOSI reconstruction. To better interpret the predictability of these driver components of \(\phi\), we perform a variance budget analysis by decomposing the variance of \(\phi\) (\({\sigma }_{\phi }^{2}\)) into contributions from three driver components and their covariances, based on Eqs. (5) and (6):
where \({\phi }_{i}\) is the driver component with \(i\) representing one of the driver variables (\({{{{{{\rm{O}}}}}}}_{2}\), \(T\), and \(S\); \(n\) = 3), and \({{{{{{\rm{Cov}}}}}}}(\phi _{j},{\phi }_{i})\) represents the covariance between any two of the different driver components. The variance decomposition of \({\phi }_{{{{{{{\rm{O}}}}}}}_{2}}\) also follows this rule, which is:
Similarly, the covariance term between \({{{{{{\rm{O}}}}}}}_{2}\) and \(T\) components can also be decomposed as:
Data availability
Output from the Community Earth System Model Decadal Prediction Large Ensemble (CESM-DPLE) and CESM model reconstruction can be downloaded at [https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.CESM1-CAM5-DP.html]. The metabolic traits data can be downloaded online at the supplementary table of ref. 10 [https://www.nature.com/articles/s41586-020-2721-y]. The GOBAI-O2 observation product35 can be downloaded at [https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.nodc:0259304]. The EN4 temperature observation product52 can be downloaded at [https://www.metoffice.gov.uk/hadobs/en4]. Source data of figures are provided with this paper and shared on Figshare53.
Code availability
The code used to calculate the normalized Metabolic Index, assess predictability, perform linear decomposition, and create all the figures are available on GitHub [https://github.com/zchenocean/viable-habitat-prediction].
References
Davis, J. C. Minimal dissolved oxygen requirements of aquatic life with emphasis on Canadian species: a review. J. Fish. Board Can. 32, 2295–2332 (1975).
Pörtner, H. O. & Peck, M. A. Climate change effects on fishes and fisheries: towards a cause‐and‐effect understanding. J. Fish. Biol. 77, 1745–1779 (2010).
Deutsch, C., Ferrel, A., Seibel, B., Pörtner, H. O. & Huey, R. B. Climate change tightens a metabolic constraint on marine habitats. Science 348, 1132–1135 (2015).
Claireaux, G. & Chabot, D. Responses by fishes to environmental hypoxia: integration through Fry’s concept of aerobic metabolic scope. J. Fish. Biol. 88, 232–251 (2016).
Long, M. C., Ito, T. & Deutsch C. Ocean Deoxygenation: Everyone’s Problem - Causes, Impacts, Consequences And Solutions (eds Laffoley, D. & Baxter, J. M.) (Gland, 2019).
Duncan, M. I., James, N. C., Potts, W. M. & Bates, A. E. Different drivers, common mechanism; the distribution of a reef fish is restricted by local-scale oxygen and temperature constraints on aerobic metabolism. Conserv. Physiol. 8, coaa090 (2020).
Keeling, R. F., Körtzinger, A. & Gruber, N. Ocean deoxygenation in a warming world. Annu. Rev. Mar. Sci. 2, 199–229 (2010).
Breitburg, D. et al. Declining oxygen in the global ocean and coastal waters. Science 359, eaam7240 (2018).
Clarke, A. & Fraser, K. P. P. Why does metabolism scale with temperature? Funct. Ecol. 18, 243–251 (2004).
Deutsch, C., Penn, J. L. & Seibel, B. Metabolic trait diversity shapes marine biogeography. Nature 585, 557–562 (2020).
Penn, J. L., Deutsch, C., Payne, J. L. & Sperling, E. A. Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction. Science 362, eaat1327 (2018).
Howard, E. M. et al. Climate-driven aerobic habitat loss in the California Current System. Sci. Adv. 6, eaay3188 (2020).
Franco, A. C. et al. Impact of warming and deoxygenation on the habitat distribution of Pacific halibut in the Northeast Pacific. Fish. Oceanogr. 31, 601–614 (2022).
Tommasi, D. et al. Managing living marine resources in a dynamic environment: the role of seasonal to decadal climate forecasts. Prog. Oceanogr. 152, 15–49 (2017).
Park, J. Y., Stock, C. A., Dunne, J. P., Yang, X. & Rosati, A. Seasonal to multiannual marine ecosystem prediction with a global Earth system model. Science 365, 284–288 (2019).
Prince, E. D. & Goodyear, C. P. Hypoxia‐based habitat compression of tropical pelagic fishes. Fish. Oceanogr. 15, 451–464 (2006).
Gilly, W. F., Beman, J. M., Litvin, S. Y. & Robison, B. H. Oceanographic and biological effects of shoaling of the oxygen minimum zone. Annu. Rev. Mar. Sci. 5, 393–420 (2013).
Mislan, K. A. S., Deutsch, C. A., Brill, R. W., Dunne, J. P. & Sarmiento, J. L. Projections of climate‐driven changes in tuna vertical habitat based on species‐specific differences in blood oxygen affinity. Glob. Change Biol. 23, 4019–4028 (2017).
Whitney, F. A., Freeland, H. J. & Robert, M. Persistently declining oxygen levels in the interior waters of the eastern subarctic Pacific. Prog. Oceanogr. 75, 179–199 (2007).
Penn, J. L. & Deutsch, C. Avoiding ocean mass extinction from climate warming. Science 376, 524–526 (2022).
Brady, R. X., Lovenduski, N. S., Yeager, S. G., Long, M. C. & Lindsay, K. Skillful multiyear predictions of ocean acidification in the California Current System. Nat. Commun. 11, 1–9 (2020).
Yeager, S. G. et al. Predicting near-term changes in the earth system: a large ensemble of initialized decadal prediction simulations using the community earth system model. Bull. Am. Meteorol. Soc. 99, 1867–1886 (2018).
Meehl, G. A. et al. Decadal climate prediction: an update from the trenches. Bull. Am. Meteorol. Soc. 95, 243–267 (2014).
Krumhardt, K. M. et al. Potential predictability of net primary production in the ocean. Glob. Biogeochem. Cycles 34, e2020GB006531 (2020).
Pozo Buil, M. & Di Lorenzo, E. Decadal dynamics and predictability of oxygen and subsurface tracers in the California Current System. Geophys. Res. Lett. 44, 4204–4213 (2017).
Payne, M. R. et al. Lessons from the first generation of marine ecological forecast products. Front. Mar. Sci. 4, 289 (2017).
O’Kane, T. J. et al. Recent applications and potential of near-term (interannual to decadal) climate predictions. Front. Clim. 5, 1121626 (2023).
Payne, M. R. et al. Skilful decadal-scale prediction of fish habitat and distribution shifts. Nat. Commun. 13, 2660 (2022).
Pinsky, M. L., Worm, B., Fogarty, M. J., Sarmiento, J. L. & Levin, S. A. Marine taxa track local climate velocities. Science 341, 1239–1242 (2013).
Dubik, B. A. et al. Governing fisheries in the face of change: Social responses to long-term geographic shifts in a US fishery. Mar. Policy 99, 243–251 (2019).
Fredston, A. et al. Range edges of North American marine species are tracking temperature over decades. Glob. Change Biol. 27, 3145–3156 (2021).
Frölicher, T. L., Ramseyer, L., Raible, C. C., Rodgers, K. B. & Dunne, J. Potential predictability of marine ecosystem drivers. Biogeosciences 17, 2061–2083 (2020).
Liu, O. R. et al. Species redistribution creates unequal outcomes for multispecies fisheries under projected climate change. Sci. Adv. 9, eadg5468 (2023).
Nye, J. A., Joyce, T. M., Kwon, Y. O. & Link, J. S. Silver hake tracks changes in Northwest Atlantic circulation. Nat. Commun. 2, 412 (2011).
Sharp, J. D. et al. GOBAI-O2: A Global Gridded Monthly Dataset of Ocean Interior Dissolved Oxygen Concentrations Based on Shipboard and Autonomous Observations (NCEI Accession 0259304) (NOAA National Centers for Environmental Information, 2022).
Garcia, H. E. & Gordon, L. I. Oxygen solubility in seawater: better fitting equations. Limnol. Oceanogr. 37, 1307–1312 (1992).
Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).
Lindsay, K. et al. Preindustrial-control and twentieth-century carbon cycle experiments with the Earth System Model CESM1 (BGC). J. Clim. 27, 8981–9005 (2014).
Moore, J. K., Lindsay, K., Doney, S. C., Long, M. C. & Misumi, K. Marine ecosystem dynamics and biogeochemical cycling in the Community Earth System Model [CESM1(BGC)]: Comparison of the 1990s with the 2090s under the RCP4.5 and RCP8.5 scenarios. J. Clim. 26, 9291–9312 (2013).
Griffies, S. M. et al. Coordinated ocean-ice reference experiments (COREs). Ocean Model. 26, 1–46 (2009).
Large, W. G. & Yeager, S. G. The global climatology of an interannually varying air-sea flux data set. Clim. Dyn. 33, 341–364 (2009).
Danabasoglu, G. et al. North Atlantic simulations in Coordinated Ocean-ice Reference Experiments phase II (CORE-II). Part II: inter-annual to decadal variability. Ocean Model. 97, 65–90 (2016).
Yeager, S. & Danabasoglu, G. The origins of late-twentieth-century variations in the large-scale North Atlantic circulation. J. Clim. 27, 3222–3247 (2014).
Yeager, S. G., Karspeck, A. R. & Danabasoglu, G. Predicted slowdown in the rate of Atlantic sea ice loss. Geophys. Res. Lett. 42, 10–704 (2015).
Willmott, C. J. & Matsuura, K. Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Clim. Res. 30, 79–82 (2005).
Bretherton, C. S., Widmann, M., Dymnikov, V. P., Wallace, J. M. & Bladé, I. The effective number of spatial degrees of freedom of a time-varying field. J. Clim. 12, 1990–2009 (1999).
Steiger, J. H. Tests for comparing elements of a correlation matrix. Psychol. Bull. 87, 245–251 (1980).
Chen, Z. et al. Seasonal prediction of bottom temperature on the northeast US continental shelf. J. Geophys. Res. Oceans 126, e2021JC017187 (2021).
Howell, D. C. Statistical Methods for Psychology (Nelson Education, 2009)
Frölicher, T. L., Joos, F., Plattner, G. K., Steinacher, M. & Doney, S. C. Natural variability and anthropogenic trends in oceanic oxygen in a coupled carbon cycle-climate model ensemble. Glob. Biogeochem. Cycles 23, GB1003 (2009).
Garcia, H. E., Locarnini, R. A., Boyer, T. P. & Antonov, J. I. World Ocean Atlas 2005 (ed. Levitus, S.) Vol. 3 (US Government Printing Office, 2006)
Good, S. A., Martin, M. J. & Rayner, N. A. EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates. J. Geophys. Res.: Oceans 118, 6704–6716 (2013).
Chen, Z. et al. Source data for figures in skillful multiyear prediction of marine habitat shifts jointly constrained by ocean temperature and dissolved oxygen. Figshare https://doi.org/10.6084/m9.figshare.24863091 (2023).
Amante, C. & Eakins, B. W. ETOPO1 arc-minute global relief model: procedures, data sources and analysis. NOAA National Geophysical Data Center. https://doi.org/10.7289/V5C8276M (2009).
Acknowledgements
This work was supported by NOAA’s Climate Program Office’s Modeling, Analysis, Predictions, and Projections (MAPP) program NA20OAR4310437 (S.S.) and NOAA’s funding project NA18NOS4780167 (C.D.). We acknowledge our participation in MAPP’s Marine Ecosystem Task Force. This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. The authors acknowledge helpful discussions and suggestions from Drs. J. Huang (WHOI) and A. C. Ross (GFDL). The authors also thank Dr. K. Krumhardt (NCAR) for her help on the CESM-DPLE forecasts. We also thank University of Connecticut Scholarship Facilitation Fund and Department of Marine Sciences for supporting the publication cost.
Author information
Authors and Affiliations
Contributions
Z.C., S.S., and M.L. assembled and analyzed the data. Z.C. wrote the manuscript. S.S. and M.L. supervised the study. Z.C., S.S., M.L., C.P., C.S., and C.D. interpreted the results and clarified the implications.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Nature Communications thanks Mark Payne and the other, anonymous, reviewers for their contribution to the peer review of this work. A peer review file is available.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Chen, Z., Siedlecki, S., Long, M. et al. Skillful multiyear prediction of marine habitat shifts jointly constrained by ocean temperature and dissolved oxygen. Nat Commun 15, 900 (2024). https://doi.org/10.1038/s41467-024-45016-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467-024-45016-5
This article is cited by
-
Opening the door to multi-year marine habitat forecasts
Nature Communications (2024)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.