Global change is fundamentally restructuring marine ecosystems, shifting distributions, phenologies, and interactions among species. Temperature (T), oxygen (O₂), and carbonate chemistry (e.g., pH, Ω A r ) naturally constrain available habitat for marine calcifiers and aerobic animals, but as ocean waters warm, become less oxygenated, and more acidified, these changes are driving habitats beyond the envelope of natural variability, resulting in major changes to species distribution and abundance, and raising the potential for major ecosystem disruptions (Howard et al., 2020; Pinsky et al., 2020). While shifts in species abundance or geographic range can be detected in historical data, local human impacts from nitrogen pollution and coastal eutrophication confound the attribution of biological changes to long-term climate trends. Effective coastal ecosystem management in the face of global change requires the means to both: (1) quantify these fundamental changes to species habitats and (2) disentangle the relative roles of climate change, natural and climatic variability, and local anthropogenic pressures in shaping those habitats. Ocean numerical models are routinely used to project the effects of climate change on shifting habitats and species distributions (Penn and Deutsch, 2022), but few coastal numerical modeling studies have investigated the potential for local coastal eutrophication to constraint marine calcifier and aerobic habitat (Bednaršek et al., 2020). As 80% of global wastewater receives no treatment before discharging to coastal waters (Rangel-Buitrago et al., 2024), such studies can help to understand whether local management of coastal eutrophication could meaningfully increase resilience of ecosystems to climate change.

In semi-enclosed seas and estuaries, the effects of eutrophication on increased primary productivity, enhanced remineralization rates, subsurface O₂ depletion, and acidification are commonly observed within the 100-m isobath (Rabalais, 2005; Rabalais et al., 2009). However, recent work has revealed that such changes can be meaningful even in upwelling-dominated coastal environments, countering the tenet that low O₂ and Ω A r that occur along eastern boundary upwelling systems are only naturally induced, without direct anthropogenic influence (Fennel and Testa, 2019). The Southern California Bight (SCB) is an open embayment situated within an eastern boundary upwelling system, located between the Baja California Peninsula and Point Conception. This region hosts a variety of marine communities and hotspots of biodiversity. The complex seabed topography and presence of islands contributes to complex circulation features. While large-scale circulation features are dominated by the California Current, submesoscale eddies accumulate and redistribute material from the coast (Dong et al., 2009). The SCB receives export from a human population of 22 million, which rivals natural upwelling in magnitude, roughly doubling available nitrogen (Howard et al., 2014). These inputs include point and non-point source discharges to the ocean from 19 ocean outfalls and 75 rivers, which release, on average, 8 million m³ d^-1 of nutrient-enriched water to the ocean (Sutula et al., 2021). The ocean outfalls all have wastewater treatment, but the majority have no nitrogen removal (Sutula et al., 2021). The effects of these inputs are increasing primary production and subsurface remineralization rates along the coast, with corresponding subsurface reductions in O₂ and aragonite saturation state ( Ω A r ) that rival or exceed that of global open-ocean O₂ loss and acidification since the pre-industrial period (Kessouri et al., 2021b).

While Kessouri et al. (2021b) quantified the change in seawater chemistry from anthropogenic nutrient inputs in the Bight, they did not document the potential for biological effects, a fundamental science gap that motivates coastal water quality managers. In the SCB, these changes in seawater chemistry can extend more than 100 km from the coast (∼30% of the Bight) (Kessouri, 2024). The region of maximum change occurs in the epipelagic zone, localized between 50 and 200 m water depth. When these declines are superimposed on areas already naturally low in O₂ and Ω A r , even small changes could be of biogeochemical and ecological significance (Levin, 2018; Roman et al., 2019). The question is whether these subsurface O₂ and acidification changes are occurring at ecologically relevant conditions, resulting in vertical compression of habitat. There are field-based demonstrated consequences of coastal acidification for shell-building zooplankton, in particular pteropods (Bednaršek et al., 2019; Feely et al., 2024). Similarly, O₂ depletion in the ocean can reduce metabolic performance of aerobic taxa (Fry, 1971; Pörtner and Knust, 2007; Seibel, 2011). Most of the literature on hypoxia focuses on acute lethal levels, but sublethal effects, even subtle ones that pose constraints on feeding times, can combine to limit growth or reproduction (Gunderson and Leal, 2016; Pinsky et al., 2020). Studies in other ecosystems have documented how short-term, low-O₂ events can give rise to immediate habitat compression of sensitive species, increasing susceptibility to overfishing [e.g. of brown shrimp and demersal fishes in the Gulf of Mexico (Craig, 2012) and of artisanal fisheries species in the Sea of Oman (Piontkovski and Al-Oufi, 2014)]. Even the behavior of smaller vertical-migrating taxa, like copepods, is shaped by seasonal O₂ and temperature (Pierson et al., 2017). On the longer-term, interactions between temperature and O₂ availability on aerobic metabolism have strong correspondence with faunal diversity, species distributions, predator-prey interactions, and changing biogeographic patterns (Deutsch et al., 2015, 2020; Seibel, 2016), and may even result in range shifts (Pinsky et al., 2020; Howard et al., 2020) and extinction (Penn and Deutsch, 2022).

In this study, we assess the degree to which modeled O₂ losses and acidification due to land-based nutrient inputs translates to changes in habitat capacity. To accomplish this, we rely on two metrics that define the habitat available for aerobic metabolism and for calcification. Aerobic habitat is determined using the Metabolic Index ( Φ ) , for which trait-based threshold varies across species. Aerobic habitat for northern anchovy, Engrualis mordax, is detailed but we also consider how aerobic habitat is modified for a range of species with metabolic traits that have differing oxygen and temperature sensitivities. Calcifying habitat is based on the saturation state of aragonite ( Ω A r ) , for which thresholds also vary among species. We evaluate calcifier habitat capacity as the thickness of the water column where Ω A r ≥ 1.4, but we also consider the sensitivity of our results to other values of Ω A r .

Our study objectives are threefold. First, we evaluate temporal and spatial patterns in aerobic and calcifier habitat capacity metrics in the SCB with output from a 20-year numerical ocean model hindcast. Second, we test how anthropogenic nutrient inputs from land-based sources alter the vertical thickness of the habitat capacity metrics. We rely on two model scenarios, the first includes natural oceanic cycles of nutrients, O₂ and carbon, to which rising global CO₂ emissions have been imposed (referred to hereafter as ‘CTRL’), and the second includes both natural oceanic cycles of nutrients and inputs from terrestrial sources, 98% of which are anthropogenic and 95% of which are point source in origin (referred to hereafter as ‘ANTH’) (Sutula et al., 2021). Third, we confirm the mechanisms by which anthropogenic nutrients contribute to the observed changes in vertical habitat capacity by analyzing changes in the biogeochemical rate processes that contribute to the O₂ and carbon cycles.

To assess whether modeled effects of land-based nutrient inputs on Bight-wide subsurface O₂ loss and acidification are biologically relevant, we employed two metrics for habitat capacity, which we adapted for this purpose. One incorporates temperature-dependent environmental O₂ as a predictor of habitat capacity for aerobic metabolism and the other incorporates carbonate chemistry as a predictor of habitat capacity for aragonite production by calcifiers. The premise of our approach is that these metrics provide information on the capacity of a specified location to provide habitat conditions that are sufficient for key processes for a species, or group of species, based on either empirical or mechanistic relationships of organismal performance with the environmental condition(s) of interest. The habitat capacity metrics are applied to model outputs from scenarios with and without land-based nutrient inputs in order to perform a difference assessment. Each metric is presented as the volume (or thickness of water-column) above a relevant threshold.

Since modeled effects of subsurface O₂ loss and acidification due to anthropogenic nutrient inputs are shown to be localized between 50 and 200 m (Kessouri, 2024), we focus our analysis on pelagic taxa. Because literature is limited on the interactive effects of O₂ and carbonate chemistry in these environments, we adapt two separate metrics to evaluate the changes in O₂ versus the changes in carbonate chemistry with an emphasis on aerobic taxa and calcifiers, accordingly.

For the effects of subsurface acidification on calcifier habitat capacity, we calculate the vertical thickness of optimal aragonite saturation state ( Ω A r ) conditions. A value of Ω A r of 1.4 is used to define the condition below which sublethal organismal responses have been documented to commonly occur (Bednaršek et al., 2019, 2021a, 2021b). One of the primary lines of evidence for this choice is derived from a synthesis of documented effects on pteropods (Bednaršek et al., 2019), in which Ω A r thresholds for a range of sublethal to lethal responses were identified and confidence in thresholds were judged with expert consensus. Pteropods are ubiquitous, holoplanktonic calcifiers that have a well-documented, specific sensitivity to ocean acidification. These calcifiers efficiently transfer energy from phytoplankton to higher trophic levels (Lalli and Gilmer, 1989; Hunt et al., 2008), and as such serve as an important prey group for ecologically and economically important fishes, bird, and whale diets (Armstrong et al., 2005; Aydin et al., 2005; Karpenko et al., 2007). Bednaršek et al. (2019) identified that Ω A r from 1.5 – 0.9 provides a risk range from mild dissolution to lethal impacts. Our selected value of Ω A r of 1.4 represents a value within recommended measurement precision for observational programs (± 0.2 (McLaughlin et al., 2015)) of thresholds where sublethal effects on calcification, growth, and severe dissolution are documented to occur, while a value of 1.0 roughly equates to lethal effects [0.9 to 0.95 (Bednaršek et al., 2019)]. In the epipelagic (0-200 m) in the Southern California Bight, conditions below saturation have not been common in the modern ocean (McClatchie et al., 2016), but are predicted to emerge as soon as the 2030’s and 2040’s (Hauri et al., 2013). However, undersaturation can occur over the benthos on the shelf (Kekuewa et al., 2022). Importantly, we perform an analysis of the sensitivity of our findings to the choice of Ω A r along this range of 1.0 to 1.4 and find the results to be largely in sensitive within this range (see Supplementary material for further details).

For the effects of subsurface O₂ loss on aerobic habitat capacity, we calculate the vertical thickness of the water column that has sufficient O₂ to provide ecological support for northern anchovy (Engraulis mordax). Northern anchovy is also holoplanktonic with greatest abundance observed in the upper 100 m (Mais, 1974). Defining sufficient O₂ for ecological support relies on the mechanistic framework of the Metabolic Index [ Φ (Deutsch et al., 2015, 2020)]. Φ is defined as the ratio of O₂ supply to resting demand. We calculate the habitat thickness for which Phi/PhiCRIT ≥ 1; this value demarcates environment's in which anchovy can sustain active energetic demands. In contrast, values below one can sustain resting but not active energetic demands, limiting population persistence. For northern anchovy, metabolic traits have been inferred from observational datasets associated with climatological O₂ and temperature conditions (Howard et al., 2020). Similarly, we evaluate a lethal threshold for northern anchovy, where Φ =   1 and O₂ supply is insufficient to meet O₂ demand. This latter analysis converts the reported active hypoxia threshold ( A o / Φ C R I T ) to the value at rest ( A o ) using Φ C R I T =   3.5 , the mean value across marine organisms (Deutsch et al., 2020); some species have Φ C R I T as low as 1.5, in which case the lethal thresholds could occur at proportionately higher values of O₂ and thus at shallower depths.

We utilize biogeochemical output from the Regional Ocean Modeling System, ROMS (Shchepetkin and McWilliams, 2005), coupled to the Biogeochemical Elemental Cycling model, BEC (Moore et al., 2004), which has been adapted for the CCS (Deutsch et al., 2021; Renault et al., 2021). BEC is a multi-element (C, N, P, O, Fe, and Si) and multi-plankton model that includes three explicit phytoplankton functional groups (picoplankton, silicifying diatoms, and N-fixing diazotrophs), one zooplankton group, and dissolved and sinking organic detritus. Remineralization of sinking organic material follows the multi-phase mineral ballast parameterization of Armstrong et al. (2001), and sedimentary processes have also been expanded. Particulate organic matter reaching the sediment is accumulated and remineralized with a time scale of 330 days, to provide a buffer between particle deposition and nutrient release. The ecosystem is linked to a carbon system module that tracks dissolved inorganic carbon and alkalinity, and an air–sea gas exchange module based on the formulation of Wanninkhof (1992).

The SCB model domain, which extends from Tijuana Mexico to Pismo Beach (U.S. Central California coast) and about 200-km offshore, is part of a nested configuration. Model nests scale from a 4-km horizontal resolution configuration spanning the entire CCS, to a 1-km resolution grid covering much of the California coast, to a 0.3-km grid in the SCB, where our investigations of local anthropogenic inputs were focused (Kessouri et al., 2021a, b). This grid, shown in Figure 1A , is composed of 1,400 × 600 grid points, with 60 σ-coordinate vertical levels using the stretching function described in Shchepetkin and McWilliams (Shchepetkin and McWilliams, 2005). The model is run with a time step of 30 s, and outputs are saved as 1-day averages. More information on the model setup and forcing is provided in other works (Deutsch et al., 2021; Renault et al., 2021; Kessouri et al., 2021a).

ROMS-BEC has been validated for atmospheric forcing, physics, and biogeochemistry including O₂, carbonate saturation state, primary productivity, and hydrographic parameters at a West Coast-wide scale (Deutsch et al., 2021) and, within the SCB, at scales at which anthropogenic nutrient inputs influence coastal eutrophication (Kessouri et al., 2021a).

We rely on two model scenarios. The first includes natural oceanic cycles of nutrients, carbon and oxygen cycles without terrestrial inputs (CTRL). The ocean carbon cycle includes the rising CO₂ signature from global observations. The second represents these same CTRL base conditions, to which inputs from terrestrial sources are added, 98% of which are anthropogenic and 95% of which are point source in origin (ANTH) (Sutula et al., 2021). Model simulations that include terrestrial inputs were forced with a monthly time series of spatially explicit inputs, including freshwater flow, nitrogen, phosphorus, silica, iron, and organic carbon representing natural and anthropogenic sources (Sutula et al., 2021). These data include ocean outfalls originating from Publicly Owned Treatment Works (POTW) and riverine discharges (1997–2017) and spatially explicit modeled estimates of atmospheric deposition. POTW effluent data were compiled from permit monitoring databases and communication with sanitary agencies. Monthly time series of surface water runoff from 75 rivers are derived from model simulations and monitoring data (Sutula et al., 2021). The model domain includes the U.S.-Mexico border and the simulations include not only U.S. land-based nutrient inputs but also Mexico transboundary flows from the Tijuana River watershed. The CTRL simulation covers the time periods of 02/1997 – 01/2001 and 08/2012 – 11/2017. The ANTH simulation covers the time period of 02/1997 – 11/2017. The overlapping time periods available for comparison between CTRL and ANTH represent a broad range of multi-year oceanographic conditions and climatic states.

Aragonite saturation state was computed with the CO2SYS algorithms (Lewis et al., 1998; Sharp et al., 2020) using daily averages of model output fields of dissolved inorganic carbon (DIC), total alkalinity (TA), temperature, salinity, and pressure. Calcifier habitat thickness was calculated as the thickness of the water column for each grid cell that was ≥ optimal Ω A r (1.4) for pteropods.

The ecological Metabolic Index ( Φ / Φ C R I T ) was computed from daily averages of model output fields of O₂ and temperature.

The metabolic traits of northern anchovy are esimated as A o / Φ C R I T =   5.4 atm^-1 (equivalent to an active hypoxia threshold of p O 2 =   0.185 atm at 15°C) and E o =   0.4 eV (the net temperature sensitivity of O₂ supply and demand) (Howard et al., 2020). p O 2 is the environmental partial pressure of O₂ and T is temperature (in K). k B is the Boltzmann constant and T r e f is the reference temperature (here, 288.15 K). Aerobic habitat thickness was calculated as the thickness of the water column where Φ / Φ C R I T ≥ 1 . An active hypoxia threshold of 0.185 atm O₂ at 15°C is equivalent to a saturation state of 88.5% and a concentration of 219 µmol kg^-1 (calculated at 15°C, S=35 [pss-78], and 0 dB gauge pressure). However, we note that organismal physiology is thought to be insensitive to the concentration of oxygen. Converting the above reference value to an in-situ concentration requires correction for local temperature and salinity, which vary over space and time. Thus, two different concentrations could refer to the same p O 2 , with no expected impact on physiological hypoxia, and vice versa.

Spatial and temporal patterns in calcifier and aerobic habitat thickness were evaluated with the ANTH simulation to identify the dominant spatial and temporal scales of variability in each. Total calcifier and aerobic habitat for the model domain was summed across all grid cells as the habitat thickness within a grid cell multiplied by the area of that grid cell.

To then evaluate how anthropogenic nutrient inputs alter calcifier and aerobic habitat thickness, we perform a difference assessment (ANTH-CTRL) where positive (negative) values represent an expansion (contraction) of habitat thickness attributable to anthropogenic nutrient inputs included in the ANTH scenario only. We focus further analyses in regions where differences in habitat thickness exceed ± 20%. We include a few sensitivity tests to some of the analytical choices. (1) To test the sensitivity of calcifier habitat capacity to the value of Ω A r , we evaluate habitat thickness for Ω A r ≥ 1.0   –   2.5 . (2) We also calculate the change in water-column thickness for a range of [O₂] from 60 – 200+ mmol m^-3. (3) Since northern anchovy, our species of focus for aerobic habitat, has its own temperature and oxygen sensitivity, we also explore how modeled aerobic habitat is changing for a range of temperature and oxygen sensitivities. We do this by estimating the aerobic habitat for both the ANTH and CTRL scenarios for the full combination of metabolic trait possibilities ( A o / Φ C R I T and E o ). The full range of traits and their probability distribution among species have been previously compiled (Penn and Deutsch, 2022). For each metabolic trait combination, available aerobic habitat is calculated as the volume of seawater where temperature and oxygen conditions translate to an ecological Metabolic Index between 1 and 1.6 (Penn and Deutsch, 2024). The change in aerobic habitat between the scenarios is derived as the difference between ANTH and CTRL. This change is then multiplied by the probability distribution of metabolic traits across species ( Supplementary Figure S1 ). The result is a distribution of aerobic habitat change across species. This approach assumes that the distribution of metabolic traits from the global-scale database corresponds with those in the California Current Ecosystem. We check this assumption by comparing the biogeographically inferred distribution of metabolic traits with that of respirometry-derived traits of those species found in the Southern California Bight and find that the two are indistinguishable ( Supplementary Figure S1 ).

For regions undergoing more than a ± 20% change in habitat thickness, we evaluate the difference in the biogeochemical rate processes from each scenario [detailed methods provided in Deutsch et al. (2021) and Kessouri (2024)]. Biogeochemical rate processes that influence the O₂ cycle include surface air-sea flux, photosynthesis, non-grazing mortality, grazing mortality, water-column remineralization, sediment-water flux, NH₄ oxidation, and nitrification [Equation A9 in Deutsch et al. (2021)]. Biogeochemical rate processes that influence dissolved inorganic carbon include air-sea flux, photosynthesis, CaCO₃ production, non-grazing mortality, grazing mortality, and water-column and sediment remineralization [Equation A11 in Deutsch et al. (2021)]. We perform a difference assessment from the monthly averages for the sum of the biogeochemical process terms. This analysis is focused between 70 and 140-m water depth to align with the depth range where habitat thickness is affected.

For both habitat metrics, there are consistent Bight-wide spatial and temporal patterns ( Figures 1A–E ), with calcifier habitat thickness greater than anchovy aerobic habitat thickness. Calcifier habitat thickness ranges from, on average, 80 to 130 m. Anchovy aerobic habitat thickness ranges, on average, from 50 to 100 m. Both habitat thickness metrics are most restricted within the Santa Barbara Channel and around San Nicolas Island. Increases generally occur along north-to-south and onshore-to-offshore gradients.

The dominant temporal scale of calcifier and anchovy aerobic habitat capacity are inter-annual and seasonal ( Figures 1C–E ). The temporal mean (across all years and seasons) of total calcifier and anchovy aerobic habitat volumes are 8.2 and 5.4 x 10³ km³, respectively, summed across the model domain. Among years, total volume for each metric can vary approximately 2-fold. In 2011, calcifier and anchovy aerobic habitat volumes were most restricted at 4.2 and 3.6 x 10³ km³; in 1998, habitat volumes were most expansive at 10.2 and 6.2 x 10³ km³, respectively. Seasonal variability also drives approximately 2-fold changes in habitat volumes. An evaluation of the annually detrended time series for both metrics shows that total habitat volumes are greatest during winter months and contract during summer, with the least amount of total habitat available during July and August ( Figures 1D, E ). Since anchovy aerobic habitat is constrained by both temperature and O₂, attribution analysis reveals that O₂ is the primary contributor to seasonal trends in total habitat volume, and that seasonal changes in aerobic habitat due to temperature counteract those changes due to that of O₂ ( Supplementary Figure S2 ).

There is not a consistent tendency for land-based nutrient inputs to result in consistent or persistent Bight-wide calcifier and aerobic habitat capacity gains or losses ( Supplementary Figure S3 ). The relative, Bight-wide differences in total habitat volume vary by ± 5% for most of the period simulated. However, at any grid cell location, the difference in calcifier habitat thickness between ANTH and CTRL can range from −34 to + 45%, and that for anchovy aerobic habitat ranges from −43 to + 80% (1^st and 99^th percentiles, respectively).

A spatial perspective of the change in calcifier and anchovy aerobic habitat thickness reveals a region of habitat loss that is expressed southeast of Santa Catalina Island ( Supplementary Figure S3 ). This region of habitat loss is a seasonal phenomenon. In seven of the nine years simulated, there is a compression event lasting approximately 2.5 months that occurs in the late summer to early fall ( Figures 2C, D ), a time period when calcifier and anchovy aerobic habitat thickness are already seasonally compressed ( Figures 1D, E ). Averaging across these seven compression events, the total spatial area experiencing recurrent calcifier habitat change is 2,364 km² (assessed as the total spatial area where the compression in calcifier habitat thickness > 20%; Figure 2A ) The equivalent total spatial area experiencing recurrent anchovy aerobic habitat change is 1,909 km² ( Figure 2B ). Despite the temporal overlap in natural seasonal and eutrophication-driven compression, there is spatial mismatch between the regions undergoing maximum habitat compression due to eutrophication with those that are naturally most restricted from broad-scale oceanographic patterns (e.g., Santa Barbara Channel; Figures 1A, B ).

Anchovy represent one ecophysiotype in a range of metabolic trait possibilities. The region of > 20% habitat compression for anchovy does not translate equally across all species ( Figure 3A ). In this region, ecophysiotypes with higher hypoxia tolerance ( A o / Φ C R I T ) gain aerobic habitat volume; while those with lower hypoxia tolerance, including anchovy, lose aerobic habitat volume. Taking a species-weighted distribution of metabolic trait possibilities, approximately two-thirds of species are losing aerobic habitat volume, and the modeled losses can be up to three times greater than the gains ( Figure 3B ). Of those subjected to aerobic habitat loss, 70% of species are subject to more habitat loss than northern anchovy. Species that are gaining aerobic habitat volume have higher hypoxia tolerance at any given temperature sensitivity ( E o ; Figure 3A ).

The spatial and vertical distribution of greatest O₂ change is occurring between approximately 40 and 100 km from the coast and between 50 and 150 meters or more below the surface ( Figures 4A, C ). The same pattern is observed for acidification, assessed as the difference in Ω A r between the two scenarios ( Figure 4B ). The lower limit of anchovy aerobic habitat coincides in space with the region of maximum O₂ loss. In the CTRL scenario, the anchovy aerobic habitat limit is at 80 m in the offshore region of maximum O₂ loss (60 km from the coast), and shoals to 60 m in the ANTH scenario. In contrast, the lower limit for calcifier habitat is deeper than the vertical region undergoing maximum acidification. Still, calcifier habitat shoals from approximately 115 m in the CTRL scenario to 87 m in the ANTH scenario. Because both subsurface acidification and O₂ loss occur across a broad depth range, our results are largely insensitive to which value of Ω A r is used to define optimal calcifier habitat ( Supplementary Figure S4 ). Similarly, loss of O₂ is occurring across a broad depth range ( Supplementary Figure S5 ). Notably, these changes to habitat capacity in the epipelagic are largely limited to sublethal effects as conditions that trigger acute lethal effects occur deeper in the water column. Conditions where Ω A r are less than 1 occur, on average, deeper than 200 m water depth ( Supplementary Figure S4 ). Similarly, acute lethal O₂ conditions ( Φ =   1 ) for northern anchovy occurs, on average, much deeper than 300 m water depth, well below the typically observed vertical distributions of anchovy and other pelagic fishes.

Of the biogeochemical rate processes contributing to habitat change within the region of 20% habitat compression, remineralization rates are exhibiting the greatest absolute change due to land-based nutrient inputs ( Figure 5 ). For the carbon cycle, remineralization rates increase from 4.26 ± 0.12 to 4.50 ± 0.12 mmol DIC m^-3 d^-1 from the CTRL to the ANTH scenario (mean ± 1 SE; N = 104 months). The increase in DIC from land-based nutrient inputs at the core of habitat compression drives the modeled decrease in Ω A r . Small changes in alkalinity (< 5 mmol m^-3) counteract the effects due to DIC ( Supplementary Figure S5 ). For the O₂ cycle, remineralization rates are also exhibiting the greatest change due to land-based nutrient inputs, changing from −5.44 ± 0.15 to −5.75 ± 0.15 mmol O₂ m^-3 d^-1 from CTRL to ANTH (mean ± 1 SE; N = 104 months).

Here, we demonstrate that eutrophication effects of land-based nutrient export to the Southern California Bight are not restricted to chemical changes in seawater acidification and O₂ loss (Kessouri, 2024), but also extend to the potential for widespread (biotic and ecological) effects on calcifier and aerobic habitat capacity. The seawater chemistry changes that occur in the epipelagic are not at levels that elicit acute, lethal effects. Rather, the sublethal habitat capacity metrics used here exhibit annually recurring compression despite large, natural seasonal and interannual cycles ( Figure 1 ). During the late summer, subsurface acidification and O₂ loss (from land-based nutrient export) routinely compress aerobic and calcifier habitat capacity ( Figure 2 ), at a time period when habitat capacity is already seasonally compressed ( Figure 1 ). Modeled habitat compression is most pronounced where excess nutrients and organic matter, which originate at the coast, are received and entrained within offshore eddies (Kessouri, 2024). These locations occur in both state and federal waters and overlap with many marine protected areas, including those around Santa Barbara Island, Santa Catalina Island, and San Clemente Island. Since seawater chemistry changes due to enhanced remineralization are occurring across a large depth range, patterns in habitat compression are largely insensitive to the value of Ω A r used to define the calcifier habitat capacity metric ( Supplementary Figure S4 ). Similarly, while we evaluate aerobic habitat capacity for northern anchovy, we confirm that this pattern is consistent among two thirds of ecophysiotypes (although at differing magnitudes of habitat loss), and that those species gaining aerobic habitat volume have higher tolerances to low O₂ ( Figure 3 ).

There is field-based evidence that the vertical structure of both Ω A r and O₂ have implications for a variety of pelagic taxa. For example, across frontal gradients in the California Current System where Ω A r between 1.0 and 1.4 can shoal by 100+ m on the scale of tens of kms, there are concurrent reductions in pteropod abundance accompanied by elevated shell dissolution (Bednaršek and Ohman, 2015). From onshore-to-offshore gradients, more severe pteropod shell dissolution and thinner shells occur close to the coast, particularly where upwelling is more intense (Feely et al., 2016; Mekkes et al., 2021). The predicted habitat compression described here coincides with the natural seasonal cycle of limited habitat availability. Multiple species of pteropods, including Limacina helicina, are present year-round in the SCB (K. McLaughlin, pers. comm.). While limited baseline information on pteropod life history characteristics exist for the Southern California Bight (Manno et al., 2017), studies suggested that spring (April-May) and fall (September-October) are periods when early life stage cohorts are most vulnerable to changing ocean conditions (Gannefors et al., 2005; Wang, 2014; Wang et al., 2017).

Ocean O₂ depletion adversely impacts marine species, assemblages, and even fisheries (Laffoley and Baxter, 2019). Long-term deoxygenation trends play a role in, for example, declines in abundance of mesopelagic fishes (Koslow et al., 2011) and shifts in zooplankton and small nekton diel migration depth (Bianchi et al., 2013). Interactive effects of sub-optimal O₂ and temperature are becoming increasingly considered (Roman et al., 2019). Sub-optimal O₂ stress depends on the O₂ supply relative to metabolic demand, and water temperature controls both chemical (O₂ solubility, diffusivity) and physiological processes (metabolic demand, ventilations rates) affecting this balance for marine ectotherms. We use the mechanistic framework of the Metabolic Index (Deutsch et al., 2015) to incorporate these dependencies into the index of aerobic habitat capacity. Specific to our focal taxa, northern anchovy have seasonal to interdecadal redistributions that correlate with aerobic habitat capacity. For example, anchovy migrate offshore during peak upwelling seasons (Mais, 1974; Laroche and Richardson, 1980), when nearshore aerobic habitat availability is lowest, even though their food supply is generally higher closer to the coast (Howard et al., 2020). Further, the southern biogeographic limit of this species is coincident with the aerobic habitat capacity threshold implied by their oscillations in time within the SCB, and vice versa (Howard et al., 2020). While reductions in this index are associated with species-specific consequences of deoxygenation at the regional scale, it remains unclear how the spatial and seasonal extent of the O₂ reduction identified here might translate to disruptions across species of varying phenologies, mobility, and ecological niches. The same can be said for the population-level consequences of subsurface acidification.

We evaluate O₂ loss and acidification effects on habitat capacity separately, as the combined effects of these stressors on biological responses are insufficiently understood (Bednaršek et al., 2020). However, studies suggest that exposure to suboptimal ranges of acidification, O₂, and temperature can make marine organisms more sensitive to O₂ loss (Breitburg et al., 2019) or less resilient to acidification (Stevens and Gobler, 2018). In this model domain, there is strong covariance between calcifier and aerobic habitat thickness (Pearson correlation coefficient = 0.79), linked to eutrophication effects on water-column remineralization, such that aerobic and calcifier habitat compress at the same time. Thus, predicted effects on habitat capacity for marine species may be underestimated (Bednaršek et al., 2020). The Metabolic Index framework does incorporate the combined effects of temperature and O₂. Biological sensitivities to all three variables – temperature, O₂, and carbonate system state – could potentially be merged through a fundamental physiological trait such as aerobic scope [i.e., a proxy for the surplus energy available for growth, reproduction, predator avoidance, etc.; en sensu Sokolova (2021)].

While we emphasize subsurface losses in pelagic habitat capacity, changes to subsurface biogeochemistry will also intersect with the complex seafloor topography of the Southern California Bight. We have yet to evaluate how anthropogenic nutrient loads are altering habitat availability for demersal organisms, particularly those that have limited mobility. Another consideration to be further explored is the enhanced productivity and food supply from anthropogenic nutrient loads. Increased food supply could potentially reduce the negative consequences of suboptimal O₂ and acidification (Breitburg et al., 2009; de Mutsert et al., 2016). As an example, total fisheries landings can remain high even if demersal species in O 2 − depleted areas decline, because nutrients can stimulate prey production in other well-mixed parts of a system (Breitburg et al., 2009; Nixon and Buckley, 2002). Abundant prey can improve stress tolerance of organisms (Marigomez et al., 2017). However, system-wide compensation through enhanced productivity will have limits as the volume of O₂- depleted waters expand (Breitburg, 2002). The catch per unit effort for selected demersal fish species along the U.S. West Coast is positively related to near-bottom O₂ concentrations, with the catch per unit effort decreasing more significantly as O₂ concentrations decrease (Keller et al., 2015). In the Humboldt Current, it is suggested that low O₂ near the coast results in a highly efficient trophic transfer and a dense anchovy population (Bertrand et al., 2011), which is beneficial for fishing activities. During these conditions, species that are less tolerant to low O₂, like sardine and jack mackerel, are restricted to offshore, well-oxygenated waters (Alegre et al., 2015; Bertrand et al., 2016).

Global climate change will further exacerbate habitat loss resulting from land-based nutrient inputs (Kessouri et al., 2021b). Strengthened stratification, from increased surface water temperatures as the global climate warms, is sufficient to worsen subsurface O₂ and acidification where it currently exists and may instigate habitat loss elsewhere (Long et al., 2019). Warming and O₂ loss by 2100 are projected to result in complete loss of aerobic habitat for northern anchovy – and thus likely extirpation – from the southern California Current System (CCS) (Howard et al., 2020). Further, the interplay of anthropogenic nutrient export and stratification where they occur could accelerate the timeline of habitat compression and potential extirpation. In this study, O₂ loss in the core of habitat compression exceeds 30 µmol kg^-1 ( Figure 4C ). This loss due to eutrophication is 1.5 times the scale of decadal O₂ loss in the southern CCS which may be occurring at around 20 µmol kg^-1 decade^-1 (Bograd et al., 2008). Acidification in the same eutrophication-induced core exceeds −0.3 units for Ω A r , and this decrease due to eutrophication is three times the decadal trend of approximately −0.1 decade^-1 (Leinweber and Gruber, 2013; Turi et al., 2016).

To conclude, we assess change in habitat capacity for pelagic calcifying and aerobic taxa due to eutrophication effects on subsurface acidification and O₂ loss from land-based nutrients. Our findings suggest that effects of land-based nutrients are not restricted to chemistry. Changes to habitat capacity defined by sublethal, ecologically relevant thresholds were pervasive during late summer, when habitat capacity is at its seasonal minimum. Despite the theoretical, experimental, and field evidence that identify the importance of the vertical structure of both carbonate chemistry and O₂ for marine pelagic communities, whether the modeled habitat compression shown here translate to population-level effects is uncertain. Species abundance and distributional data from existing monitoring programs, like the California Cooperative Oceanic Fisheries Investigations and the Southern California Bight Regional Monitoring Program, could be assessed through the lens of eutrophication-driven habitat compression predicted here. These same monitoring programs are developing and implementing biological metrics for acidification. Coupling region-wide modeling with experimental and field programs will be a robust foundation for study of highly urbanized coastlines.