Biotic and thermal drivers alter zooplankton phenology in western Lake Erie

Environmental change, particularly warming and eutrophication, can alter phenology in aquatic systems. Understanding controls on zooplankton phenology is important due to their central role in food webs. While patterns in zooplankton phenology have been well documented, we lack an understanding of how abiotic and biotic drivers influence lake zooplankton phenology during the summer. We examined the phenology of four common zooplankton taxa (Daphnia retrocurva, Skistodiaptomus oregonensis, Mesocyclops spp., Dreissenid veligers) in western Lake Erie during 1995–2022, a period with increasing eutrophication and Bythotrephes longimanus biomass. Many phenology metrics varied by 3 weeks or more from 1995 to 2022. The dominant controls of phenology were temperature and biotic factors, especially grazer‐defended phytoplankton (cyanobacteria and other colonial or filamentous taxa) and the invasive predator B. longimanus, which frequently interacted. Our results show that aspects of environmental change interact to shape zooplankton phenology, which can influence phytoplankton biomass and energy flow to higher trophic levels.

Climate warming is altering the timing of physicochemical and biological events in aquatic systems (i.e., Phenology; Ji et al. 2010;Woods et al. 2022).Yet, organismal phenology can also be influenced by biological interactions.For example, zooplankton phenology is shaped by changes in the availability of high-quality food and potential predation (Winder and Schindler 2004;Feuchtmayr et al. 2010;Uriarte et al. 2021).Predicting how zooplankton phenology responds to environmental change is challenging because zooplankton occupy a central food web position between autotrophs and predators, which often exhibit strong and potentially variable seasonality.However, understanding patterns and controls on zooplankton phenology is critical due to the important role zooplankton play in phytoplankton control (Ger et al. 2016) and fish recruitment (Dettmers et al. 2003;May et al. 2021).
Zooplankton phenology is, ultimately, shaped by factors influencing demographic parameters (e.g., growth, fecundity, and mortality) such as temperature, predation, and food availability (He et al. 2022).Water temperature shapes zooplankton emergence from diapause (Jones and Gilbert 2016), growth and reproduction rates (Vijverberg 1980;Gillooly 2000), and prey and predator population demographic rates (Gillooly et al. 2002;Savage et al. 2004).Predation increases zooplankton mortality and thus alters multiple aspects of biomass dynamics and likely phenology.The timing and magnitude of phytoplankton blooms (i.e., food availability) can shape zooplankton phenology during spring (Winder and Schindler 2004;Feuchtmayr et al. 2010;Pothoven and Vanderploeg 2021).During the summer, predation and food quality may be a more important determinant of zooplankton dynamics than temperature (Sommer et al. 2012), but phenological patterns during this period have received less attention in lakes (Vadadi-Fülöp and Hufnagel 2014).
During summer periods, eutrophication often influences zooplankton biomass dynamics; however, isolating the mechanism can be challenging due to covariation between planktivores and blooms of colonial or filamentous phytoplankton that have grazer deterrents (Jeppesen et al. 2005;Ger et al. 2014).In response to phytoplankton blooms, selective grazers like calanoid copepods may reach peak biomass later because they can avoid phytoplankton with grazer deterrents (DeMott and Moxter 1991;Bundy et al. 2005;Moe et al. 2022) while omnivorous copepods may benefit from an increase in microzooplankton (Pace 1986;Bundy et al. 2005).By contrast, generalist grazers, like Daphniids, may reach peak biomass earlier because poor quality diets (i.e., cyanobacteria) and predation depress population growth rates (Moe et al. 2022).Given the central role zooplankton play in lake food webs, we need a better understanding of (1) how much zooplankton phenology varies during summer periods and (2) the factors shaping phenology during this period (Vadadi-Fülöp and Hufnagel 2014).
To address these questions, we characterized the seasonality and phenology of four zooplankton taxa in the western basin of Lake Erie during 1995-2022.Our focal taxa were common and represent functionally and taxonomically different groups.First, we characterized each taxon's biomass seasonality and asked whether that varied in space and time.We hypothesized that seasonality would vary interannually and spatially due to environmental change (see below).Second, we calculated four phenological metrics: the day of year of peak biomass (DOY PB ) for taxa with changing seasonality as well as the start (DOY start ), middle (DOY middle ), and end (DOY end ) of the season.Then, we examined how those phenological metrics were related to water temperature, eutrophication, food availability, and predation.We hypothesized that during the summer, biological interactions (i.e., phytoplankton with grazer deterrents [PP GD ] and planktivory) would be the most important drivers of zooplankton phenology.

Study site and focal taxa
The western basin of Lake Erie has been altered by environmental change, including eutrophication and cyanobacteria blooms (Watson et al. 2016), declining ice cover (Leshkevich et al. 2018), and invasive species including Dreissenid mussels (Karatayev et al. 2014) and the predatory cladoceran Bythotrephes longimanus (hereafter: Bythotrephes) (Barbiero and Tuchman 2004).Air temperature also increased in this region (Dobiesz and Lester 2009).

Data sources
Zooplankton, phytoplankton, and temperature data were provided by the Lake Erie Plankton Abundance Study, which has monitored plankton and physicochemical variables biweekly at eight sites in the western basin of Lake Erie from May to September since 1995 (Supporting Information Fig. S2).Sampling and enumeration methods are described in detail by O'Donnell et al. (2023a;2023b) (summary provided in Supporting Information Appendix S1).Retrocurva and Oregonensis were identified as species, while Mesocyclops and Dreissenid veligers were identified as genera.Copepod immatures were not identified as species or included in Oregonensis or Mesocyclops biomass estimates.Since Oregonensis is rare at western sites, we only used data from eastern sites (sites 16, 14, 8, and 3).
To better understand the ecological and physicochemical factors shaping zooplankton phenology, we used time series representing patterns in phytoplankton biomass with and without grazer deterrents (PP GD and PP NGD , respectively), invertebrate and vertebrate planktivore biomass or density, and water temperature.Eutrophication at these sites is reflected in PP GD , which includes large, filamentous, or colonial taxa (O'Donnell et al. 2023a) and was, during 2008-2022, dominated by Aulacoseira in May-June and cyanobacteria (especially Microcystis spp.) in July-September (Supporting Information Fig. S3).Bythotrephes data were only available from 1995 to 2022 for two sites (27 and 37), so the average was used for the remaining sites.Data sources and postprocessing are described in the Supporting Information Appendix S1.All data and code required to reproduce this analysis are provided in Hood and Bailey (2023).

Statistical analyses
Our approach for characterizing trends in the predictor time series, biomass seasonality for the focal taxa, phenological metrics, and determining the most likely predictors of those phenological metrics is fully described in the Supporting Information Appendix S1.Briefly, to examine the long-term trends in the predictor variables, we used biased random walk models, which evaluate trends while accounting for interannual autocorrelation.We used model selection to determine whether interannual patterns in predictors were best characterized by a basin-wide trend or site-specific trends.Then, to characterize seasonal patterns of our focal taxa and determine whether those patterns varied in space or changed during 1995-2022, we used Bayesian generalized additive mixed hurdle models (GAMM hurdle models).These models are robust to zero-inflated data and combine a presence-only model with a presence-absence model.
For each taxon, we quantified four metrics of zooplankton phenology: DOY PB , DOY start , DOY middle , and DOY end .We estimated DOY PB from the GAMM hurdle models while DOY start , DOY middle , and DOY end were quantified as the 15 th , 50 th , and 85 th percentile of biomass accumulation (Greve et al. 2005), which was predicted using locally weight polynomial regression models (LOESS).We characterized interannual patterns in these phenological metrics with GAMM models and used model selection to determine how patterns differed among sites.Then, to better understand how population growth and mortality contributed to phenology, we used linear models to examine how zooplankton biomass varied with each phenology metric.
Finally, we used generalized linear mixed models and model selection to identify the factors that shaped zooplankton phenology.We evaluated four common predictor variables: water temperature, phytoplankton with grazer deterrents (PP GD ), Leptodora kindtii biomass, Bythotrephes biomass, and planktivorous fish density.To characterize the biomass of preferred food we used phytoplankton without grazer deterrents (PP NGD ) for Retrocurva and Oregonensis (recognizing that Oregonensis also consumes microzooplankton), phytoplankton < 10 μm (PP L10 ) for veligers, and zooplankton biomass (ZP BM ), excluding predacious cladocerans, for Mesocyclops (Williamson 1986).Mesocyclops consumes L. kindtii, but we used L. kindtii and ZP BM as separate predictors.All predictor variables were log e -transformed (excluding temperature), centered, and scaled to their standard deviation within each month, transforming these variables into scaled monthly anomalies.We paired each phenology metric with the mean of predictor values during the prior 30 d (Uriarte et al. 2021).

Interannual trends in predictors
We observed no interannual trend in water temperature, PP NGD , and PP L10 ; whereas PP GD , ZP BM , Bythotrephes, L. kindtii, and planktivorous fish showed a trend in at least 1 month (Fig. 1; Supporting Information Fig. S4).Interannual patterns in water temperature, PP NGD , and PP L10 were similar among sites (a basin-wide pattern; Supporting Information Table S1).During 1995-2022, PP GD showed a basin-wide pattern in June, July, and September and increased at one site, at least in all months except July.Interannual patterns in ZP BM differed among sites and increased in May (two sites), August (six sites), and September (six sites).Bythotrephes showed a basin-wide pattern during all months except June and increased during August and September.L. kindtii also showed a basin-wide pattern during all months except September and increased in May and decreased at one site in September.Planktivorous fish showed a basin-wide pattern during all months but August and increased in July and decreased at one site in August.

Daphnia retrocurva
Retrocurva exhibited an unimodal biomass seasonality, which did not vary among sites but did vary interannually (Fig. 2; Supporting Information Tables S2-S4).While the seasonality of Retrocurva did not vary among sites, peak biomass was lower at sites 27 and 29 (Supporting Information Fig. S5).
In general, Retrocurva phenology advanced (moved earlier) during 1995-2022 (Fig. 4 The line depicts the biased random walk model prediction (ribbon is AE 1 standard error).We used model selection to determine if interannual patterns were similar among sites (all; black lines) or unique to each site (colored by site).Thick lines indicate that there was a significant temporal trend (i.e., 95% confidence interval for bias term did not overlap zero), while thin lines indicated a lack of a temporal trend.Biased random walk model bias terms (AE 95% CI) are provided in Supporting Information Fig. S4.
during 2000-2014 (Supporting Information Fig. S6).By contrast, Retrocurva biomass decreased with DOY middle and DOY end , suggesting that these phenological advances were associated with high biomass and potentially faster population growth.
Retrocurva phenology generally advanced when PP GD and Bythotrephes biomass was high and L. kindtii biomass and planktivore density was low (Fig. 5; Supporting Information Tables S5-S9).Specifically, DOY start advanced with PP GD , Bythotrephes, and was influenced by a PP NGD Â L. kindtii interaction; increasing PP NGD biomass advanced DOY start when L. kindtii biomass was low and delayed it when L. kindtii biomass was high.DOY middle was advanced by PP GD and Bythotrephes, while DOY PB was advanced by Bythotrephes and by a PP GD Â planktivore interaction.DOY end was advanced by Bythotrephes and, when L. kindtii biomass was low by PP GD ; when L. kindtii biomass was high, DOY end was delayed by PP GD .
During 1995-2022, Veliger DOY PB did not change while DOY start , DOY middle , and DOY end generally advanced (Fig. 4).Generally, Veliger phenology was advanced by PP GD and predators (Fig. 5; Supporting Information Tables S5, S12-S14).DOY start advanced with PP GD biomass and planktivore density.DOY middle was influenced by PP GD Â Bythotrephes and PP L10 Â L. kindtii interactions.At low levels of invertebrate  predation, DOY middle advanced with PP GD and was delayed by PP L10 , while at high invertebrate predator biomass, DOY middle was not influenced by PP GD and was delayed by PP L10 .DOY end was advanced by increases in Bythotrephes and PP GD biomass.
DOY start and DOY middle varied similarly among sites, while each site showed a unique pattern in DOY end (Fig. 4).During 1995-2022, DOY start was delayed by 2 d, but was highly variable among sites (inner quartile range: 24 d).DOY middle and DOY PB were delayed by 16 and 29 d, respectively.Patterns in DOY end varied among sites; DOY end showed no trend at two sites, while it was delayed by 36-72 d during 1996-2005 at the other two sites.Oregonensis biomass decreased with DOY s- tart and DOY middle , but not DOY end ; biomass was higher during 2003-2022 than 1995-2002 (Supporting Information Fig. S6).
Oregonensis phenology was influenced by temperature, PP GD , and invertebrate predators (Fig. 5; Supporting Information Tables S5, S17-S20).DOY start advanced with warm temperatures and with PP GD when Bythotrephes biomass was low (a PP GD Â Bythotrephes interaction); when Bythotrephes biomass was high, PP GD delayed DOY start .DOY middle was delayed by L. kindtii and was influenced by a PP GD Â Bythotrephes interaction; increasing PP GD advanced DOY middle when Bythotrephes biomass was low and delayed DOY middle when Bythotrephes biomass was high.DOY PB was delayed by PP GD and Bythotrephes.DOY end was influenced by a temperature Â Bythotrephes interaction; warm temperatures advanced DOY end when Bythotrephes biomass was low and delayed DOY end when Bythotrephes biomass was high.
While Mesocyclops DOY PB did not vary during 1995-2022, DOY start , DOY middle , and DOY end varied considerably among sites and through time (Fig. 4).DOY start was delayed by 5 d at four sites while it advanced 17-23 d at the other four sites.DOY middle was delayed by 17 d while DOY end was, depending on the site, delayed by 4-26 d.Mesocyclops biomass decreased with DOY start , DOY middle , and DOY end ; biomass was generally higher during 2000-2022 than 1995-2003 (Supporting Information Fig. S6).
Mesocyclops phenology was related to a combination of temperature and trophic factors, which often interacted (Fig. 5; Supporting Information Tables S5, S22-S25).DOY start advanced with warmer temperatures and responded to a PP GD Â Bythotrephes interaction.When Bythotrephes biomass was low, DOY start advanced with PP GD while it was delayed by PP GD when Bythotrephes biomass was high.DOY middle was mediated by a similar PP GD Â Bythotrephes interaction.DOY end was advanced by PP GD , delayed by planktivores, and influenced by a ZP BM Â Bythotrephes interaction.An increase in ZP BM advanced DOY end when Bythotrephes was low and delayed DOY end when Bythotrephes was high.

Discussion
Zooplankton phenology in western Lake Erie was shaped by multiple aspects of environmental change, which often interacted.Specifically, changes in phenology were related to temperature as well as PP GD (93% of cases), a symptom of recent eutrophication, and the invasive predator Bythotrephes (93% of cases); PP GD and Bythotrephes often interacted (58% of cases).While warming and food availability are wellestablished determinants of zooplankton phenology (Vadadi-Fülöp and Hufnagel 2014; Uriarte et al. 2021), evidence that summertime eutrophication and predation, particularly from invasive species, can influence phenology is limited (Woods et al. 2022).
In the absence of an interaction, PP GD generally advanced zooplankton phenology, which was often associated with high biomass and likely rapid population growth.PP GD was a better predictor than metrics of preferred prey (PP NGD , PP L10 , and ZP BM ); so the cause of this pattern for four very different zooplankton taxa is unclear.PP GD may provide a predator refuge (Briland et al. 2020); however, high PP GD and predator biomass were often associated with a phenological delay, which was frequently associated with low biomass.Alternatively, the influence of PP GD might be attributed to an increase in food quantity.While filamentous and colonial phytoplankton can depress individual growth rates of Daphnia and Dreissenid veliger in lab experiments (DeMott et al. 1991;Boegehold et al. 2019), PP GD could represent an increase in the abundance of poor-quality food that increases population growth rates.Oregonensis and Mesocyclops selectively feed on motile prey and thus can avoid filamentous and colonial taxa (DeMott and Moxter 1991; Bundy et al. 2005) while consuming small zooplankton which may increase during cyanobacteria blooms (Pace 1986) and may not be wellcharacterized by prey metrics.Food availability has been associated with advancing phenology in aquatic ecosystems (Greve et al. 2001;Pothoven and Vanderploeg 2021;Uriarte et al. 2021); however, our understanding of how summertime eutrophication influences phenology is more limited (Vadadi-Fülöp and Hufnagel 2014).
While understanding of how predation influences phenology is limited (but see Dur et al. 2022;Manca et al. 2007) our results illustrate its potential importance.Bythotrephes, which became common post-2014, frequently interacted with PP GD to advance phenology when PP GD was low and delay phenology when PP GD was high.Trophic interactions between each focal taxa and Bythotrephes differ; Bythotrephes preys intensely on Retrocurva and veligers, but likely not on Oregonensis (Schulz and Yurista 1999; Barbiero and Tuchman 2004;Berges et al. 2020).While Bythotrephes does not consume mature Mesocyclops, Mesocyclops density often declines following Bythotrephes invasion either because Bythotrephes consumes immature Mesocyclops or competes with adults for prey (Yan et al. 2001;Barbiero and Tuchman 2004).Bythotrephes is rare before July; thus, Retrocurva and Veligers (at western sites), can escape intense predation through an advance in phenology.Our study contributes to the limited understanding of how invasive species influence prey phenology.
Temperature played an important role in shaping DOY start and DOY end .Positive temperature anomalies advanced DOY s- tart for Mesocyclops and Oregonensis while it delayed DOY end for Oregonensis.Directionally similar temperature coefficients were also included in Veliger and Retrocurva models with ΔAIC < 5. Similarly, in marine systems, warming often advances the timing of spring zooplankton and delays the timing of fall zooplankton (Mackas et al. 2012).
Our results are based upon analyses of monitoring data which play an important role in tracking ecosystem state and determining how environmental change influences natural ecosystems.However, experiments and examination of demographic rates, is required to identify the underlying mechanisms shaping the phenological patterns we report.Furthermore, since zooplankton monitoring in temperate lakes is commonly limited to $ April-October, we need a better understanding of how winter conditions contribute to plankton succession and phenology (Ozersky et al. 2021).
Zooplankton phenology, particularly characteristics associated with biomass dynamics, is frequently highly variable in freshwater and marine systems.Supporting previous marine and freshwater research (e.g., Feuchtmayr et al. 2010;Uriarte et al. 2021;Dur et al. 2022;Moe et al. 2022), we show that zooplankton phenology is shaped by temperature as well as the direct and indirect effects of eutrophication and predation.While thermal conditions likely control the boundaries of phenology, our work contributes to a growing understanding of the importance of trophic factors in shaping zooplankton phenology.The complex interplay between these variables helps clarify why zooplankton phenology is so variable and illustrates the potential for phenology to shape seasonal succession, top-down control of phytoplankton, and energy flow to higher trophic levels, particularly important fisheries.
Fig. 1.Timeseries of centered, scaled, and logged (log e ) water temperature (Temp.;C), phytoplankton without grazer defenses (PP NGD ; μg L À1 ),phytoplankton less than 10 μm (PP L10 ; μg L À1 ), phytoplankton with grazer defenses (PP GD ; μg L À1 ), zooplankton biomass (ZP BM ; μg L À1 ), Leptodora kindtii biomass (μg L À1 ), Bythotrephes longimanus biomass (μg L À1 ), and planktivorous fish density (fish; catch hectare À1 ) in May-September during 1995-2022.The line depicts the biased random walk model prediction (ribbon is AE 1 standard error).We used model selection to determine if interannual patterns were similar among sites (all; black lines) or unique to each site (colored by site).Thick lines indicate that there was a significant temporal trend (i.e., 95% confidence interval for bias term did not overlap zero), while thin lines indicated a lack of a temporal trend.Biased random walk model bias terms (AE 95% CI) are provided in Supporting Information Fig.S4.
DOY start advanced 24 d at all sites during 1995-2015 and then was delayed by 9 d during 2016-2022.DOY middle advanced 8 and 33 d at sites with unimodal and bimodal seasonality, respectively.DOY end was more variable among sites; it advanced 26 d at four sites, 24 d at two sites during 2005-2022, and at the remaining two sites was delayed by 19 d during 1996-2012 and then advanced through 2022.Like Retrocurva, Veliger biomass decreased with DOY start , DOY m- iddle , and DOY end and was higher during 2003-2022 than 1995-2002 (Supporting Information Fig. S6).

Fig. 4 .
Fig. 4. Generalized additive mixed model estimates of interannual patterns in the start (DOY start ), middle (DOY middle ), and end (DOY end ) of the season for Daphnia retrocurva, Dreissena spp.veligers, Skistodiaptomus oregonensis, and Mesocyclops spp.Model selection was used to determine whether splines for the year were similar across all sites (all; black line) or were unique to each site or group of sites.Final models are provided in Supporting Information Table S5 and interannual patterns in phenology metrics for each taxa and site are shown in Supporting Information Fig. S7.

Fig. 5 .
Fig. 5. Generalized linear model estimates of how predictors (centered and scaled to standard deviation) influenced the day of year (centered) of peak biomass (DOY PB ) as well as the start (DOY start ), middle (DOY middle ), and end (DOY end ) of the season for Daphnia retrocurva, Dreissena spp.veligers, Skistodiaptomus oregonensis, and Mesocyclops spp.In each facet, we provide the predictor terms in each model, the level of significance (" " = p > 0.1; "Á" = 0.05 > p < 0.1; "*" = 0.01 > p < 0.05; "**" = 0.001 > p < 0.01; "***" = 0.001 < p), and the pseudo R 2 (bottom right-hand corner).Interactions are colored based on temperature or prey and are shown for low and high levels of predator values.Predictor values are the average of monthly anomalies (centered and scaled for each month) during the month prior to the phenology metric.