Detecting changes in statistical indicators of resilience prior to algal blooms in shallow eutrophic lakes

. Algal blooms in lakes and reservoirs can be considered regime shifts from a clear-water to algae-dominated state that often occurs abruptly. Under experimental conditions, these regime shifts have been predicted from rises in variance and autocorrelation (generic resilience indicators) of state variables monitored at a high frequency. The goal of this study was to evaluate the behavior of resilience indicators prior to a critical transition in lakes that naturally experience algal blooms. Ambient lake conditions provide several potential hurdles that could inhibit the detection of meaningful changes in resilience indicators prior to a critical transition such as stochastic nutrient loading, spatial complexity, and decreased resilience due to higher baseline nutrient concentrations. We compiled ﬁ ve lake-years of high-frequency monitoring of chlorophyll a , phycocyanin, dissolved oxygen, and pH from four hypereutrophic lakes. Despite the fac-tors that might hinder detecting statistical indicators of changing resilience in hypereutrophic ecosystems, we found that a rise in resilience indicators did occur prior to a critical transition in three out of four possible lake-years, with rise beginning between 5 and 33 d prior. In one lake-year, a critical transition occurred soon after the monitoring began, preventing detection of rising variance or autocorrelation signals which are calculated using a 21-d rolling window. These results add to the growing body of evidence that rises in resilience indicators can be detected in ecosystems prior to a regime shift if monitoring programs are properly designed to capture the dynamics; however, continued research is needed to better understand the conditions under which resilience indicators may be useful as an early warning detection tool for lake management.


INTRODUCTION
Freshwater algal blooms can be harmful to aquatic ecosystems and ecosystem services and may be increasing in frequency, scope, duration, and intensity across the globe (O'Neil et al. 2012, Glibert 2017, Ho et al. 2019. Aquatic ecosystems can be negatively affected by algal blooms resulting in a loss of biodiversity, development of anoxic conditions, and the production of toxins by organisms such as cyanobacteria, which dominate many blooms (Landsberg 2002). Human exposure to the toxins produced by cyanobacteria can result in nausea, vomiting, pulmonary consolidation, or liver damage (Chorus et al. 2000, Codd et al. 2005, Corbel et al. 2014. Furthermore, blooms have economic consequences due to the loss of ecosystem services such as drinking water provisioning and reduced revenue from recreational services (Dodds et al. 2009, Harris andGraham 2017). In short, algal blooms can endanger human health, ecosystem health, and ecosystem services. As such, being able to predict when blooms are going to occur is valuable for preventing or mitigating their negative effects.
Algal blooms can seem unpredictable as they appear suddenly with only a short period of time between low phytoplankton biomass and nuisance levels. The abrupt transition from a clearwater state to an algal-dominated state in lakes is a regime shift (Scheffer et al. 2009). Specifically, a bloom is often observed to be a transition from small cycles in phytoplankton biomass, or a stable point, to large cycles in biomass (Batt et al. 2013). This abrupt transition from a stable point to cycles, a Hopf bifurcation, is preceded in many systems by a rise in variance and a rise in autocorrelation of the state variables (Biggs et al. 2009). The rise in variance and autocorrelation of state variables over time and in space is an indicator of changing resilience (Scheffer et al. 2009, Carpenter et al. 2011 and can be combined with regular monitoring and used as early warning indicators of an impending regime shifts (Carpenter et al. 2011, Dakos et al. 2012, Gsell et al. 2016, Butitta et al. 2017, Nijp et al. 2019.
A rise in the resilience indicators, both over time and across space over time, has been detected prior to the regime shift to algal dominance in experimental whole-lake manipulations (Butitta et al. 2017. In these experiments, cyanobacteria blooms were induced through continual nutrient additions to oligotrophic lakes. Prior to the critical transition, the resilience indicators increased, providing an early warning of the impending shift . While these experiments provided important proof-of-concept that resilience indicators can be detected and were predictive of an oncoming regime shift, questions remain regarding the effectiveness of this approach for predicting algal blooms in a nonexperimental context. These experimental manipulations purposefully induced the regime shift to the algal-dominated state with regular additions of known quantities of nitrogen and phosphorus. However, in non-experimental ecosystems, nutrient inputs are episodic in time and space, which might interfere with detecting dynamics in the resilience indicators (Carpenter and Brock 2006). What's more, the experimental lakes were oligotrophic prior to the nutrient additions and likely had a higher resilience at the beginning of the experiment compared to a nutrient-rich lake. The increased resilience in an oligotrophic system could come from mechanisms such as lower internal phosphorus loading (Orihel et al. 2017) and large-bodied zooplankton grazers. A higher initial resilience could increase the likelihood of observing rises in autocorrelation and variance prior to the critical transition because it lengthens the period of time to detect these dynamics.
Although the current evidence suggests that observing the dynamics of resilience indicators may be a useful tool for predicting algal blooms , the conditions under which the rise in variance and autocorrelation can be detected in a non-experimental setting remain unclear. Careful evaluation of the efficacy of resilience indicators as a predictive tool in non-experimental lakes is needed before further development for management purposes. We collected five lake-years of high-frequency time series from ecosystems with a known history of seasonal algal blooms to (1) determine if a critical transition occurs during a bloom in each case, (2) evaluate if rises in variance and autocorrelation could be detected prior to a critical transition from a clear-water state to an algaldominated state, and if so, (3) to characterize the conditions in which these signals are detectable including the length and frequency of data collection and the severity of the bloom. From this post hoc analysis, we make recommendations for future algal bloom and regime shift monitoring and for the detection of resilience indicators in lakes.

METHODS
We collected high-frequency time series during the growing seasons of 2014, 2015, and 2018 in four shallow lakes that experience seasonal algal blooms for a total of five lake-years of data. The study lakes spanned a range of surface areas and v www.esajournals.org depths, and they are all hypereutrophic largely due to the prevalence of agricultural land use in the watersheds (Arbuckle and Downing 2001;  Table 1). All the lakes in our study were also on the Clean Water Act 303(d) list of impaired waters at the time of sampling for excessive algal growth or bacterial impairment.

Data collection
High-frequency time series were collected using YSI EXO2 (monitoring in 2014 and 2015) and YSI EXO3 (monitoring in 2018) multiparameter sondes (Yellow Springs Instrument, Yellow Springs, Ohio, USA) equipped with sensors for Total Algae, which measures both chlorophyll a and phycocyanin, dissolved oxygen percent saturation (DO), pH, and temperature (Ortiz et al. 2019). While chlorophyll a and phycocyanin concentrations are direct measures of the ecosystem state with regard to an algal bloom, DO saturation is also a state variable that responds to the shift to algal dominance due to changes in primary production (Batt et al. 2013, Sirota et al. 2013. Another variable that responds to primary production is pH, particularly in poorly buffered systems. As dissolved carbon dioxide is removed from the water column by photosynthesis, the concentration of carbonic acid decreases, thereby causing pH to rise. The sondes were deployed approximately 0.5 m below the surface at the historical deep site in each lake during the ice-free period. The sensors were calibrated every two weeks from 2014 to 2015 and regularly in 2018 when indicated by the KorEXO software which monitors instrument drift during deployment. The length of the collected time series ranged from 99 to 153 d of continuous high-frequency (15 min) data collection (Table 1). From the high-frequency data, daily mean values for each state variable-chlorophyll a, phycocyanin, DO, and pH-were calculated and used in the resilience indicator analysis.
Synoptic water quality samples, which included total nitrogen and phosphorus and the algal pigments, chlorophyll a and phycocyanin, were taken regularly during the sonde deployment periods (Table 1). Total nitrogen was measured as the sum of total Kjeldahl nitrogen (US EPA 1993c, method 351.2 v2) and nitrate + nitrite measured using the cadmium reduction method (US EPA 1993a, method 4500-NO 3 -F). Total phosphorus was measured using the ascorbic acid method (US EPA 1993b, method 365.1 v2).

Statistical analysis
In order to evaluate the dynamics of the resilience indicators prior to a critical transition, we first needed to identify if and when a critical transition occurred in each lake-year. Critical transitions are generally associated with algal blooms, while there being no standard definition for what constitutes an algal bloom state vs. a non-bloom state in lakes. Studies that have provided a definition based on a threshold algal biomass have used long-term data with frequent observations of the water body in a non-bloom state (Smayda 1997) which is not available for the lakes in this study. Therefore, in order to determine if a critical transition occurred, we adapted the online dynamic linear modeling analysis presented in (Taranu et al. 2018; and the code therein) to estimate the eigenvalue of each pigment time series in each lakeyear. Time-varying autoregressive models were fit to the time series of the pigments using a discount factor of 0.9, and the optimal order was Table 1. Summary of water quality parameters and lake characteristics for each lake-year including mean total nitrogen (TN) and total phosphorous (TP), and agricultural land use (percentage watershed ag) in their respective HUC12 watershed. v www.esajournals.org selected using AIC. When a critical transition occurs, the eigenvalues cross one from below, consistent with a loss of resilience (Biggs et al. 2009, Scheffer et al. 2009. If the eigenvalue does not cross one from below, a critical transition from the clear-water state to the algaldominated state is unlikely to have occurred. We evaluated the behavior of the resilience indicators in relation to the day of year (DOY) that the eigenvalue crossed one from below in each pigment time series as well as the behavior if there was no evidence of a critical transition. The resilience indicators were calculated using a rolling window analysis of log (x) transformed daily averages of pigment concentrations and DO, leaving pH untransformed. Note that the overall pattern of the dynamics in resilience indicators calculated using untransformed data was similar to the patterns in the transformed data presented here (Appendix S1). While there are many candidate resilience indicators (Dakos et al. 2012), in this study we examined lag-1 autocorrelation (hereafter, AC) and a measure of variance, standard deviation (SD), as these statistics have successfully detected the regime shift in question in experimental manipulations (Batt et al. 2013. In this analysis, we used a 21-day rolling window which strikes a balance between detecting signals in the resilience indicators due to an impending regime shift, the timing of the critical transition in relation to the length of the time series, while minimizing false positives in the resilience indicator behavior. Rolling window analysis of AC and SD of the chlorophyll a, phycocyanin, DO and pH time series, and sensitivity analyses were done in R (version 3.4.4; R Core Team 2018) with the earlywarning package (Dakos et al. 2012). Rolling window length sensitivity was tested for a range of 10-21 d. Here, we present the results for a 21-d window length analysis of the entire dataset (see Appendix S1 for 15 day (d) rolling window results).
The rolling window time series of the AC and SD resilience indicators were examined for increases in the indicator values prior to the critical transition (if a transition occurred). We considered a rise in the resilience indicator as beginning when there were five consecutive days of increase in the rolling window statistic. We selected the five-day minimum as a balance between detecting ephemeral rises that would result in a false positive and being able to detect a signal prior to the transition. Additionally, in the experimental manipulations, changes in SD and AC relative to a non-blooming reference ecosystem were approximately five days or longer prior to a bloom ). The critical transition in question is the shift in phytoplankton biomass from a stable point to stable cycles. As such, the pigment time series which are a proxy for phytoplankton biomass were evaluated for the evidence of a critical transition, and the date of that critical transition was applied to the pH and DO time series which are influenced by phytoplankton biomass, but not a direct measure of that biomass. All four of the state variables were examined for rises in the resilience indicators prior to the date of the critical transition (i.e., early warning indicators). Additionally, the dynamics of the resilience indicators were examined for false positives (a five-day rise without evidence of a critical transition).

RESULTS
All five lake-years had high pigment concentrations, indicative of blooms, at some point in the time series. Of the five lake-years, Black Hawk 2015 and Green Valley 2014 had evidence of a critical transition in both pigments based on the eigenvalues time series crossing one from below ( Fig. 1). Black Hawk Lake had the largest bloom with a maximum daily average concentration of 169 lg/L (Fig. 2). The eigenvalues remained at or near 1 early in the season (DOY 156), but clearly crossed one from below on DOY 182 (Fig. 1). Both AC and SD of chlorophyll a and phycocyanin began and continued to rise for approximately two weeks prior to the critical transition on DOY 182 which coincided with an increase in the algal biomass (Table 2; Fig. 2). While there was also a rise in SD of the DO and pH state variables prior to the critical transition, AC did not respond. There was evidence of a second critical transition in phycocyanin on DOY 237 in Black Hawk 2015 which was preceded by a short, but sharp, rise in SD and AC of phycocyanin prior to the transition ( Table 2).
The critical transitions in Green Valley 2014 occurred as the bloom collapsed (Fig. 3). In this v www.esajournals.org instance, only the SD resilience metric rose prior to the transition, and only for the phycocyanin, DO, and pH variables (Table 2). During the following year in the same lake, Green Valley 2015, the critical transition also occurred as chlorophyll a concentrations were falling and were preceded by a rise in AC and SD of chlorophyll a for a week ( Table 2). The SD of DO also rose consistently for 33 d prior to the transition, although no other resilience metrics or state variables responded to the critical transition. Interestingly, while chlorophyll a was highest in the spring and continued to decline over the course of the summer in both lake-years, the phycocyanin concentrations in Green Valley Lake had a bloom and bust pattern in both years, mirrored by the dissolved oxygen and pH values (Figs. 3, 4).
The critical transition in the chlorophyll a time series in South Twin Lake in 2018 occurred on DOY 171 as pigment concentrations were rising (Table 2; Fig. 5). This transition occurred within 21 d of the sensors being deployed in that lakeyear, and therefore, the dynamics of the resilience metrics are unable to be examined for this transition. Finally, there was evidence of two critical transitions early in the season in the Swan Lake 2018 chlorophyll a time series as a bloom occurred. The first transition on DOY 156 was too early in the sensor's deployment to observe rolling window metrics; however, the second transition on DOY 170 was late enough to observe the dynamics of the resilience metrics (Fig. 6). Neither of the resilience metrics in any of the state variables rose prior to the critical transition on DOY 170.
Black Hawk Lake in 2015 had the highest percentage of potential resilience indicators respond by rising prior to the critical transition (Table 2), which aligned with a bloom in the lake (Fig. 2). Interestingly, Black Hawk 2015 is also the longest time series (Table 1) and a later bloom, which helped to clearly capture the dynamics of the resilience indicators prior to the transition. Conversely, the early-season blooms and short observation windows in South Twin 2018 and Swan Lake 2018 prevented or hampered detection of rising resilience indicators prior to the critical transition (Table 1; Figs. 5,6). In the two lake-years where the critical transition occurred as the blooms collapsed, Green Valley 2014 and 2015, only 37.5-50% of the potential resilience indicators rose prior to the transition (Table 2; Figs. 3,4).
Across all lake-years, 14 resilience indicators of the possible 30 rose prior to the critical transition ( Table 2) v www.esajournals.org pigments were more responsive prior to a transition than DO and pH. Additionally, detecting a rise in SD prior to a critical transition was more than twice as common as a rise in AC among all lake-years. While the results for the 21-d rolling window are presented here, we also evaluated the dynamics of the resilience indicators for a 15d rolling window length. Decreasing the window length from 21 to 15 d resulted in slightly fewer rises in the resilience indicators being detected; however, the pattern of rising resilience indicator among lakes and among state variables remained the same (Appendix S1).
To assess the accuracy of the rolling window metric being detected prior to the critical transition, we tallied the rates of true and false positives and false negatives among all lake-years. The true positive rate, or the number of times the resilience indicators rose prior to a critical transition, was 46.7%. Conversely, the false negative rate, or the number of times the resilience indicators did not rise prior to a critical transition, was 53.5% (Table 2). A true false positive rate, or the number of times there was a rise in the resilience indicators without an impending critical transition, is difficult to calculate with truncated time v www.esajournals.org series as the dynamics of the resilience indicators could be heralding another impending transition that occurred after monitoring ceased in that lake-year. However, examining the rolling window time series, the AC of DO, and SD of DO and pH both rose steadily after the critical transition in Black Hawk 2015 (Fig. 2) and Green Valley 2014 (Fig. 3). Additionally, there was an increase in the SD of chlorophyll after the transition in Green Valley 2015 (Fig. 4) and sharp spike in chlorophyll SD in Swan Lake 2018 (Fig. 6) shortly after the transition occurred. All of these dynamics may be interpreted as false positives if the information were being examined in real time instead of post hoc.

DISCUSSION
In this study, we analyzed five lake-years of data from shallow, hypereutrophic lakes to assess the dynamics of resilience indicators prior to the critical transition to algal bloom formation or collapse. We found that three out of the four lake-years where rolling window resilience indicators could be observed, expressed a distinct rise in AC and/or SD in at least one state variable between 5 and 33 d prior to the transition. However, the number of state variables that rose prior to the critical transition within a given lake-year, congruent with our hypothesis of their behavior prior to a regime shift, was quite variable. As such, using the dynamics of resilience indicators as early warning indicators of an impending regime shift to algal dominance should be employed cautiously in the types of ecosystems evaluated in this study.
Both time series length and the monitoring scheme played an important role in being able to detect rising dynamics in the resilience indicators prior to critical transitions. For example, Blackhawk Lake in 2015 had the highest number (75%) of resilience indicators with a distinct rise prior to the first critical transition on DOY 182. This lake also had the longest continuous time series at 153 d, and there was a longer period of Note: A dash indicates a rise in resilience indicators could not be evaluated as the transition occurred within 21 d of measurements beginning in that lake-year.
† A critical transition did not occur in this state variable during this lake-year.
v www.esajournals.org observation prior to the critical transition than in other lakes. Under these circumstances, a rise in at least one of the resilience indicators of all four state variables was observed, lending a high amount of confidence that the rise in the indicators was an early warning of the impending critical transition. In contrast, both Swan Lake and South Twin Lake in 2018 were monitored for shorter periods of time with the critical transition occurring early in the monitoring period. This likely prevented observing clear rises in the resilience indicators due to the short observation period.
The algal pigment state variables (chlorophyll a and phycocyanin) had the highest incidence of the resilience indicators rising prior to a critical transition across lake-years. The pigment concentrations are a more direct measure of algal biomass compared to DO or pH, which may account for the responsiveness of these state variables (Batt et al. 2013). In lake-years where a critical transition occurred, the transition was generally in advance of the peak in algal biomass v www.esajournals.org bloom was collapsing, the resilience indicators did rise for some variables prior to the transition, providing evidence that the online dynamic linear modeling can detect when lakes both enter and exit a blooming state and that the resilience indicators are responsive to these transitions, as hypothesized.

Resilience indicators and bloom prediction
Rises in rolling window resilience indicators examined in this study and other statistical moments (Dakos et al. 2012) have been interpreted as early warning indicators of an impending shift to algal dominance in other lakes (Butitta et al. 2017) and for other regime shifts (Carpenter et a. 2011). However, the mixed results from the post hoc analysis in the five lake-years presented here demonstrate that if the dynamics of resilience indicators are going to be used as early warnings of a regime shift and a lake management tool, sensor-based monitoring program needs to be specifically designed for this goal. Most importantly, extended, continuous monitoring is needed, particularly before the onset of the bloom season. Typically, sensor-based lake monitoring occurs during portions of the ice-off period in temperate regions, which does not lend itself to fully capturing all of the dynamics in the ecosystem. There can be substantial portions of the annual primary production occurring outside the traditional monitoring window, such as a spring bloom period or under ice (Hampton et al. 2017). Longer periods of pre-bloom conditions relative to the rolling window length are necessary to capture the dynamics in resilience indicators, particularly if they are to be used in an early warning context. Given these constraints, happenstance high-frequency data collection is unlikely to be suitable for evaluating resilience indicator dynamics to predict critical transitions (Gsell et al. 2016).
Other components of a monitoring program can also affect the efficacy of detecting rises in resilience indicators prior to a critical transition. In our study, the sensors were measuring at a fixed point, approximately 0.5 m below the surface. We used the mean daily pigment values in the analysis, partially to account for any mixing or phytoplankton migration that may have caused variability at a sub-daily timescale (Klemer et al. 1996, Cui et al. 2016). However, phytoplankton biomass may still have been over-or undersampled with fixed depth measurements, particularly if there was a surface scum-forming bloom. The sensors were also positioned at a v www.esajournals.org fixed location over the deepest point in the lake which is a common lake monitoring scheme. However, lakes can be spatially complex and sampling at a fixed location can result in important algal dynamics being missed. Positioning sensors in algal hotspots that consistently bloom first or for the longest period may yield more robust detection of rises in resilience indicators (Dokulil andTeubner 2000, Alexander andImberger 2009). Finally, phytoplankton may not be the only source of primary production in the ecosystem, making the use of DO and pH to detect rises in resilience indicators difficult to interpret. For example, in Swan Lake in 2018, the DO and pH dynamics became decoupled from the pigment dynamics in the later portion of the season, coinciding with the growth of large macrophyte beds in the lake. As such, relying solely on these secondary proxies of primary production may be misleading in a regime shift context.
Overall, when there were rises in resilience indicators prior to the critical transition, they began to rise 5-33 d prior (median = 16 d). In an early warning context, this size of a warning period may allow time for lake managers to act proactively, or at minimum communicate with recreational users about impending hazards, such as potential increase in cyanotoxin concentrations. If resources are available, these warning v www.esajournals.org periods would also allow lake managers to mitigate the bloom's potential impact using techniques such as aerators, algaecides, and harvesting (Lilndenschmidt 1999, Carmichael et al. 2000, Ma and Liu 2002. However, in order to intervene, decision makers need to have confidence in the signal they are interpreting. In our post hoc analysis, it is easier to determine a meaningful rise in SD or AC of a state variable than if the analysis were happening in real time (e.g., Pace et al. 2017). While there are tools such as quickest detection (Carpenter et al. 2014) that provide a more objective evaluation of resilience indicator dynamics, the method requires information about the dynamics in both regimes which might not be available in all cases. Further effort is needed to develop methods for decreasing false positives and increasing confidence in resilience indicators as regime shift prediction tools. Additionally, evaluating and comparing the performance of other algal bloom prediction models with resilience indicator dynamics would be helpful for further understanding the performance of algal bloom prediction models.
The evaluation of resilience indicator dynamics under non-experimental conditions in this study has expanded our understanding of the efficacy of generic resilience indicators as early warning indicators of regime shifts under non-experimental conditions. These results add to the growing body of evidence that rises in resilience indicators can be detected in ecosystems prior to a regime shift (Weijerman et al. 2005, Scheffer et al. 2009, Nijp et al. 2019, Schmitt et al. 2019; however, the ability to detect these dynamics is predicated upon properly designed monitoring programs as demonstrated in these hypereutrophic lakes. While this study is among the first to evaluate resilience indicator dynamics prior to an algal bloom in non-experimental conditions, further work is needed across a broad array of trophic states and lake types to better define the circumstances under which rises in resilience indicators precede the transition to an algal bloom state.