Anoxia begets anoxia: A positive feedback to the deoxygenation of temperate lakes

Declining oxygen concentrations in the deep waters of lakes worldwide pose a pressing environmental and societal challenge. Existing theory suggests that low deep-water dissolved oxygen (DO) concentrations could trigger a positive feedback through which anoxia (i.e., very low DO) during a given summer begets increasingly severe occurrences of anoxia in following summers. Specifically, anoxic conditions can promote nutrient release from sediments, thereby stimulating phytoplankton growth, and subsequent phytoplankton decomposition can fuel heterotrophic respiration, resulting in increased spatial extent and duration of anoxia. However, while the individual relationships in this feedback are well established, to our knowledge, there has not been a systematic analysis within or across lakes that simultaneously demonstrates all of the mechanisms necessary to produce a positive feedback that reinforces anoxia. Here, we compiled data from 656 widespread temperate lakes and reservoirs to analyze the proposed anoxia begets anoxia feedback. Lakes in the dataset span a broad range of surface area (1–126,909 ha), maximum depth (6–370 m), and morphometry, with a median time-series duration of 30 years at each lake. Using linear mixed models, we found support for each of the positive feedback relationships between anoxia, phosphorus concentrations, chlorophyll a concentrations, and oxygen demand across the 656-lake dataset. Likewise, we found further support for these relationships by analyzing time-series data from individual lakes. Our results indicate that the strength of these feedback relationships may vary with lake-specific characteristics: For example, we found that surface phosphorus concentrations were more positively associated with chlorophyll a in high-phosphorus lakes, and oxygen

Declines in bottom water DO concentrations are often attributed to climate change and/or increased nutrient inputs (Bartosiewicz et al., 2019;Jane et al., 2023;Jenny, Francus, et al., 2016).Increased air temperatures have been shown to drive increased duration of thermal stratification (Foley et al., 2012;North et al., 2013;Oleksy & Richardson, 2021;Woolway et al., 2021), which reduces or inhibits mixing of oxygen to the bottom waters (hypolimnion).Consequently, increases in stratification duration may provide more time for hypolimnetic DO depletion to occur, resulting in lower late-summer DO concentrations and increased duration of anoxia.Changes in stratification duration appear to be a particularly important driver of DO declines in recent decades (ca. 1950-2020;Jane et al., 2023).
However, historical nutrient inputs have likely also played a role in deoxygenation by increasing phytoplankton biomass and consequently oxygen demand (Jenny, Francus, et al., 2016;Jenny, Normandeau, et al., 2016).The relative importance of these two pathways to deoxygenation (i.e., greater stratification duration due to climate change and greater oxygen demand due to eutrophication) likely varies both among lakes and within lakes over time.Consequently, understanding interannual DO dynamics across many lakes may be critical to disentangling the independent effects of stratification duration and eutrophication amidst ongoing changes in global climate and land use (e.g., Moss, 2011;Parmesan et al., 2022).
Here, we analyze a positive feedback, derived from decades of aquatic research, by which anoxia (i.e., DO at or near 0 mg/L) during a given year begets increasingly frequent and severe occurrences of anoxia in subsequent years.In this "anoxia begets anoxia" (ABA) feedback, anoxic conditions promote internal phosphorus release, thereby stimulating phytoplankton growth and subsequent decomposition, which in turn fuels increased heterotrophic respiration and further accelerates hypolimnetic DO declines over time (Figure 1).
As long-term limnological data have become increasingly accessible (e.g., Jane et al., 2021;Pilla et al., 2020), we now have the opportunity to test the strength and ubiquity of this feedback on a multi-continental scale.
While the individual relationships in the ABA feedback cycle (Figure 1) are well established, these relationships occur over multiple timescales and amidst numerous other interacting factors (e.g., climate variation) that could prevent the detection of the overall feedback.Hypolimnetic anoxia has been shown to enhance internal loading of phosphorus from sediments (e.g., Mortimer, 1941;Nürnberg, 1984;Orihel et al., 2017;Figure 1a).However, while redox-controlled phosphorus release fluxes have received significant attention, sediment characteristics, microbial processing, and demand had a stronger influence on the extent of anoxia in deep lakes.Taken together, these results support the existence of a positive feedback that could magnify the effects of climate change and other anthropogenic pressures driving the development of anoxia in lakes around the world.

K E Y W O R D S
air temperature, anoxia, chlorophyll a, dissolved oxygen, feedback, hypolimnion, lake, oxygen demand, phosphorus, residence time F I G U R E 1 The proposed positive feedback through which "anoxia begets anoxia" (ABA).Hypolimnetic anoxia results in internal hypolimnetic total phosphorus (TP) loading (a), which in turn increases epilimnetic TP (b) and stimulates phytoplankton growth, resulting in increased chlorophyll a (chl a; c).Phytoplankton decomposition fuels increased oxygen demand rates (d), which further drive hypolimnetic oxygen declines (e).This feedback can be externally influenced by increased air temperatures (gray dashed lines), among other factors.
Increases in hypolimnetic total phosphorus (TP) are expected to increase surface water (epilimnetic) TP concentrations within a summer stratified period through both biological and physical processes (e.g., organism-mediated transport, diffusion, and internal seiche dynamics; Carpenter et al., 1992;Cottingham et al., 2015;Haupt et al., 2010;Kamarainen et al., 2009) or during autumn mixis when epilimnetic and hypolimnetic waters homogenize (e.g., Nürnberg & Peters, 1984;Wetzel, 2001;Figure 1b).Higher epilimnetic TP concentrations in turn can stimulate phytoplankton growth in many lakes, thereby increasing chlorophyll a (chl a, Figure 1c; Schindler, 1974), though many other important factors, including nitrogen concentrations, climate, and light availability, also contribute to phytoplankton growth (e.g., Paerl & Huisman, 2008;Reinl et al., 2023).Increased phytoplankton biomass and subsequent decomposition may fuel increased biological oxygen demand (Figure 1d; Ladwig et al., 2021;Müller et al., 2019;Pace & Prairie, 2005) and result in earlier onset of anoxia (Figure 1e), although climate can also play an important role in driving DO dynamics in many lakes, as discussed above.Given the substantial complexity to each of these relationships, all operating on different timescales, it remains unclear the extent to which the full positive feedback plays a role in controlling DO dynamics within lakes around the world.Lake characteristics including size and residence time could potentially mediate the strength of the ABA feedback across lakes, though these relationships remain largely untested because they can only be characterized with long-term monitoring data across many diverse lakes.Lakes with longer residence time or larger sediment area:volume ratios may have greater sediment-water interactions, increasing the influence of oxygen demand on hypolimnetic DO, as well as the influence of hypolimnetic DO on hypolimnetic TP (e.g., Jagtman et al., 1992).Likewise, lake size may control the importance of mixing dynamics between the epilimnion and hypolimnion, and residence time may affect the extent to which chl a and hypolimnetic TP influence biogeochemical dynamics the following year (Wetzel, 2001).While many of these expected relationships have not been assessed across lakes, an empirical analysis of data from 2849 lakes suggests that the impact of phosphorus concentrations on chl a may be stronger in shallow lakes relative to deep lakes, potentially due to differences in light availability and macrophyte cover (Zhao et al., 2023).Characterizing the effect of lake characteristics on the ABA feedback relationships is needed to identify which lakes are most susceptible to the feedback, enabling managers to prioritize conservation efforts across lakes.
In this study, we analyzed data from 656 widespread temperate lakes to study the drivers and consequences of interannual changes in hypolimnetic DO.Our research had three primary goals: First, we assessed the extent of support for each of the hypothesized relationships between anoxia, hypolimnetic phosphorus concentrations, epilimnetic phosphorus concentrations, epilimnetic chl a, and oxygen demand across and within lakes (Figure 1).Second, we analyzed records of air temperature at each lake to assess how the ABA feedback may interact with changes in climate (Figure 1).We focused on climate as an external driver of the ABA feedback in lieu of accessible nutrient loading records for the study lakes.Third, we analyzed whether the strength of ABA relationships may vary with lake characteristics including lake depth and residence time.While our multi-lake approach precluded detailed consideration of external nutrient inputs and use of causal inference methods within a lake over time, analyzing data from many lakes was essential to testing the proposed relationships in this study and disentangling lake-specific effects amidst substantial heterogeneity.

| Overview of data compilation and analysis
Analyzing the ABA feedback required time-series data for hypolimnetic DO, hypolimnetic TP, epilimnetic TP, epilimnetic chl a, hypolimnetic oxygen demand, and climate records across numerous lakes (Figure 1).We compiled in-lake data from 656 geographically widespread stratified lakes to enable these analyses (Section 2.2).
We used linear mixed models, including relevant lags and climatic data when appropriate (Section 2.3.2) to assess support for the ABA feedback relationships across all lakes.We then ran the same linear models within individual lakes when sufficient data were available to assess whether the strength of ABA relationships may vary with lake characteristics (Section 2.3.3).All data compilation and analyses are described in detail below.

| In-lake data
We synthesized data from a total of 656 temperate, seasonally stratified lakes (Figure 2; Appendix S1: Text S1.1).Data were collated from Jane et al. (2021;n = 316 unique lakes not also available in the other datasets described here), the U.S. Wisconsin Department of Natural Resources (DNR; n = 163), the U.S. New Hampshire Volunteer Lake Assessment Program (VLAP; n = 93), the U.S. Lake Stewards of Maine (LSM) Volunteer Lake Monitoring Program (n = 48), the U.S. Adirondack lakes database (Leach et al., 2018;Winslow et al., 2018;n = 17), and members of the Global Lake Ecological Observatory Network (GLEON; n = 29).Chl a data from Filazzola et al. (2020) were added for n = 15 lakes for which we did not have any other chl a data.
Data availability and collection methods differed substantially among sites (documented in Lewis et al., 2023).For each site, we collated available data for DO, water temperature, TP, and chl a, as well as lake metadata including geographic coordinates, depth (mean and maximum), surface area, and elevation (Lewis et al., 2023).Total nitrogen (TN) and dissolved organic carbon (DOC) were also compiled, but were more limited in availability (n = 111 lakes for DOC and n = 119 lakes for TN), motivating us to primarily focus on TP in our analyses below.To harmonize multiple datasets, quality control was performed on all data, as described in the data publication (Lewis et al., 2023).
In sum, the complete dataset consisted of 108,736 distinct water temperature and DO profiles across 656 lakes during 1938-2022 (Appendix S2: Figure S2.1).The median data duration was 30 years at each lake (range: 3-81 years).Lakes in the dataset had a median depth of 14 m (Z max ; range: 6-370 m), median surface area of 100 ha (range: 1-126,909 ha), and median elevation of 264 m (range: −215-2804 m).The lakes were located in 18 countries across five continents, with latitudes ranging from −42.6 to 68.3 (Lewis et al., 2023).

HydroLAKES
We collated additional metadata for each lake using HydroLAKES, a global database of 1.4 million lakes (with surface area ≥10 ha; Messager et al., 2016).For lakes with missing mean or maximum depth (i.e., the depths were not reported with the data; n = 43), we used HydroLAKES data to fill in these values (Lewis et al., 2023).We also compiled residence time estimates from HydroLAKES to assess whether the strength of ABA feedback relationships may vary with differences in residence time across lakes.

Profile interpolation
We interpolated all temperature and DO profiles to a 1-m resolution following Jane et al. (2021).Briefly, we selected all profiles with at least three depths, then used the pchip() function of the pracma R package (Borchers, 2022) to interpolate measurements from the surface to the deepest sampled depth.
To account for variation and error in sampling procedures, we implemented a standardized screening protocol to remove temperature and DO profiles that were substantially shallower or deeper than the reported maximum depth of the lake (Appendix S3).

Mean concentrations
We averaged data for all focal variables to an annual timestep using data from the entire stratified period and, separately, the late-summer period at each lake (Appendix S1: Text S1.2).The late summer (i.e., mid-July through August in the northern hemisphere, following Jane et al., 2021) is when DO concentrations are likely to approach their lowest value (Wetzel, 2001), and may consequently be a critical time period for some processes in the ABA feedback.Conversely, other processes occurring throughout the entire summer stratified period (e.g., oxygen demand, hypolimnetic temperature) can also be critical to the ABA feedback, motivating the study of both periods within a year.
For each profile during either the entire summer stratified period or the late-summer period, we calculated the depths of the top and bottom of the metalimnion (the middle thermal layer of the lake) using the rLakeAnalyzer R package (Winslow et al., 2019).We used mean metalimnion depths to estimate the bottom of the epilimnion and top of the hypolimnion for each lake year.We then averaged all hypolimnetic and epilimnetic water quality measurements throughout the time period of analysis, using interpolated profiles for temperature and DO and all measurements for TP, chl a, TN, and DOC.To estimate the strength of stratification at the thermocline, F I G U R E 2 Data were compiled from a total of 656 widespread temperate lakes, with data availability differing across sites.(a) Map of all sites included in this dataset.Note that due to overlapping data points, many sites are not visible.More detailed maps of the United States and Europe are provided in Appendix S2: Figures S2  and S3.Map lines delineate study areas and do not necessarily depict accepted national boundaries.(b): Summary of data availability for water temperature, dissolved oxygen (DO), total phosphorus (TP), and chlorophyll a (chl a) in the epilimnion and hypolimnion of lakes in this study.
we calculated maximum buoyancy frequency using rLakeAnalyzer (Read et al., 2011;Winslow et al., 2019) for each temperature profile.
Maximum buoyancy frequency was averaged throughout the stratified period for each lake year (Table 1).

| Volume-weighted hypolimnetic oxygen demand
We calculated volume-weighted hypolimnetic oxygen demand (VHOD; hereafter oxygen demand) within each lake year, following Wetzel and Likens (2000).Briefly, we used measured or modeled bathymetric contours and interpolated DO profiles to calculate the volume-weighted hypolimnetic DO concentration for each sampling date, then used linear regression models to calculate the rate of decline in volume-weighted hypolimnetic DO concentrations within the summer stratified period (Burns, 1995;Håkanson, 2005;Quinlan et al., 2005;Wetzel & Likens, 2000;Appendix, S4).We calculated an oxygen demand rate based on the raw data, as well as a temperature-corrected oxygen demand rate following Pace and Prairie (2005).Detailed methods for both calculations are provided in Appendix S4.

TA B L E 1
Explanatory variables used for mixed model regression.We tested several possible explanatory variables for each response variable using a mixed model approach.The time period over which mean values were calculated for each lake year is provided for all water column variables.For information on lags used, see Appendix S7: Figures S1-S5.Epilimnion and hypolimnion are abbreviated as epi.and hypo.throughout.

Anoxic factor
Oxygen

| Anoxic factor
Anoxic factor (AF) describes the spatial and temporal extent of anoxia within a lake and is therefore a useful metric of deoxygenation in lakes that experience hypolimnetic anoxia (Nürnberg, 1995(Nürnberg, , 2019)).AF is expected to increase with increased oxygen demand, and can predict internal TP loading in lakes that experience hypolimnetic anoxia (Nürnberg, 1995(Nürnberg, , 2019; Figure 1).Here, we calculated AF within each lake year following Nürnberg (1988) and Nürnberg et al. (2019), modified to address limited data availability across and within lakes (Appendix S5).Briefly, we estimated the duration of anoxia using oxygen profiles, oxygen demand, and modeled turnover dates, and we used modeled or measured bathymetry to quantify the spatial extent of anoxia within each lake year.The DO threshold for anoxia was defined operationally, as described below (Section 2.3.3), with detailed methods provided in Appendix S5.

| Climate data
To disentangle the roles of changing climate and in-lake processes on DO dynamics in stratified lakes, we collated monthly air temperature and precipitation data for every lake in our dataset from the ERA5 climate reanalysis.ERA5 is a fifth-generation product from the European Centre for Medium-Range Weather Forecasts (ECMWF) and provides meteorological data from 1959 to 2022 on a 0.25-degree global grid (Hersbach et al., 2019).For our analysis, we used the monthly 2-m air temperature and total precipitation ERA5 data products, and found the closest gridded values for every lake in our dataset.We summarized "seasonal" air portion of the ice-covered period in lakes that experience seasonal ice cover (Magnuson et al., 2000).

| Data analysis
To analyze the proposed ABA relationships, we used lag analysis (Section 2.3.1),mixed effects modeling (Section 2.3.2), and within-lake regressions (Section 2.3.3).All data analyses were performed in R, version 4.2.1 (R Core Team, 2022).Analysis code is archived as a Zenodo repository for reproducibility (Lewis & Lau, 2023).

| Lag analysis
Several of the relationships in the proposed ABA feedback may operate across years, rather than within 1 year.To assess the appropriate lag for each step, we calculated the Spearman correlation between each variable of analysis and the preceding variable in the feedback cycle (e.g., between oxygen demand and chl a; Figure 1) with 0-, 1-, and 2-year lags.These correlations were calculated separately for each lake with at least 10 years of paired data for the target parameters.Across all lakes, we calculated whether the mean of the resulting distribution of correlations was significantly different than zero using Wilcox tests with α = .05.

| Mixed effects modeling
To assess the proposed mechanisms by which anoxia could create a positive feedback that promotes subsequent anoxia (Figure 1), we used linear mixed models to estimate the magnitude and direction of effect for drivers of AF, epilimnetic and hypolimnetic TP, epilimnetic chl a, and oxygen demand among lake years.To assess the relationship between oxygen demand and hypolimnetic DO concentrations in lakes that did not experience anoxia (i.e., AF = 0 days throughout the entire time series), we conducted an additional regression analysis for oxygen demand and late-summer hypolimnetic DO concentrations, rather than AF (Appendix S6).Lake ID was included as a random effect on the intercept in all models.
Mixed effect models were run using the package lme4 in R (Bates et al., 2023).
For each response variable, we filtered all data to only include lake years with complete data for all proposed explanatory variables (Table 1).We log-transformed chl a and TP concentrations due to the substantial positive skew of these data, and we Zstandardized all explanatory variables.We fit linear mixed models for all possible combinations of explanatory variables and identified the best model using corrected Akaike information criterion (AICc).We report all selected models within two AICc units of the best model (Burnham & Anderson 2002).We assessed the multicollinearity of all models using the variance inflation factor, which we calculated using the vif() function from the package car in R (Fox et al., 2022).
We plotted the coefficient estimate for all fixed effects in the selected models to visually compare the magnitude of effect for each explanatory variable.For these visualizations, we calculated 95% confidence intervals of the fixed effects using the confint.merMod()function from lme4 (Bates et al., 2023).

| Operational definition of anoxia
We used an operational DO threshold to define hypolimnetic anoxia, following other studies on anoxia in lakes (e.g., Elshout et al., 2013;LaBrie et al., 2023;Nürnberg et al., 2019).To identify this threshold, we performed a breakpoint analysis and piecewise regression for hypolimnetic DO and TP using the package segmented in R (Muggeo, 2023; Appendix S6: Text S6.1).
We then added slope difference (U) and change point (G 0 ) parameters for the breakpoint relationship, and used the resulting breakpoint as a threshold value for our calculation of AF (Appendix S5).

| Within-lake regressions
To assess whether the across-lake trends identified using mixed models were observable within individual lakes, we performed linear regressions separately at each lake.For each of our focal response variables (see Table 1), we used the same model formulations from the across-lakes analysis (i.e., the explanatory variables from Table 1 that were selected via AICc) to perform regressions within a lake.We saved the resulting coefficient estimates for each explanatory variable used to predict this focal response.We then plotted the distribution of coefficient estimates for all explanatory variables across all lakes, and we compared the median of these distributions to the mixed effect model coefficient estimates.For each response variable, we only included lakes that had at least 10 years of paired data for the response variable and all selected explanatory variables.We removed n = 81 lakes that never experienced anoxia (i.e., AF = 0 throughout the timeseries) from the within-lake analysis of the drivers of AF.

Driver analysis
The coefficient estimates for explanatory variables included in the ABA feedback (e.g., the coefficient of epilimnetic TP for predicting epilimnetic chl a) indicates the magnitude of the response, while accounting for other drivers (Table 1).As an exploratory analysis to assess which lakes are most susceptible to the ABA feedback, we analyzed whether there were significant differences in these coefficients based on differences in lake characteristics.For this analysis, we developed linear models predicting the coefficient estimate for each focal variable in the ABA feedback (Table 1) based on (individually) maximum depth, surface area, mean depth, residence time, dynamic ratio (square root of lake area divided by mean depth; Håkanson, 1982), and mean concentrations of focal (ABA) variables (i.e., hypolimnetic DO, epilimnetic and hypolimnetic TP, epilimnetic chl a, and oxygen demand).We then used AICc to select the model(s) with the greatest explanatory power.We did not assess more complicated model structures (e.g., multiple drivers and interaction effects) due to the relatively small sample size for some of these analyses (e.g., n = 35 lakes for oxygen demand).

| Climate effects
To summarize the effects of climatic variation on oxygen dynamics, we analyzed monthly and annual air temperature data.First, we calculated correlations between monthly air temperatures and, separately, hypolimnetic temperature, oxygen demand, AF, and late-summer DO concentrations (Appendix S8).Then, we summarized the effects of high and low annual air temperature anomalies on AF and late-summer oxygen concentrations (Appendix S8).

| Operational definition of hypolimnetic anoxia
We identified a breakpoint relationship whereby hypolimnetic TP increased substantially after DO decreased below a threshold of 1.8 mg/L (56 μmol/L), averaged throughout the hypolimnion (Figure 3).Subsequently, we used 1.8 mg/L as our DO threshold for anoxia in all analyses.Of the 356 lakes with at least 10 years of hypolimnetic DO data, 146 lakes (34%) crossed the threshold of 1.8 mg/L during their time series (i.e., had at least 1 year with hypolimnion-averaged DO <1.8 mg/L and at least 1 year with DO

F I G U R E 3
Piecewise mixed model regression identified a breakpoint in the relationship between hypolimnetic dissolved oxygen (DO) and total phosphorus (TP) at 1.8 mg/L DO.Here, points represent individual lake-years.

| Regression analyses support expected relationships within and across lakes
Our analyses across 656 lakes provided support for the ABA feedback.Of the explanatory variables used in our model selection process (Table 1), all variables that were predicted to promote the ABA feedback were found to be statistically significant drivers of their predicted responses (Figure 4), with expected temporal lags as applicable (0-1 years; Appendix S7).High AF was associated with high hypolimnetic TP (Figure 4a), and high hypolimnetic TP was associated with high epilimnetic TP, both within and between years (i.e., both Hypo TP and Hypo TP t−1 had positive coefficients; Figure 4b).
High epilimnetic TP was in turn associated with high chl a within a year (Figure 4c), and high chl a was associated with high oxygen demand (both VHOD and VHOD std 10°C ) the following year (Figure 4d; Appendix S10).Lastly, high oxygen demand was associated with greater AF in the lakes that experienced hypolimnetic anoxia (Figure 4e).For the lakes that did not exhibit anoxia during their time series, high oxygen demand was associated with low late-summer DO concentrations (Appendix S6).
All of the ABA relationships observed to be significant across hundreds of lakes (n = 111-386; Figure 4) were also supported by regression analyses conducted within individual lake time series (with n = 35-157 lakes for each analysis; Figure 5).The direction of each of the ABA relationships was identical within and across lakes (Figure 5).The magnitude of the median coefficient estimates for The proposed "anoxia begets anoxia" (ABA) feedback (bottom right) was supported by linear mixed model results across all variables (panels a-e; see Table 1).Here, panel titles indicate the response variable for each panel and y-axis labels indicate explanatory variables.X-axes indicate the magnitude and 95% confidence interval of the parameter estimate for each explanatory variable presented on the y-axis.The black vertical lines in panels (a-e) denote a parameter estimate of zero.Blue rectangles highlight drivers in the hypothesized ABA feedback (bottom right).Explanatory variables are ordered by the magnitude of the parameter estimate within each panel.chl a, chlorophyll a; epi., epilimnion; hypo., hypolimnion; TP, total phosphorus.
ABA explanatory variables within lakes (e.g., the coefficient for chl a in the multiple linear regression with oxygen demand as a response variable) tended to be slightly smaller than the mixed model coefficient estimate (Figure 5) for each relationship, except for oxygen demand as a predictor of AF (Figure 5e).
While the hypothesized ABA feedback was supported by regression analyses, variability in the focal response variables (i.e., AF, TP, chl a, and oxygen demand; Table 1) was also modulated by additional driving factors, as expected (Figure 1; Appendix S8).Specifically, climatic variables were selected as part of the optimal model for nearly all focal variables: spring air temperatures were important drivers of AF and chl a, spring and summer precipitation were significant drivers of epilimnetic TP, and winter precipitation was a significant driver of hypolimnetic TP (Figures 4 and 5).Water temperature also played a role in explaining variation in several focal responses: Hypolimnetic temperatures were a significant predictor of both AF and oxygen demand (Figures 4 and 5).For all responses, we found substantial variability in the random intercept of the mixed-model regressions among lakes (Table 2) and variability in within-lake regression coefficients (Figure 5), indicating external lake-specific factors that influence the state of each response variable at a given lake.Random effects were largest for AF, and residual standard deviation from mixed-model analyses was highest for oxygen demand and epilimnetic chl a (Table 2).Across lakes, our analyses indicate that the relative strength of ABA relationships varied with lake characteristics.Specifically, the coefficient for the effect of epilimnetic TP on chl a was larger for lakes with high mean epilimnetic TP values; the coefficient for the effect of oxygen demand on AF was larger for lakes with deep mean depth; and the coefficient for the effect of chl a on oxygen demand was larger for lakes with long residence time (Figure 6).
The other ABA feedback relationships were not significantly mediated by any one of our candidate predictors (see "Driver analysis" section).

| DISCUSS ION
In analyzing ABA relationships both across and within 656 lakes, we found support for all linkages in the hypothesized ABA feedback (Figures 4 and 5).These results provide empirical support for the existence of a positive feedback mechanism that could intensify the development of anoxia in lakes around the world.
Furthermore, our results indicate that the strength of these relationships likely varies with lake characteristics, including mean depth, TP concentrations, and residence time.To our knowledge, our work is the first to quantitatively document all of the relationships that enable anoxia to beget increasingly frequent or more intense anoxia in future years across a large, multi-continental dataset of lakes.

| Decades of research facilitate identification of ABA feedback
Individual relationships in the ABA feedback have been the subject of substantial research inquiry over the past century or longer (e.g., Sachs, 1874;Schindler, 1974;Thienemann, 1928).While these previous studies primarily focused on examining biogeochemical dynamics within one lake, they provided support for the individual relationships in the ABA feedback (Figure 1).Modeling studies provided a means of simultaneously considering all ABA relationships, and have shown mechanistic support for the existence of an ABA feedback in seasonally stratified lakes (Carpenter, 2003;Carpenter & Lathrop, 2008).However, model simulations have indicated that the susceptibility of individual lakes to a trophic regime shift, as a result of the ABA feedback, depends on multiple lake-specific parameters (i.e., macrophyte presence, temperature, mean depth; Genkai -Kato & Carpenter, 2005), highlighting the need for a multilake empirical approach.
By synthesizing data across many lakes, our mixed model approach allowed us to identify biogeochemical dynamics that likely would have been difficult to detect in individual lakes.The strength of this approach is reflected in the fact that coefficient estimates from our mixed model regressions, which integrate data from many lakes, were typically slightly larger in magnitude than the median coefficient estimates of regressions run within individual lakes (Figure 5), although both approaches showed support for the The strength of "anoxia begets anoxia" feedback relationships may be modulated by lake characteristics.(a) The coefficient for the effect of epilimnetic total phosphorus (epi.TP) on chlorophyll a (chl a) was most positive in lakes with high mean epi.TP.(b) The coefficient for the effect of the previous year's chl a on volume-weighted hypolimnetic oxygen demand (VHOD) was most positive in lakes with long residence times.(c) The coefficient for the effect of VHOD on anoxic factor (AF) was most positive in lakes with deep mean depths.This relationship was robust to including all data (solid regression line) and excluding disproportionately influential points (i.e., Cook's distance greater than 3× the mean, n = 12 lakes; shown as a dashed line).Linear regressions are presented as solid lines.
existence of the ABA feedback.Across-lake regressions included a larger range of variation for predictor variables than is typically observed within individual lakes, which likely facilitated the detection of more substantial predictor-response effects.Through the study of the hypothesized ABA feedback, we found support for several individual limnological relationships, some of which had not been previously analyzed on a widespread scale.Below we discuss our findings for each ABA relationship and their implications in the context of previous work (Sections 4.1.1-4.1.5).
4.1.1| Effect of anoxia on hypolimnetic TP (Figure 1a) In this study, we observed a strong positive relationship between hypolimnetic anoxia and TP concentrations both within and across lakes.
Across lakes, our breakpoint analysis detected a threshold relationship whereby hypolimnetic DO had a stronger effect on TP when DO concentrations decreased to levels approaching anoxia (<1.8 mg/L; Figure 3).Our results reinforce previous research affirming that AF (the duration and spatial extent of anoxia) may be strongly positively correlated with hypolimnetic TP concentrations (Figures 4 and 5; e.g., North et al., 2014;Nürnberg et al., 2019).A threshold relationship between DO and TP is well supported by previous research across sediment core incubations, in situ sediment chamber measurements, and mass-balance whole ecosystem analyses (e.g., Anderson et al., 2021;Einsele, 1936;Mortimer, 1942;Orihel et al., 2017).Here, our threshold value of 1.8 mg/L DO, averaged throughout the entire hypolimnion, likely reflects DO conditions of ~0 mg/L near the sediment-water interface (which inherently is challenging to quantify empirically), resulting in enhanced TP loading (Nürnberg, 2019).We note that our identified breakpoint of 1.8 mg/L is also remarkably similar to those identified in previous sediment incubation work (Doig et al., 2017;Matisoff et al., 2016;Orihel et al., 2017).Overall, this analysis indicates that the ABA mechanism may require hypolimnetic DO concentrations to decrease to low levels (i.e., <1.8 mg/L) before a feedback effect will occur.
In our dataset, it was common for lakes to cross the threshold of 1.8 mg/L (34% of n = 356 lakes).Lakes where oxygen concentrations declined below 1.8 mg/L had lower DO concentrations in the year following the onset of anoxia than in the year prior to the onset of anoxia (Appendix S9: Figure S9.1).While our dataset was not a random or fully representative sample of global lakes, the large number of lakes which crossed the 1.8 mg/L threshold in this study suggests that the ABA feedback may be prevalent.
4.1.2| Effect of hypolimnetic TP on epilimnetic TP (Figure 1b) We found moderately strong support for an effect of hypolimnetic TP on epilimnetic TP both within 1 year and between years (i.e., hypolimnetic TP influences epilimnetic TP the following year).
While the directionality of this relationship can be difficult to identify in the absence of detailed nutrient input data (i.e., epilimnetic TP can affect hypolimnetic TP, vice versa, or a third driver may simultaneously influence both), existing research provides strong support for this effect.Elevated hypolimnetic TP concentrations can increase epilimnetic TP concentrations within a summer stratified period through organism-mediated transport, diffusion, and internal seiche dynamics (e.g., Carpenter et al., 1992;Cottingham et al., 2015;Haupt et al., 2010;Kamarainen et al., 2009;Nürnberg, 2009;Soranno et al., 1997).
At the onset of autumn mixing, the concentration of TP in the hypolimnion fundamentally determines the amount of potential TP input to the epilimnion, which can have legacy effects throughout the subsequent autumn, winter, and spring (e.g., Nürnberg & Peters, 1984;Wang et al., 2019).
4.1.3| Effect of epilimnetic nutrients on epilimnetic chl a (Figure 1c) We found a strong positive association between surface water TP concentrations and surface water chl a, both within and across lakes, likely reflecting the fact that interannual variability in phosphorus concentrations can play an important role in regulating phytoplankton growth in lakes (Figures 4 and 5).Our study follows many decades of data that illustrate the positive effect of TP on phytoplankton biomass (MacKeigan et al., 2023;Schindler, 1974;Smith, 1982).In this study, we were unable to identify an effect of epilimnetic TN concentrations on chl a, suggesting that in these lakes, TP may play a more important role in regulating phytoplankton growth.However, we note that data availability was substantially greater for TP (n = 387 lakes) than for TN (n = 86 lakes), and complexities of nitrogen forms (not considered here) may hinder the detection of a nitrogen effect.Previous research has documented the importance of nitrogen for limiting or co-limiting phytoplankton growth in some lakes, over multiple timescales (Elser et al., 2007;Lewis et al., 2020;Lewis & Wurtsbaugh, 2008;Paerl et al., 2016;Scott et al., 2019).Consequently, our study highlights the need for long-term, speciated nitrogen data to disentangle the role of nitrogen in the ABA feedback.
4.1.4| Effect of epilimnetic chl a on oxygen demand (Figure 1d) Support for the relationship between epilimnetic chl a and oxygen demand was relatively weaker than for the other ABA relationships, although still consistent within and across lakes.We expected that this relationship would be more challenging to detect than the other ABA relationships due to high levels of spatio-temporal heterogeneity in chl a and uncertainty associated with oxygen demand calculations (e.g., modeled bathymetry and the assumption of a closed system).Interestingly, the effect of chl a appeared to occur at least as strongly between years as within a year.Legacy effects of chl a on oxygen demand are intuitive and expected, as decomposition of sediment organic matter (including settled phytoplankton biomass) may constitute the majority of the total hypolimnetic oxygen demand in many lakes (Steinsberger et al., 2020).Likewise, limited sampling of early-season bloom events could have partially obscured the role of within-year chl a on oxygen demand.Regardless, our analyses provide support for both within-year and betweenyear effects of phytoplankton blooms in perpetuating anoxia.
4.1.5| Effect of oxygen demand on hypolimnetic anoxia (Figure 1e) The positive relationship between oxygen demand and AF is well supported by this study, and is also intuitive: as biological and chemical demand for oxygen increases, the onset of anoxia is likely to occur earlier in the stratified period, increasing the total duration of anoxia (Figures 4 and 5).Furthermore, in lakes that did anoxia (e.g., Jane et al., 2023;Nürnberg, 1995).Previous work has identified that the duration of summer stratification is increasing across many lakes (Woolway et al., 2021), driving decreased latesummer oxygen concentrations (Jane et al., 2023).However, the factors that control oxygen demand are changing less consistently: Temporal trends in hypolimnetic temperature are highly variable across lakes (Pilla et al., 2020;Richardson et al., 2017), as are trends in chl a from 1980 to present (Kraemer et al., 2022).
Consequently, it is not surprising that trends in oxygen demand appear to be inconsistent across lakes (Jane et al., 2023).In this study, our focus on annual and subannual timescales allowed us to more precisely investigate the mechanisms at play within and across 386 lakes (Figure 4e), identifying that variability in oxygen demand has the potential to drive a feedback effect in some lakes that experience hypolimnetic anoxia.

| Lake characteristics can increase susceptibility to the ABA feedback
Through our cross-lake analyses, we identified that the ABA feedback may be stronger in some lakes than others.In particular, mean epilimnetic TP concentrations, mean depth, and residence time each modulated ABA feedback relationships (Figure 6).
First, the effect of TP on chl a was strongest in lakes with high mean epilimnetic TP concentrations, especially for lakes with TP concentrations greater than ~10 μg/L (Figure 6a).These mesotrophic to eutrophic/hypertrophic lakes also tended to experience substantial variability in epilimnetic TP concentrations, which likely made the effect of changing TP concentrations more detectable in our standardized linear regression analyses (Appendix S11: Figure S11.1).
Ultimately, our finding that TP and chl a are more closely correlated at high TP concentrations may provide some resistance to the initiation of the ABA feedback in oligotrophic lakes, while further accelerating the ABA feedback as eutrophication proceeds due to external or internal nutrient loading.
Second, the effect of the previous year's chl a on oxygen demand was strongest in lakes with long residence times (Figure 6b).In these lakes, decomposing chl a and autochthonous organic carbon may have more time to settle and accumulate on the hypolimnetic sediments, fueling oxygen demand the following year.Conversely, the effect of the previous year's chl a on oxygen demand was negligible in lakes with residence time less than ~100 days (Figure 6b), as chl a may be quickly flushed and exported downstream from these lakes.Consequently, lakes with longer residence time may be more susceptible to the ABA feedback.
Third, the magnitude of the effect of oxygen demand on AF generally increased with increasing mean depth of the lake (Figure 6c).
Mechanistically, deeper lakes often have relatively lower oxygen demand due to low sediment area to hypolimnetic volume ratios (Livingstone & Imboden, 1996;Müller et al., 2012;Steinsberger et al., 2020).Consequently, variation in oxygen demand can substantially affect the amount of time it takes to reach anoxia in these deep lakes.Conversely, in shallow lakes, hypolimnetic DO concentrations may be more strongly impacted by factors other than oxygen demand, including hypolimnetic primary production, stratification phenology, and mixing events (Wetzel, 2001).Ultimately, deep lakes (i.e., mean depth >5 m; Figure 6) appear to have a particularly strong coupling between oxygen demand and AF, strengthening the ABA feedback in these lakes.
Combined, these results suggest that deep mesotrophic or eutrophic lakes with long residence times are particularly likely to be susceptible to the ABA feedback, though more data are needed to test these hypotheses.Importantly, our identification of factors that may affect the strength of the ABA feedback across lakes would not have been possible without the use of a multi-lake dataset like the one analyzed in this study.et al., 2021).Specifically, the time period during which temperatures fall in the historical range of spring temperatures is shortening across Northern Hemisphere mid-latitudes, which are representative of most of the lakes in this study (Wang et al., 2021).Conversely, the time period during which temperatures fall in the historical range of summer temperatures is lengthening (Wang et al., 2021;Woolway, 2023).Our work highlights the importance of accounting for these differential changes in seasonal air temperatures, not just annual means, when anticipating how changes in climate may affect hypolimnetic DO dynamics.Furthermore, as spring air temperatures continue to increase across many lakes, our work suggests that these climatic changes may play a role in causing hypolimnetic oxygen concentrations to decline, potentially initiating the ABA feedback.

| Strengths and limitations of this analysis
Using regression models within and across lakes, we were able to simultaneously analyze the extent of support for each of the relationships in the hypothesized ABA feedback.Lakes analyzed in this study span five orders of magnitude in surface area and two orders of magnitude in maximum depth (Z max ; Lewis et al., 2023).
Amidst these substantial differences, we found consistent support for the ABA feedback relationships within and across lakes.
While the dataset analyzed here is larger than those used in previous studies, data limitations continued to constrain our analysis.

| CON CLUS IONS AND G LOBAL CHANG E IMPLIC ATIONS
We found widespread empirical support for the ABA feedback in temperature and precipitation values by averaging across multiple months for each lake year, with southern hemisphere data offset by 6 months.Spring values were calculated as the average of March and April air temperature or precipitation (following, e.g.,Williamson et al., 2015).While stratification onset varies across latitudes and lakes, these spring months are the most likely to correspond to ice melt and spring mixing across the temperate lakes in this study(Woolway et al., 2021; Appendix S1: Figure S1.2).Summer values were calculated as the average of July and August air temperature or precipitation, as these summer months most closely correspond with our late-summer in-lake data and were the warmest 2 months on average across the dataset (Appendix S2: Figure S2.4).Winter temperature and precipitation were calculated as the average of January and February air temperature and precipitation.These winter months were, on average, the coldest months in our dataset (Appendix S2: Figure S2.4), and likely constituted a significant

≥1. 8
mg/L).Lakes that crossed this threshold (n = 146) were more common than lakes that had consistently anoxic (n = 120) or consistently oxic (n = 90) hypolimnia.Furthermore, lakes that crossed the threshold of 1.8 mg/L had lower DO concentrations in the year following the first year of anoxia than in the year prior to the first year of anoxia (Appendix S9: Figure S9.1).

F
Linear regressions analyzing time-series data within individual lakes provide further support for the proposed "anoxia begets anoxia" (ABA) feedback.Here, panel titles indicate the response variable for each panel (a-e) and y-axis labels indicate explanatory variables.Individual points represent regression coefficients from within one lake.Density distributions describe the distribution of parameter values across lakes, with colors delineating the quartiles of the distribution (purple: 0%-25%, blue: 25%-50%, green: 50%-75%, and yellow: 75%-100%).Black and white circles at the bottom of each distribution mark the parameter estimate from the mixed model analysis (Figure4).Gray vertical lines denote a parameter estimate of zero.Blue rectangles highlight drivers in the hypothesized ABA feedback.Explanatory variables are ordered by the magnitude of the mixed-model parameter estimate for consistency with Figure4.All x-axes range from −1 to 1 to enable comparison among panels.chl a, chlorophyll a; epi., epilimnion; hypo., hypolimnion; TP, total phosphorus.TA B L E 2 Random and residual variation from linear mixed models.Model structure and fixed effects are summarized in Figure4.
not experience anoxia throughout the time series of data used in this study, oxygen demand was negatively associated with latesummer DO concentrations (Appendix S6: Text S6.2), supporting that oxygen demand and DO concentrations are closely coupled in both oxic and anoxic lakes.Across the dataset, the effect of oxygen demand on hypolimnetic oxygen conditions occurred simultaneously with an additional positive effect of spring air temperatures (Figures 4 and 5; Appendix S6: Text S6.2), and in anoxic lakes AF was further regulated by autumn air temperatures (Figures 4 and 5).Positive associations between anoxia and spring and autumn air temperatures may highlight the important role that stratification duration (i.e., both onset in spring and end in autumn) can play in driving the spatial and temporal extent of Specifically, we were unable to analyze the effects of external nutrient loads, or DOC concentrations on the ABA feedback due to lack of data, and we were unable to use causal inference methods to study ABA dynamics within individual lakes over time.Moreover, the majority (82%) of lakes analyzed here are temperate lakes located in the United States; consequently, results may not be fully generalizable to global lakes, and more research is needed to characterize DO dynamics in a broader, representative range of ecosystems, especially in tropical and southern hemisphere lakes.Our calculated AF values have substantial uncertainty, particularly with respect to stratification end dates, though we have done our best to minimize these uncertainties through detailed methodological testing (Appendix S5).To standardize across a wide range of lakes and sampling regimes, our analysis considered the entire hypolimnion as one homogenized layer, averaging over potentially meaningful variation in DO dynamics across a depth gradient in the hypolimnion (e.g.,LaBrie et al., 2023).Given the promising results we observed here, further exploration of depth-resolved DO declines across lakes likely has substantial potential to further our understanding of biogeochemical processing in lakes.
analyzing time-series data across 656 diverse lakes.Relationships were particularly strong between oxygen demand and AF; AF and hypolimnetic TP; and epilimnetic TP and chl a. Conversely, the effect of epilimnetic chl a on oxygen demand was comparatively less strong, though still detectable both within and across lakes.As oxygen concentrations are decreasing in many lakes around the world, accounting for the ABA feedback may help effectively prioritize restoration and conservation efforts.Notably, our work suggests that catchment-scale nutrient management may be particularly critical for preventing deterioration of water quality in lakes with late-summer hypolimnetic DO concentrations just above 1.8 mg/L that have not yet crossed this threshold.These lakes are less likely to currently experience feedback effects of anoxia, but may cross this threshold in the future, thereby initiating an ABA feedback that, once triggered, will make water quality management more challenging.As climate and land use continue to change on a global scale, understanding and accounting for the ABA feedback may enable more effective conservation of culturally, economically, and ecologically important lake ecosystems.Conceptualization; data curation; formal analysis; funding acquisition; investigation; methodology; project administration; software; visualization; writing -original draft.Maximilian P. Lau: Conceptualization; data curation; investigation; methodology; project administration; software; writing -review and editing.Stephen F. Jane: Conceptualization; data curation; writing -review and editing.Kevin C. Rose: Conceptualization; data curation; writing -review and editing.Yaron Be'eri-Shlevin: Data curation; writing -review and editing.Sarah H. Burnet: Data curation; validation; writing -review and editing.François Clayer:

Climate change has the potential to trigger the ABA feedback
(Jane et al., 2023;Shatwell et al., 2019;Woolway et al., 2021)nship between climate variation and deoxygenation.Importantly, this climate variability may have the potential to push hypolimnetic of the hypolimnion(Jane et al., 2023;Shatwell et al., 2019;Woolway et al., 2021).While mean air temperatures are increasing around the world as a result of anthropogenic climate change, these impacts are not consistent across seasons or locations (Masson-Delmotte