DOC, grazers, and resilience of phytoplankton to enrichment

Phytoplankton blooms often follow nutrient enrichment. Differences among lakes in light‐absorbing dissolved organic carbon (DOC) may shift bloom thresholds to higher nutrient loads and thereby increase resilience of lakes to enrichment. To explore this idea, we measured resilience to experimental enrichment of two lakes with contrasting DOC concentrations. We compared bloom thresholds in both lakes using a model of phytoplankton response to DOC and nutrients, a dynamic time series indicator of resilience, and two empirical measures of stochastic resilience, mean exit time and median survival time. For the dynamic indicator and ecosystem model the lake with higher DOC was more resilient to enrichment. However, the distributions overlapped for stochastic indicators of resilience of the two lakes. These analyses show that DOC interacts with mixing depth and zooplankton biomass to affect resilience. Strong contrasts in DOC and many observations are needed to discern effects of DOC on resilience to enrichment.

states are shaped by the interaction of nutrient enrichment with grazers and physical properties of the water column (Soranno 1997;Kosten et al. 2012;Isles et al. 2015).
The contrasting states are separated by an unstable threshold (Scheffer 2009). Each state is locally stable because equilibrium is restored after small disturbances. Large disturbances, however, may cause the ecosystem to cross the threshold to the alternate state (Carpenter et al. 2022).
Many lakes are simultaneously experiencing increases in water color and nutrient loading (Leech et al. 2018). Water color from chromophoric dissolved organic carbon (DOC) is mainly terrestrial in origin (Wilkinson et al. 2013) and concentrations change in response to hydrogeology, vegetation, and precipitation in the watershed (Gergel et al. 1999;Zwart et al. 2016;Carpenter and Pace 2018). Inputs of DOC affect stratification, phytoplankton growth, and response to nutrient enrichment by altering light and nutrient availability (Jones et al. 2005;Rinke et al. 2010;Solomon et al. 2015). Higher loads of DOC decreased stability of mesocosm food webs and increased their sensitivity to nutrient pulses (Jones and Lennon 2015). Terrestrial DOC is accompanied by nitrogen and phosphorus some of which supports phytoplankton growth (Qualls et al. 1991;Kissman et al. 2017;Corman et al. 2018). Primary production has a hump-shaped response to increasing DOC concentrations due to tradeoffs between nutrient enhancement and shading effects (Kelly et al. 2018;Olson et al. 2020). This response seems dependent on carbonnutrient stoichiometry , and correlations of DOC, color, and limiting nutrients are variable among watersheds and lakes (Lapierre et al. 2021;Stetler et al. 2021).
Resilience measures the tendency of a state, such as the lowphytoplankton state of a lake, to persist despite changes in its environment, such as nutrient supply (Holling 1973). While stability measures the rate of recovery from local perturbations (Ives and Carpenter 2007), resilience considers effects of repeated and large disturbances, environmental trends, and the possibility that the ecosystem could cross thresholds to alternate states (Holling 1996;Scheffer et al. 2015). If the low-phytoplankton state is resilient to a specified nutrient input then it does not cross the threshold to the high-phytoplankton state. Thus resilience is the response of a specified state of the ecosystem to a specified environmental perturbation . Resilience of phytoplankton to nutrient inputs depends on nutrient stoichiometry, herbivory, and food-web structure (Elser et al. 1998;Carpenter et al. 2001).
Since DOC affects phytoplankton stability and response to nutrient enrichment, we hypothesized that the resilience of phytoplankton to nutrient enrichment may depend on DOC as well as grazing. This study compares resilience to nutrient enrichment of phytoplankton biomass in two experimental lakes with contrasting DOC for which we have detailed data on phytoplankton response to enrichment (Wilkinson et al. 2018). We apply three approaches to assess resilience: (1) a model of chlorophyll responses to nutrient enrichment, DOC, and grazers; (2) comparative resilience of phytoplankton to experimental enrichment of two contrasting lakes while monitoring dynamic indicators of resilience (Ives and Dakos 2012;Scheffer et al. 2015); and (3) comparisons of stochastic resilience of the low-phytoplankton state to nutrient enrichment Carpenter et al. 2022). We find that higher DOC may increase resilience of chlorophyll to enrichment but this effect is complicated by thermocline depth, grazer biomass, and stochasticity.

Lake descriptions
Whole-lake enrichments were conducted on Peter and Tuesday Lakes located in Gogebic County, Michigan, USA (46 250 N, 89 500 W). Before enrichment Peter and Tuesday lakes had relatively low primary production (Carpenter and Kitchell 1993;Carpenter et al. 2005). Both lakes are surrounded by fringing bogs and hardwood-conifer forests. A third lake, Paul Lake, served as undisturbed reference ecosystem.
Although the lakes are < 1 km apart they differ in some key limnological characteristics. Prior to this experiment we measured DOC, total phosphorus (TP), total nitrogen (TN), and their ratios during years with no enrichment, 2003-2012. Tuesday Lake had higher concentrations of DOC, TP, and TN during these unenriched years (Supporting Information Fig. S1).
During this study, the mean thermocline during summer stratification was deeper in Peter Lake (1.9 m) than in Tuesday Lake (1.2 m). DOC and chlorophyll mean summer concentrations were higher in Tuesday Lake (9.7 mg L À1 and 17.0 μg L À1 , respectively) than in Peter Lake (5.4 mg L À1 and 7.9 μg L À1 ). Zooplankton biomass was relatively low in Tuesday Lake (1.8 μg C L À1 ) vs. Peter Lake (12.2 μg L À1 ) consistent with food web contrasts established during earlier ecosystem experiments Pace et al. 2013). Overall, nutrient concentrations in this study and prior manipulations  indicate phytoplankton limitation by nutrients, particularly phosphorus.

Limnological methods
During 2013-2015 we made weekly measurements of DOC (Carpenter et al. 2017a) and nutrients (Carpenter et al. 2017b). A portion of filtered water sample was preserved with 200 μL of 1 mol L À1 H 2 SO 4 per 20 mL of sample and analyzed for DOC with a Shimadzu 5050 TOC analyzer (Shimadzu). Zooplankton biomass was estimated daily in 2015 by duplicate daytime verticals tows of an 153 μm mesh conical net, towed through the upper two-thirds of the water column over the deepest point of the lake Zooplankton samples were concentrated onto preweighted filters, dried, and reweighted to determine biomass .
Phycocyanin, dissolved oxygen, and pH were measured at 5-min intervals using Hydrolab.
DS5X sensors (OTT Hydromet) sensors suspended at a depth of 0.75 m from a centrally located buoy in Paul, Peter, and Tuesday lakes (Pace et al. 2020b). Daily phycocyanin was the mean sensor value between 22:00 h to 04:00 h to avoid quenching (Rousso et al. 2021). Daily ranges (maximum À minimum for each day) were calculated for dissolved oxygen saturation and pH. A manual sample for chlorophyll analysis was taken near the buoy each day between 1000 and noon and stored in a dark cooler. The sample was filtered (pore size = 0.7 μm) and extracts from the filters were analyzed fluorometrically for chlorophyll a (Chl a) concentration (Holm-Hanson 1978;Pace et al. 2020a). We used thermistor chains with temperature sensors every half meter. The 5-min temperatures were averaged to daily values to determine mixed layer depth. Nutrient additions were made by dissolving ammonium nitrate (NH 4 NO 3 ) and phosphoric acid (H 3 PO 4 ) iñ 20 L of lake water and then distributing the mixture into the epilimnion of each lake in the propeller wash of an underway electric motor. Nutrient additions occurred between 11:00 h and 13:00 h daily. Descriptions of enrichments for 2013-2015 are presented by Wilkinson et al. (2018).

Overview of stability and resilience analyses
We report three analyses of stability and resilience for two enriched lakes with contrasting DOC: (1) analysis of a model of phytoplankton response to DOC, grazers and enrichment fit to daily data from 2015, (2) multivariate time series analysis of daily data during 2015 to determine the time of transition from a low-phytoplankton state to a high-phytoplankton state, and (3) calculation of two indices of stochastic resilience, mean exit time and median survival time, using high-frequency data from 2013 to 2015. Below we summarize each resilience method and then the data used for the analyses.

Model of phytoplankton response to DOC and enrichment
We developed a simple ecosystem model to assess DOC and enrichment effects on phytoplankton (see section "Phytoplankton response to DOC, enrichment, and grazing" in Supporting Information S1). Daily observed time series were phytoplankton concentration A (μg C L À1 ), zooplankton biomass H (μg C L À1 ), DOC concentration C (mg L À1 ), and thermocline depth Z T (m) (Supporting Information Fig. S3). C, A and Z T are necessary to estimate m, the mean irradiance in the mixed layer using extinction coefficients measured for these lakes. A is measured as chlorophyll concentration (μg L À1 ) and converted to carbon units (μg L À1 ) using the measured C : Chl mass ratio of 60 (Carpenter et al. 2016). For model projections the time step Δt was 0.1 d. We fit the model to data and calculated 1-d projections using the difference equation: DOC affects growth through mean irradiance in the mixed layer m (Supporting Information Eqs. S1-S4) and daily P load (P) determines phytoplankton carrying capacity through parameter p. The model rationale, parameter definitions and estimates, and methods of stability analysis are presented in section "Phytoplankton response to DOC, enrichment, and grazing" in Supporting Information S1.

Effect of DOC on thermocline depth
The model of phytoplankton response to DOC, nutrients, and grazing requires an equation to predict thermocline depth from DOC and P load rate. Thermocline depth affects the phytoplankton model through the mean irradiance of the mixed layer (Supporting Information Eqs. S4-S6) and dilution of areal P loads to calculate P concentration. We estimated the effect of DOC, P load rate, and lake fetch on depth of the thermocline by regression analysis of eight nearby lakes for years between 1999 and 2016 when all variates were measured between 15 July and 15 August (see section "Estimating thermocline depth from lake area, DOC and P load" in Supporting Information S1).

Time series indicators of resilience
We assessed stability at each daily time step of 2015 using multivariate time series models (Ives and Dakos 2012) estimated by Bayesian updating (Pole et al. 1994). Multivariate time series models (Supporting Information Eqs. S13-S15) were fit to daily series of four observed variates (Supporting Information Fig. S4): log 10 phycocyanin, log 10 Chl a, delta dissolved oxygen saturation (daily maximum À daily minimum), and delta pH (daily maximum À daily minimum). Models are described in section "Multivariate time series analysis to assess stability" in Supporting Information S1 and worked examples with data and R scripts are presented by Carpenter et al. (2022).

Stochastic resilience
We computed mean exit time and median survival time ) to measure stochastic resilience. Mean exit time from a state is the mean time until the ecosystem crosses the threshold to an alternate state, and median survival time is the median time that the ecosystem occupies a specified state. Stochastic resilience accounts for the random fluctuations of the data, whereas resilience indicators from the multivariate time series analysis and fits of the ecosystem model provide point estimates of the threshold. For details of the calculations see "Stochastic measures of resilience" in Supporting Information S1.

Data for analyses
The phytoplankton model and multivariate time series analysis used daily observations from 2015 (Pace et al. 2020a) including chlorophyll concentration, zooplankton biomass, mixed layer depth, and DOC. 2015 was the only enriched year with daily zooplankton biomass data needed for the phytoplankton model. During 2015, daily enrichments to Peter and Tuesday lakes were identical (3 mg P m À2 d À1 with N: P ratio 15 : 1) starting on day of year 152. In 2015, Peter bloomed earlier and enrichment was ended on day 180 (Pace et al. 2017). Tuesday bloomed later in the summer and enrichment ended on day 240 (Wilkinson et al. 2018).
High-frequency measurements of phycocyanin from all three enriched years (Pace et al. 2020b) were used to compute stochastic resilience (exit time) because of the large number of observations required by the method (Carpenter et al. 2022).

Phytoplankton model
The phytoplankton model uses chlorophyll to indicate biomass because we have extensive data for the light extinction coefficient of chlorophyll (Carpenter et al. 1998). Nonetheless chlorophyll and phycocyanin are highly correlated at the daily scale of the dynamic model (Supporting Information Fig. S8; for Peter Lake r = 0.944, Tuesday Lake r = 0.891) indicating the importance of Cyanobacteria in the observed blooms. Covariates (DOC, thermocline depth, and zooplankton biomass) appear in (Supporting Information Fig. S2).
One-day predictions of chlorophyll concentration track the observations (Fig. 1A,B) and time series of predicted and observed 1-d change in chlorophyll are similar (Fig. 1C,D). The largest differences in predicted and observed 1-d change, and the greatest variability, occur near the unstable thresholds (approximately day 175 in Peter Lake and day 219 in Tuesday Lake). Analyses of residuals are presented in section "Phytoplankton response to DOC, enrichment, and grazing" in Supporting Information S1.
Model predictions are consistent with alternate states for conditions that we have observed in each lake (Fig. 1E,F). The net growth line of each panel crosses zero at three points, stable equilibria at low and high phytoplankton biomass and an intermediate unstable threshold. We used the model to compare effects of DOC and zooplankton biomass on the critical threshold for daily phosphorus load. The critical P load to a given lake is the lowest P load where a high-pigment equilibrium appears. At the critical P load a single low equilibrium for chlorophyll concentration transitions to three equilibria. The calculation of the critical threshold for P load required a regression model to predict thermocline depth from DOC and P load (see section "Estimating thermocline depth from lake area, DOC, and P load" in Supporting Information S1). For specified values of DOC, zooplankton biomass, and thermocline depth the critical threshold for P load is calculated from the phytoplankton model (see section "Equilibria and critical points" in Supporting Information S1).
Using Peter Lake as an example thermocline depth decreases as DOC and P load increase ( Fig. 2A). Thermocline depths are within the range observed for Peter Lake (Supporting Information Fig. S2). Critical P load increases slightly with DOC for constant zooplankton biomass (Fig. 2B). In contrast the critical P load increases substantially with zooplankton biomass for constant DOC (Fig. 2B). P load rate (mg P m À2 d À1 ) and (B) critical P load rate vs. DOC and zooplankton biomass (mg C m À3 ) for the area of Peter Lake.

Multivariate time series analyses
Stability and loss of stability were estimated daily from the parameter matrix of multivariate time series models (see section "Multivariate time series analysis to assess stability" in Supporting Information S1). Modulus of the maximum eigenvalue indicates stability if the eigenvalue is less than 1, or loss of stability when the eigenvalue exceeds 1 (Fig. 3). The unenriched reference lake, Paul Lake, had eigenvalues below 1 and appeared to be stable throughout 2015 (Fig. 3A). Peter Lake (Fig. 3B) was destabilized on day 175, with Cyanobacterial blooms apparent by day 170 (Pace et al. 2017;Wilkinson et al. 2018). Although Tuesday Lake received an identical nutrient load, the first indication of instability occurred 44 d later on day 219 (Fig. 3C). Cyanophytes dominated during the blooms in both lakes in 2015 (Wilkinson et al. 2018).

Stochastic resilience
Stochastic indicators of resilience show no consistent differences between Peter and Tuesday lakes (Fig. 4). Resilience to enrichment is illustrated by mean exit time and median survival time of the low-pigment state (Fig. 4A,C). For the low pigment basin, Peter Lake has slightly shorter mean exit time but slightly longer survival time. This difference may reflect the greater influence of extreme values on the mean than the median. Resilience of the Cyanobacterial bloom is illustrated by mean exit time and median survival time of the highpigment state (Fig. 4B,D). Peter Lake has shorter mean exit time and slightly shorter median survival time, suggesting that the bloom of Peter Lake is less resilient than that of Tuesday Lake. The relatively pronounced shift of the mean suggests a stronger influence of extreme fluctuations in Peter Lake.

Discussion
The ecosystem model fitted to the data tracks the delayed bloom in Tuesday Lake (Fig. 1). Analysis of the fitted model suggests that the critical P load, an indicator of resilience, rises steeply with zooplankton biomass and weakly with DOC (Fig. 2B).
The multivariate time series analysis shows that, under constant and equal rates of enrichment, Peter Lake was destabilized 44 d before Tuesday Lake (Fig. 1). This contrast shows greater resilience of Tuesday Lake to enrichment. However, the lakes differed in DOC, thermocline depth, and zooplankton biomass and perhaps other dimensions. Given these differences, we cannot conclude that higher DOC concentrations caused the higher resilience of Tuesday Lake, but it is a possible mechanism for the difference in resilience.
The multivariate time series analysis (Fig. 3) and the critical P load (Fig. 2) provide point estimates of resilience. These indicators show declining resilience and mark the time of critical transition but do not account for the shape of the stability basin or the chance that a random shock may knock the ecosystem over the threshold (Scheffer et al. 2015). Stochastic indicators of resilience use the random fluctuations of the ecosystem to sample the stability basins and thereby reconstruct their shape from a large sample of data . This more complete analysis revealed a more complex picture (Fig. 4). Differences in resilience to enrichment and resilience of the Cyanobacterial bloom are sensitive to random events, and differences between the lakes are smaller than their ranges of variability.
Due to the large number of observations required to compute the stochastic indicators of resilience, we used the highfrequency sensor data from all 3 yr of nutrient enrichments. During 2013-2015 Tuesday Lake had higher DOC but zooplankton biomass, thermocline depth, nutrient enrichment rates and weather fluctuated in each lake each year (Wilkinson et al. 2018;Pace et al. 2019). This variation potentially contributed to overlap in stochastic resilience of the lakes. More frequent pigment measurements (e.g., each minute) could have allowed calculations of stochastic resilience during 2015 alone , providing a comparison of the stability landscapes under constant enrichment. Lakes are exposed to multiple, interacting stressors including eutrophication and brownification (Leech et al. 2018). Using a set of enrichment experiments in lakes with contrasting DOC and other factors, we showed that DOC may affect influence resilience but the effect is complicated by thermocline depth, grazer biomass, and randomness. Further whole-lake experiments are needed to assess the effects of DOC on resilience to enrichment. High-frequency datasets can be collected from a greater diversity of lakes, leading to comparison of stability landscapes under different conditions. A direct manipulation of DOC while collecting high-frequency data before, during, and after the change in DOC is needed as a direct test of the effect of DOC on resilience of phytoplankton. Such studies would lead to a broader understanding of the factors that affect lake ecosystem resilience to enrichment.