Zooplankton recovery from a whole-lake disturbance: Examining Zooplankton recovery from a whole-lake disturbance: Examining roles of abiotic factors, biotic interactions, and traits roles of abiotic factors, biotic interactions, and traits

Community assembly following disturbance is a key process in determining the composition and function of the future community. However, replicated studies of community assembly at whole-ecosystem scales are rare. Here, we describe a series of whole-lake experiments, in which the recovery of zoo-plankton communities was tracked following an ecosystem-scale disturbance, that is, application of the piscicide, rotenone. Using a before-after-control-impact design, 14 lakes in eastern Washington were studied: Seven lakes were treated with rotenone, while seven lakes acted as reference systems. Each lake was monitored up to 6 months before and 1 – 2 years after the rotenone treatments. Zooplankton samples and environmental measurements were collected approximately monthly from each lake. Community responses following disturbance were assessed using metrics of abundance, diversity, and community composition, as well as taxonomic group abundance. Zooplankton recovery was also assessed using species traits related to habitat, feeding mode, trophic level, body size, and life history. In addition to patterns of recovery, potential mechanisms were explored relating to abiotic conditions, biotic interactions, and traits. There were steep declines in the abundance (average across years: 99%) and diversity (average across years: 75%) of the zooplankton community following rotenone treatment. Although abundance had


INTRODUCTION
Disturbances are a key factor in shaping communities and ecosystems. The intensity, frequency, size, and return interval of disturbance regimes impact community and ecosystem structure and function (Turner, 2010;White & Pickett, 1985). Species may become adapted to these events (Dayton, 1971;Gunderson, 2000); however, anthropogenic disturbances may pose novel threats that species are not adapted to. For instance, anthropogenic disturbances can shift communities to new equilibria (Brasher, 2003;Folke et al., 2004). The factors that influence community responses to disturbance are still uncertain, but remain a vital concern for resource managers whose understanding of recovery processes helps to dictate their management strategies (Fraterrigo & Rusak, 2008;Holling & Meffe, 1996).
Recovery to disturbance evolves over time. Immediately following severe disturbances, community assembly occurs, which is dependent on immigration and emergence from dormancy, and can be affected by the severity of the disturbance and system characteristics (Brock et al., 2003;Chase, 2003;Myers et al., 2015). The first species that appear following disturbance tend to be those that are resistant to the disturbance, as well as species that develop quickly, representing high relative fitness (HilleRisLambers et al., 2012;Jentsch & White, 2019). Early colonists may show priority effects, where their early colonization leads to short-term dominance of the community (Frost et al., 2006;Jenkins & Buikema Jr., 1998;Urban & De Meester, 2009). Additional stressors on the community, such as predator introductions, may reduce the importance of priority effects and lead to secondary succession, where the community evolves to consist of slower developing species (DeMott, 1989;Koch, 1974). Secondary succession is dependent upon species interactions, which are affected by food availability, abiotic conditions, and species traits (Frost et al., 2006;Lepori & Hjerdt, 2006). Additionally, legacy effects resulting from historical disturbances can have a lasting influence on community recovery to future disturbances (Johnstone et al., 2016). However, the relative importance of the factors that drive recovery in communities following disturbance is still unclear.
Studies on community recovery from disturbance have limitations. Disturbances are frequently unplanned, and thus, there are rarely data to determine baseline conditions prior to impact. Long-term study sites have demonstrated slow-developing successional processes from unplanned disturbances, but these studies often lack replication and appropriate reference systems (Peterson et al., 2003;Turner et al., 2003), making it difficult to ascertain whether recovery trends are synchronous across space. Scale is also a factor, as small-scale experiments may not accurately capture larger scale processes, upon which management decisions depend (Schindler, 1998). Additionally, many studies rely on coarse structural attributes, such as abundance, richness, or diversity, which may miss certain recovery patterns.
Lakes are a model system for studying recovery from disturbance because of their defined boundaries (Post et al., 2007). Whole-lake manipulations have often, though not always, resulted in clear ecological responses, such as eutrophication following the addition of phosphorous (Schindler, 1974). Zooplankton occupy a central role in lake food webs (Walsh et al., 2016) and are especially useful subjects due to their small size, short life span, fast generation times (<10 days; Gillooly, 2000), ease of collection, and sensitivity to environmental perturbations. Zooplankton can experience both top-down (i.e., predation) and bottom-up (i.e., resource availability) pressures given their position in the middle of lake food webs. Additionally, zooplankton can exhibit synchronous dynamics at a variety of timescales (Vasseur et al., 2014), making them an ideal system to study recovery from disturbance.
Zooplankton functional traits have been described in the literature (Barnett et al., 2007;Hébert et al., 2016), which may provide key insights into how communities recover from perturbations. In particular, Litchman et al. (2013) developed a framework for trait-based approaches to better represent zooplankton functions in communities and ecosystems. Zooplankton body size and motility were identified as key ecological traits that transcend multiple functions and represent important biological trade-offs, and are therefore representative of ecological strategies (Litchman et al., 2013). The application of this framework may help advance the study of community recovery from disturbance by disentangling individual species responses from a more generalizable functional trait response.
Rotenone, a pesticide used to manage nuisance fish populations, is an example of a severe, short-lived disturbance, due to its high toxicity to obligate gill-breathing organisms, as well as crustacean zooplankton (Melaas et al., 2001), and its ability to rapidly breakdown to nontoxic forms (Dalu et al., 2015;Finlayson et al., 2014). Previous studies reviewed by Vinson et al. (2010) on the recovery of zooplankton communities from rotenone treatments have been equivocal, varying in the time to recovery and taxon-specific sensitivity. However, these results were largely based upon short-term, unreplicated studies with no standard definition of recovery. Though previous studies have yielded important insights (Duggan et al., 2015;Melaas et al., 2001), there has been a lack of generalizable trends resulting from the application of this powerful pesticide. Further, there is a clear management need to understand both short-and long-term effects of rotenone on nontarget organisms using well-designed studies (Vinson et al., 2010).
Here, we describe patterns and mechanisms of crustacean zooplankton recovery in replicated whole-lake experiments. Our objectives were to examine overall patterns of zooplankton community structure following rotenone, and examine the factors that influence recovery. We assessed the recovery of rotenone-treated lakes to nontreated lakes by sampling before and after impact, using a before-after-control-impact (BACI) design (Underwood, 1991). BACI studies are well suited to assess community responses to disturbance, as the effect of treatment on the impacted sites is compared to baseline changes in the reference sites during the same time period (Christie et al., 2019). Here, we define recovery as no significant differences in zooplankton community structure, taxonomic groups, and community composition in treatment lakes compared to reference lakes (Xiang et al., 2014). Lastly, we use spatial synchrony, defined as synchronous changes among lakes within a region (Vasseur et al., 2014), as well as within-lake synchrony of zooplankton, to examine the consistency of zooplankton recovery over time.
We hypothesize that there will be three major drivers structuring recovery: functional traits (e.g., life history, morphology, and physiology), abiotic conditions, and biotic interactions. We present four testable predictions that stem from these drivers, either acting alone or potentially interacting (Table 1, references therein), summarized briefly here. (1) We predict that cyclopoid copepods will be the first taxonomic group to recover, as they can undergo dormancy in an advanced stage, may have priority effects advantages, and, due to their mobility, may experience a greater benefit from the loss of planktivorous fishes compared to other taxonomic groups. Although we did not quantify dormancy, crustacean zooplankton exhibit fairly consistent diapausing strategies within taxonomic groups (reviewed in Gyllström & Hansson, 2004). (2) We predict that zooplankton in warmer, more productive lakes will recover the fastest. Some possible mechanisms include greater emergence from diapaused eggs, faster swim speeds (and therefore higher encounter rates, reducing Allee effects), faster developmental rates, and greater food resources. (3) We predict that fish stocking following rotenone treatment will promote species coexistence via predation on the most abundant taxa (i.e., early colonists), and thus, we will observe increased zooplankton diversity following fish stocking. (4) We predict that there will be asynchrony between taxonomic groups within lakes, inferring that competition between groups will affect recovery patterns. Ultimately, our research may inform fish management in lakes and, more broadly, the interplay of disturbance factors, traits, abiotic conditions, and biotic interactions in structuring community recovery from disturbance.

Study sites and design
In our study, the impacts on seven rotenone lakes were compared to seven reference lakes, with sampling occurring both before and after treatments occurred in both sets of lakes. Rotenone was applied to lakes in the fall of either 2014 (n = 4) or 2015 (n = 3). All of the lakes in our study are located in eastern Washington in three ecoregions: (1) Columbia Plateau (Amber, Badger, Dry Falls, Lower Hampton, Rat, Upper Hampton, and Widgeon); (2) Okanogan (Big Twin and Lost); and (3) Canadian Rocky Mountains (Bayley, Browns, Cedar, McDowell, and No Name) (Figure 1). The Columbia Plateau ecoregion is a semiarid landscape with a mix of channeled scabland and coulee areas, sagebrush vegetation, and is highly impacted by large irrigation agriculture, including the vast Columbia Basin Project. The Okanogan and Canadian Rocky Mountain regions are mountainous landscapes of coniferous forests. The climate of these regions consists of hot, dry summers and cold, snowy winters with increasing precipitation at higher elevations.
Lakes were chosen to be as similar as possible across categories (i.e., reference vs. rotenone). The study lakes are relatively small (4.5-100 ha), low elevation (305-1300 m), shallow to moderate depth (6-30 m), and mostly mesotrophic (total phosphorus [TP] from 8-33 μg/L) ( Table 2). All of our treatment lakes, but one, have been previously treated with rotenone within the previous 13 years. There were no significant differences between reference and rotenone lakes for each of these variables prior to treatment using a Welch's t test, which assumes unequal variances between groups (Table 2). However, both categories of lakes span relatively large ranges across variables and include some outliers (e.g., years posttreatment and TP). We took the additional step of comparing lakes across multiple response variables using an analysis of similarity, which analyzes how similar groups of sampling units are. Using the R library vegan (v. 2.5-7), we first performed a Z score transformation of the key lake variables (Table 2) and then analyzed the data using Euclidean distances with n = 999 permutations. There was no significant difference between reference and rotenone lakes (R = 0.078, p = 0.138). Thus, we can cautiously conclude that lake categories are relatively similar and will indicate instances where outliers may be influential.
These lakes are mostly hydrologically disconnected, except for seasonal connections to small ponds and wetlands. Both reference and rotenone lakes were typically stocked annually with fingerlings (5.1-7.6 cm) and/or catchable-sized (27.9-33.0 cm) fish by Washington Department of Fish and Wildlife (WDFW; Appendix S1: Table S1).

Sampling methodology and zooplankton enumeration
Lakes were sampled monthly when conditions permitted from June 2014 to September 2016 (n = 230 total samples; Appendix S1: Table S2). Our sampling partially coincided with required sampling as a part of WDFW's National Pollutant Discharge Elimination System (NPDES) permits conditions for applying rotenone (permit number WA0041009). Winter sampling was inconsistent due to ice conditions and accessibility, which created safety concerns. Forest fires also limited access to two of our lakes (Lost and Rat). Rotenone-treated lakes were sampled more intensively in the weeks before and after rotenone exposure. Rotenone treatments were applied in October 2014 and 2015 using powder and/or liquid forms, at a final concentration that ranged from 2.0 to 4.0 ppm (median = 3.7; Appendix S1: Table S3). We are aware that some lakes have been rotenoned multiple times in the past (WDFW, unpublished), but do not have comprehensive data on prior treatments.  Elgmork (1980), Gerritsen (1980), Naess and Nilssen (1991) Traits (physiology) Motile, raptorial feeders will experience weakened predation risk, increasing feeding efficiency Barnett et al. (2007), Litchman et al. (2013) Biotic Priority effects, where early colonists monopolize resources and exclude later colonizing species Jenkins and Buikema Jr. (1998), Urban and De Meester (2009) 2. Recovery will be fastest in warmer, more productive lakes Abiotic + traits (life history) Greater emergence from diapause at higher temperatures Arnott and Yan (2002) Abiotic + traits (physiology) Temperature increases swim speeds, increasing encounter rates with potential mates; temperature increases development rates Gillooly et al. (2002), Kramer et al. (2011) Abiotic + biotic More food resources available in productive lakes will hasten recovery Lepori and Hjerdt (2006) 3. Fish stocking after rotenone will increase species diversity Biotic Most abundant taxa (early colonists) will be preyed on by fish, allowing other species to establish Tucker and Fukami (2014) 4. Asynchrony between taxonomic groups within lakes Biotic Competition for resources Vasseur et al. (2005) Traits (life history) Faster developing species will appear early in succession, especially spring/early summer, as there will be more food available; slower growing species will outcompete early successional species later in the summer Koch (1974), Sommer et al. (1986), DeMott (1989) Physical and chemical characteristics of lakes were measured on each sampling visit at the site of maximum depth. Different multiparameter meters were used to measure physical and chemical variables throughout the study; however, the same meter was consistently used on each lake, with a few exceptions (Appendix S1: Table S4). Meters were calibrated prior to each monthly sampling event. Data points were removed when anomalous values were detected that were thought to be errors. Temperature, pH, conductivity, and dissolved oxygen were recorded at 1-m intervals. Water clarity was measured at the deep site using a Secchi disk. In July of each sampling year, an integrated water sample of the epilimnion was collected using a 2.5-cm-diameter tube sampler for nutrient analysis. These unfiltered water samples were stored on ice in the field and frozen upon return to the laboratory at the end of the day until later analysis. Samples were analyzed with a Shimadzu UV-1800 Spectrophotometer (Shimadzu, Kyoto, Japan)  methods with a persulfate digestion. Method detection limits for TN and TP were 10 and 2 μg/L, respectively. Vertical plankton samples were taken using an 80-μm mesh, 1.2-m-length net with a 30 cm diameter. An 80-μm mesh net is routinely used in lake monitoring (Yan et al., 1996) and is able to capture even very small crustaceans and nauplii (A. Strecker, unpublished). Each lake was sampled at three sites: shallow, middle (intermediate between the deep and shallow sites), and the deepest spot in the lake, determined by a bathymetric map. This sampling methodology allowed assessment of the full diversity of crustacean zooplankton in the lake, including littoral taxa (Walseng et al., 2006). We chose to focus on crustacean zooplankton (excluding rotifers) for several reasons: (1) to be consistent with the sampling required by the NPDES permit and (2) limitations on gear availability across four WDFW districts (80-μm mesh nets). The shallow site was ≥4 m deep to account for the length of the net. The same locations were sampled at each visit, confirmed with Global Positioning System (Appendix S1: Table S5). Zooplankton were preserved at a final concentration of 70% ethanol.
Identification of zooplankton was conducted with a Leica M165C microscope (Leica Microsystems Inc., Buffalo Grove, IL). To reduce the high number of samples to enumerate (230 lake-dates Â 3 depths = 690 total samples), composite samples were made of each lake by volume weighting the deep, middle, and shallow sites. Here, volume weighting refers to the volume of the water filtered by the plankton net during sampling. In brief, each sample was gently homogenized and poured into a volumetric cylinder. These samples were then topped off with deionized water to a volume representing that samples' contribution to the overall volume of all three samples. For example, if the deep spot = 50 L of water filtered, middle = 30 L, and shallow = 20 L, the representative volumes could be 50, 30, and 20 ml. All three volumetric flasks were then combined into a single composite sample. Subsamples from the composite sample for enumeration were obtained with a plankton splitter.
T A B L E 2 Summary of lake historical, physical, and chemical metrics, with means and SD for reference and rotenone lakes and results of Welch's two-sample t-tests (df = 12 for all categories except year posttreatment, which was df = 11) Lake Years post Abbreviations: cond, specific conductance (water column); pH, surface pH; TN, total nitrogen; TP, total phosphorus; years post, years since last rotenone treatment; Z max , maximum depth.
The enumeration procedure followed Strecker and Arnott (2005), which included counting at least 50 adult individuals of each species, 25 juveniles of each order, until reaching 250 individuals per sample. This protocol is designed to search more of the sample for rare species. Adult individuals were identified to species level when possible, and copepod juveniles (nauplii and copepodids) were identified to order or subclass using taxonomic keys (Haney et al., 2013;Thorp & Covich, 2010). Early developmental stages of cladocerans (e.g., red eye and curled tail spine) were identified simply as juvenile cladocerans (Toyota et al., 2016).
Metrics of zooplankton community structure were calculated for each sample, including abundance, Shannon-Wiener diversity, and abundance of each of the three major taxonomic groups (calanoids, cyclopoids, and cladocerans). Categorizing taxa within these groups was necessary because many species were only found in a few lakes. These metrics were averaged monthly to compare baseline trends in reference lakes to trends in rotenone lakes.

Statistical analysis
We first present our methods for describing zooplankton responses to rotenone and then follow with methods used to test specific predictions (Table 1). We used linear mixed effects (LME) models on community metrics (total density, diversity, and taxonomic group abundance) to test for patterns of recovery from rotenone. LME models are common to BACI experiments with repeated measures data testing the effects of disturbance (Underwood, 1991). LME models are a valid alternative to repeated measures analysis of variance (ANOVA), because they can account for non-normalcy and nonindependence (Harrison et al., 2018).
Lakes that were treated with rotenone in 2014 and 2015 were analyzed in separate models (hereafter referred to as 2014 rotenone and 2015 rotenone), in order to account for the differences in treatment years. Two LME models were run for 2014 rotenone lakes: Year 1 and Year 2 following treatment. The 2014 rotenone Year 1 model contrasted the four rotenone lakes with the seven reference lakes from June 2014 to September 2015, whereas the Year 2 model again used pre-impact data from June to September 2014 contrasted with postimpact data from October 2015 to September 2016. A single LME model was used for 2015 rotenone lakes, contrasting reference and rotenone treatment lakes from May 2015 to September 2016. Additionally, preexisting differences between reference and rotenone lakes before any treatment was applied were tested with LME models for both 2014 rotenone and 2015 rotenone lakes, contrasting metrics from June to September of each respective treatment year.
Fixed effects in LME models included the following: (1) treatment (reference vs. rotenone); (2) period (before vs. after, where after is separated into Year 1 and Year 2 postimpact for 2014 rotenone lakes); and (3) the interaction of period and treatment. A significant interaction was taken as evidence that the treatment had a significant effect and that recovery had not yet occurred (Bro et al., 2004). To control for nonindependence in data over time, month was included as a random effect (Harrison et al., 2018). Lake was also included as a random effect to capture site-to-site variation. To identify the optimal random effects structure, we compared models with no random effects, random slopes, or random slopes and intercepts, selecting the model with the lowest Akaike information criterion (AIC) score (Zuur et al., 2009). Models were fit using restricted maximum likelihood. Model assumptions were examined by visual inspection of QQ plots and distribution of residuals (Harrison et al., 2018). Both total zooplankton and taxonomic group abundance were log 10 (x + 1) transformed to meet assumptions. Given that we were interested in the interaction term (period Â treatment) because of its relevance to understanding recovery, we did not perform model simplification and report on all fixed effects. We also tested for differences in water quality variables (specific conductance, pH, Secchi, dissolved oxygen, and nutrients) using LME models with the same approach described above. Spatial synchrony of zooplankton abundance was examined across four categories of lakes and periods: (1) reference-before; (2) treatment-before; (3) referenceafter; and (4) treatment-after (sample sizes in Appendix S1: Table S2). Although examinations of spatial synchrony are more typically performed on long-term datasets, we focus our question on whether zooplankton recovery was similar across lakes, with the caveat that small sample sizes may limit inference. Spatial synchrony was measured with r i , the intraclass correlation coefficient, as generated from two-way ANOVAs without replication, following Rusak et al. (1999). This method is effective at detecting synchrony over time without the confounding effects of interlake differences. Zooplankton data were detrended prior to analysis to remove potentially confounding effects of linear changes resulting from rotenone treatment. To test whether r i differed from 0, we used a restricted randomization test (n = 9999 permutations) that shuffled zooplankton abundance across months, but preserved data within each lake. These procedures were repeated with each taxonomic group and performed separately for the 2014 and 2015 treatment years. Only the first year following treatment was used for 2014 lakes. We used the false discovery rate test to control for multiple comparisons (Benjamini & Hochberg, 1995) within each study year taxonomic group combination. Given the small number of lakes, especially for rotenone treatment (n 2014 = 4, n 2015 = 3), we did not test for synchronous subsets.
We examined changes in zooplankton community composition with nonmetric multidimensional scaling (NMDS) of all 14 lakes. Patterns of dissimilarity between rotenone-treated and reference lakes before and in the first year after treatment were compared using NMDS with the Bray-Curtis distance metric (Faith et al., 1987). Species abundances were square root transformed to adjust for potential skewed effects from overly abundant species. One time point (McDowell Lake, 11 June 2015) was removed because there were no live zooplankton. We included only those species present in >1 lake and >5% of samples in order to reduce the influence of rare taxa (Cao et al., 2001). Juveniles (nauplii, copepodids, and cladocerans) were excluded from this analysis as they provide little insight on compositional changes, as was Leptodora kindtii, which may not have been accurately sampled with daytime tows. All time points and lakes were run together to compare results across the entire study. Increased variation in community composition of rotenone lakes and reference lakes before and after treatment was tested using a beta dispersion test (M. J. Anderson et al., 2006). PER-MANOVA could not be used because of differences in group dispersion (M. J. Anderson & Walsh, 2013).

Prediction 1
We sought to understand whether between-lake differences in recovery patterns were driven by combinations of functional traits, abiotic, and/or biotic conditions. To test our prediction that cyclopoids would be the first group to recover, we took a multipronged approach, assessing patterns between taxonomic groups with LME models (above), analyzing functional traits, and comparing observed densities to critical densities needed for sexual reproduction.
To assess the role of traits in understanding recovery, we compiled a list of species-level traits from several sources (Barnett et al., 2007, Hébert et al., 2016Appendix S2). As zooplankton species trait information can be sparse for certain taxa and traits, we were able to find comprehensive data on habitat, feeding mode, trophic level, body size, body mass, and egg development rate (details in Appendix S2). From this dataset, we generated a species Â trait matrix, which was populated by 0 or 1 (binary variables), numerical categories (ordinal variable), or numerical values (continuous variables). The species Â trait matrix was multiplied by a site Â species matrix (with relative abundances of species from each lake on each sampling date) to generate a site Â trait matrix, representing abundances of traits, following the approach of Strecker et al. (2011). Traits were standardized to a Z score. With the site Â trait matrix, we performed a redundancy analysis (RDA) to relate multivariate trait abundances to three predictor variables representing recovery: period (0 = before treatment and 1 = after treatment), months post rotenone (discrete), and month of the year (discrete, representing seasonality). Sampling dates prior to rotenone were included and coded as 0 for the variable months post rotenone. We restricted our analyses to treatment lakes only. RDA was appropriate, given the relatively short gradient observed. Predictor variables were not collinear (variance inflation factor < 10). Forward selection of predictor variables was performed following Blanchet et al. (2008).
In order to evaluate the potential importance of Allee effects in limiting recovery, we compared densities of commonly encountered taxa in rotenone lakes to published values of critical densities required for sexual reproduction (Gerritsen, 1980;Gray & Arnott, 2011). Critical densities were compared for Bosmina longirostris, Daphnia pulicaria, diaptomid copepods (Aglaodiaptomus leptopus, Leptodiaptomus novamexicanus, or Skistodiaptomus oregnonensis, depending on which species was present in each lake), and Mesocyclops edax. Copepods are obligately sexual, whereas B. longirostris and D. pulicaria are cyclically asexual. Though some D. pulicaria populations are known to exhibit obligate asexuality, populations in Washington State are thought to primarily exhibit cyclic asexuality (Hebert & Finston, 2001). Temperate bosminids are thought to be cyclically asexual (De Melo & Hebert, 1994). While the timing and outcomes of sexual reproduction differ between cladocerans and copepods (resting eggs and either resting or subitaneous eggs, respectively), in order for populations to reestablish following rotenone, the critical density threshold must be surpassed. We note that there can be significant life history variation within species, but posit that failure to achieve critical densities could be one (of many) competing hypotheses for differential recovery from disturbance. A similar approach to evaluating zooplankton recolonization following acidification was taken by Gray and Arnott (2011). We compared species densities to both high and low published critical densities to encompass a range of possible reproductive outcomes.

Prediction 2
We evaluated the influence of water temperature and Secchi depth (as a proxy of algal biomass) on zooplankton recovery. To do this, we compared zooplankton taxonomic group densities in the summer after rotenone treatment to densities in reference lakes using % change: We chose July and August to compare across rotenone and reference lakes, as these months had the most comprehensive sampling. Mean water column temperature was averaged across July and August, as was Secchi depth. LME models were used to test our prediction, with treatment year as a random variable. Models were compared with the corrected AIC (AIC c ) for small sample sizes, and the importance of water temperature and Secchi depth was tested in separate models because of the small sample size (n = 7). Since all of the models without the random effects had the lowest AIC c , we present the results from these simplified generalized least squares models.

Prediction 3
The effects of fish stocking on recovery patterns of zooplankton diversity across rotenone and reference lakes were assessed using an analysis of covariance (ANCOVA), with the change in Shannon-Wiener (ΔSW) diversity as the response variable: where SW after represents average zooplankton diversity in the 2 months following fish stocking and SW before represents average zooplankton diversity in the 2 months prior to fish stocking. Fish stocking density (number of fish per hectare) was the explanatory variable (Appendix S1: Table S1). Fish stocking typically occurred in the spring following rotenone treatment, between March and June, and both rotenone and reference lakes were stocked. Although the timing of fish stocking could potentially confound changes in diversity with seasonal planktonic succession (i.e., communities are becoming more diverse simply because of seasonal changes in spring), we found that there was no relationship between fish stocking date and ΔSW (linear model: F 1,15 = 0.007, p = 0.932). Reference lakes sampled in both years following treatment were treated as independent replicates, in addition to rotenone lakes (n = 17).
Prediction 4 Lastly, we tested our prediction that there will be asynchrony between taxonomic groups following rotenone. We infer that asynchrony between taxonomic groups is a signal of competition (Vasseur et al., 2005), acknowledging that both synchrony and asynchrony can represent other ecological processes (e.g., predation and environmental fluctuations). We estimated synchrony in treatment lakes with Kendall's W in the R package synchrony v. 0.3.8 (Gouhier & Guichard, 2014). Abundance of each group was detrended prior to analysis (as in spatial synchrony analyses). Kendall's W ranges between 0 and 1, with the former representing lack of synchrony and the latter representing perfect synchrony. Spearman's ranked correlation (ρ) was used to distinguish asynchrony (negative value) from lack of synchrony (zero) (Gouhier & Guichard, 2014). Monte Carlo randomizations (n = 999) shuffled taxonomic groups to determine significance with correct error rates.

Water quality
There were no major changes in physical and chemical parameters following treatment (Figure 2). Secchi depth, surface pH, surface dissolved oxygen, and specific conductance did not show any significant interaction effects (with two exceptions, below) (Table 3, Figure 2). Additionally, there were no significant interaction effects for TP or TN, indicating that nutrient concentrations stayed relatively stable following rotenone (Appendix S3: Tables S1, S2). There was a significant interaction effect for Secchi depth for the first year after treatment in 2015 rotenone lakes and the second year after treatment in 2014 rotenone lakes. Water clarity in the 2015 rotenone lakes decreased by an average of 1.5 m following treatment, while remaining stable in reference lakes. In the second year following treatment for 2014 rotenone lakes, sampling was sporadic, but on average, water clarity increased by 2.2 m following treatment, while reference lakes exhibited a modest mean increase of 0.4 m. Thus, water clarity changes were not consistent across lakes treated in different years. Water clarity increases following rotenone have been observed previously, though this change may be more prevalent in lakes that were initially eutrophic and may be short-lived (Prejs et al., 1997).

Patterns of zooplankton recovery
Prior to the rotenone treatment, zooplankton communities from reference and rotenone lakes were largely similar across metrics of abundance, diversity, and taxonomic group abundance, with the exception of cyclopoid copepods (Tables 4 and 5). Cyclopoids were significantly more abundant in rotenone lakes compared to reference lakes in 2014 but not 2015.
Zooplankton abundance declined immediately following rotenone treatment in 2014 rotenone lakes when compared to relatively stable reference lake abundances in October (99% decline; Figure 3a). Zooplankton abundance in rotenone lakes then showed a rapid increase to pretreatment levels by February of the following year ($4 months after treatment), which exceeded pretreatment levels for 2-3 months, and then remained stable for the duration of the study (Figure 3a). There was a significant interaction between period Â treatment for 2014 rotenone lakes in the first year (Table 4), indicating that rotenone had an overall significant negative effect on zooplankton abundances and that recovery had not yet occurred. This effect was no longer significant in the second year following treatment for the 2014 rotenone lakes ( Table 4), suggesting that recovery had occurred. The same trend was apparent for zooplankton abundance of 2015 rotenone lakes, which showed a decrease in abundances following treatment (99% decline compared to reference lakes in winter) and then a return to pretreatment levels by May of the following year ($7 months after treatment) (Table 4). There was no significant interaction for abundance in the 2015 rotenone lakes, suggestive of recovery.
Shannon-Wiener diversity decreased immediately following treatment by 79% and 71% in 2014 and 2015 rotenone lakes, respectively, with a general increasing trend until June of the following year ($8 months after treatment; Figure 3b). There was a significant interaction between treatment Â period for diversity in the 2014 rotenone lakes in both Year 1 and Year 2, as well as the 2015 rotenone lakes, indicating that diversity still had not fully recovered to pretreatment levels by the end of the study (Table 4).
When zooplankton were split into taxonomic groups, patterns emerged in the rates of response to rotenone  Note: Random effects are reported with marginal R 2 (i.e., fixed effects; R 2 M ) and conditional R 2 (i.e., fixed and random effects; R 2 C ) for each separate model. When "na" is supplied for R 2 M and R 2 C , a model with fixed effects only was selected as the optimal model structure (see Section "Methods" for details). Numerator df is 1 for all effects; denominator df (reported below) varies due to the use of the Satterthwaite approximation. *p < 0.10, **p < 0.05. treatment (see "Prediction 1"). Calanoid copepod (100% decline both years), cladoceran (2014: 99.7%; 2015: 99.6%), and cyclopoid copepod (2014: 89%; 2015: 93%) abundances all declined immediately following rotenone treatment compared to reference lakes (Figure 4). For cyclopoids, which were more abundant in rotenone lakes compared to reference lakes prior to treatment, the decline was even more steep when comparing samples immediately before and after rotenone treatment (2014: 98%; 2015: 99%). The concentration of rotenone used had no effect on % changes in abundance immediately following treatment (generalized least squares, all p > 0.05).
Calanoid copepod abundance increased at the slowest rate (Figure 4a). There was a significant interaction in the 2014 rotenone lakes in both Year 1 and Year 2, as well as 2015 rotenone lakes, suggesting that rotenone had a consistent negative impact, with incomplete recovery even 2 years after treatment (though abundance appeared recovered near the end of the second year for 2014 rotenone lakes; Table 5).
Cladoceran abundance increased at a modest rate following the disturbance, returning to pretreatment abundances by April or May of the following year (6-7 months after treatment) (Figure 4b). There was a significant treatment Â period interaction for cladocerans in the first year for the 2014 rotenone lakes, and a marginally Random effects are reported with marginal R 2 (i.e., fixed effects; R 2 M ) and conditional R 2 (i.e., fixed and random effects; R 2 C ) for each separate model. When "na" is supplied for R 2 M and R 2 C , a model with fixed effects only was selected as the optimal model structure (see "Methods" for details). Numerator df is 1 for all effects; denominator df (reported below) varies due to the use of the Satterthwaite approximation. *p < 0.10, **p < 0.05. significant interaction for 2015 rotenone lakes, but no interaction for Year 2 of the 2014 rotenone lakes, suggesting that recovery had occurred by Year 2 (Table 5).
Cyclopoid abundances increased the fastest following treatment of all rotenone lakes, with returns to pretreatment abundances between February and April of the following year (4-6 months after treatment) (Figure 4c). There was a significant treatment Â period interaction for cyclopoid copepods in the first year of both 2014 and 2015 rotenone lakes, but no interaction for the second year ( Table 5), suggesting that recovery had occurred by Year 2.
Trends in juvenile zooplankton (copepodids, nauplii, and juvenile cladocerans) abundance following rotenone treatment were qualitatively similar to adults (Appendix S3: Figure S1). Calanoid copepodids exhibited steep declines and the slowest recovery, whereas cyclopoid copepodids did not decline to the same extent and rebounded relatively quickly. Notably, nauplii exhibited a rapid recovery following treatments, particularly in the 2014 rotenone lakes, achieving pretreatment densities $January 2015.
Spatial synchrony was not observed in treatment lakes, either before or after rotenone treatment, in both 2014 and 2015 (Table 6). Positive and significant synchrony was observed in cladocerans and total zooplankton in reference lakes after treatment in 2014, as well as calanoid copepods in reference lakes prior to treatment in 2015. We note that relatively small sample sizes and short temporal windows may have limited our ability to detect synchrony across lakes.
Analysis of zooplankton community composition showed increasing variability in treatment lakes following rotenone treatment when compared to relatively stable community composition in reference lakes (Figure 5, Appendix S3: Figure S2). This change in community variability between treatment and reference lakes after rotenone treatment was significant (beta dispersion test: F 3,226 = 32.546; p < 0.001). The stress value of the NMDS was 0.2, which is borderline as values >0.2 should be interpreted with caution; however, this is expected due to ordinating drastically different posttreatment communities, in conjunction with stable reference communities (Clarke, 1993). The increased variability posttreatment can be attributed to changes in the abundance of cyclopoid (e.g., Acanthocyclops robustus and M. edax) and calanoid (e.g., L. novamexicanus and Skistodiaptomus oregonensis) copepods, as well as the cladoceran Chydorus spp. (Figure 5i). Examining each lake individually, most lakes shifted to the bottom left quadrant of the ordination plot immediately following treatment (October) (Figure 5a-g), corresponding with the littoral  Table S1 for sample sizes Chydorus and the cyclopoid A. robustus, but no calanoid copepods. Community composition generally returned to pretreatment community composition by June of the year following rotenone ($8 months after treatment).

Prediction 1
There was support for our prediction that cyclopoid copepods would recover fastest (Figure 4c, Appendix S3: Figure S1). To better understand potential mechanisms, we examined the relationship between zooplankton traits and recovery using RDA ( Figure 6). RDA axes 1 and 2 explained a modest 15.1% of variance in traits, and the overall RDA was significant (F = 9.70, p = 0.001, n = 999 permutations). Significant predictor variables included period and months post rotenone (all p < 0.01). The period after treatment was positively associated with larger body mass, longer egg development times, omnivore-carnivore trophic levels, and raptorial feeding strategies-many of the traits associated with cyclopoid copepods (Appendix S2). The period after treatment was negatively correlated with omnivores and herbivores, filter feeding, and suspension feeding. The number of months post rotenone was positively correlated with pelagic taxa, as well as larger body sizes, and negatively correlated with abundance of littoral taxa. We compared densities of four common taxa in rotenone lakes to critical densities for sexual reproduction Year 2 F I G U R E 4 Monthly averaged log 10 adult abundance of (a) calanoids, (b) cladocerans, and (c) cyclopoids (individuals per cubic meter) for reference, 2014 rotenone, and 2015 rotenone lakes. Error bars represent AE1 SE, red lines indicate date of treatment, and gray arrows represent study periods obtained from the literature. We found that some species (e.g., M. edax and D. pulicaria) rapidly exceeded their critical densities following rotenone and sustained relatively stable populations after reaching their critical density (Figure 7). However, other taxa, such as B. longirostris, never reached their critical densities following treatment and densities remained low or at zero for the rest of the study. Diaptomid calanoid copepods demonstrated variable trends: S. oregonensis (McDowell), A. leptopus (No Name), and L. novamexicanus (Upper Hampton) established relatively stable populations following rotenone in some lakes, though at lower densities compared to pretreatment. Other populations of these diaptomids did not recover (e.g., L. novamexicanus in Widgeon and Lower Hampton lakes and A. leptopus in Rat Lake).

Prediction 2
The influence of abiotic variables on zooplankton recovery was examined in the summer following rotenone. Secchi depth, a proxy of algal biomass, had no effect on zooplankton taxonomic group recovery, measured as % change in abundance between rotenone and reference lakes (calanoids: R 2 = 0.538, p = 0.061; cladocerans: R 2 = 0.059, p = 0.600; Figure 8). Given that cyclopoid copepods were significantly more abundant in rotenone compared to reference lakes before treatment occurred, we instead calculated % change in abundance for this taxonomic group as the difference between the year before and after treatment: We found no relationship between Secchi depth and % change for cyclopoids (R 2 = 0.006, p = 0.865). Water temperature had no effect on recovery of calanoid (R 2 = 0.513, p = 0.070) or cyclopoid (R 2 = 0.210, p = 0.301) copepods, but did have a significant positive effect on cladocerans, with the highest temperatures having the largest % change in abundance (R 2 = 0.615, p = 0.037) (Figure 8).

Prediction 3
In order to test the influence of fish stocking on recovery, we analyzed zooplankton diversity before and after fish stocking that occurred in the spring following rotenone treatment (Figure 9). There was no difference between rotenone and reference lakes in the slope of the relationship between fish stocking density and ΔShannon-Wiener diversity using ANCOVA (F 1,14 = 0.353, p = 0.562). However, we observed a significant overall effect of fish stocking density (F 1,14 = 11.066, p = 0.005). Given the lack of difference between rotenone and reference lakes, we pooled all lakes and found that ΔShannon-Wiener diversity increased significantly with stocking density, indicating that zooplankton communities became more diverse, on average, following fish restocking (R 2 = 0.435, p = 0.004).

Prediction 4
Lastly, we examined synchrony of zooplankton taxonomic groups within each lake as an indicator of competition. Contrary to expectations, we failed to observe asynchrony T A B L E 6 Spatial synchrony of zooplankton across lakes, as estimated by the intraclass correlation coefficient (r i ) with p values from randomizations of zooplankton taxonomic groups in rotenone lakes. We observed significant synchronous patterns in the all but two of the rotenone lakes (Lower Hampton and Widgeon), in addition to nearly all reference lakes in 2014 and 2015 (Table 7), suggestive of seasonal changes.

DISCUSSION
We found that a combination of species traits, abiotic, and biotic factors likely influenced the recovery of zooplankton communities to rotenone. Overall zooplankton abundances recovered relatively quickly to pretreatment levels, but recovery of diversity and community composition took longer. Species traits were likely key determinants in the recovery process, as we found that taxonomic groups recovered at varying rates, including calanoid copepods, which failed to recover within 2 years following rotenone treatment. We failed to observe patterns of spatial synchrony in zooplankton following rotenone treatment, suggesting that patterns of recovery from disturbance may be highly site-specific. Using different metrics to analyze recovery trajectories can lead to various interpretations of a "recovered" community and may be beneficial in describing recovery patterns. These findings give insight into the major drivers of recovery to disturbance and may help guide management of lakes with rotenone treatments.

Species traits as indicators of recovery
Species traits were valuable to understand different recovery processes between taxonomic groups in our study (Table 1). Relative fitness differences between species are a significant mechanism driving patterns of community assembly (HilleRisLambers et al., 2012). Zooplankton species differ in a number of traits including dormancy strategy, reproduction, feeding, and development. We observed that pelagic species and larger taxa (body length) were correlated with time since disturbance, whereas taxa with larger body mass, longer egg development times, raptorial feeding, and omnivorecarnivore trophic levels were associated with the posttreatment period ( Figure 6). Declines in littoral taxa might be expected, as rotenone can persist in lake sediments (discussed below). This potential resiliency of large-bodied taxa (both length and mass) to environmental disturbance contradicts theory that smaller taxa should dominate following disturbance (Odum, 1985). Pulse disturbance theory suggests that traits that are resistant to the disturbance will be selected for, as will traits related to rapid colonization (i.e., the colonization-competition trade-off; Jentsch & White, 2019). In our study, large-bodied D. pulicaria and M. edax recovered quickly from rotenone exposure (Figure 7), suggesting that these taxa were resistant to the disturbance and/or effective colonists. Body size is a "master" trait for zooplankton, representing a host of functions that determine ecological strategies (e.g., feeding, growth, and metabolism; Litchman et al., 2013). Thus, traits related to resistance to disturbance (e.g., diapause at advanced stages) and fast growth (e.g., large size) may be keys to understanding the recovery process. We explore these both below. Dormancy strategies are a key determinant of early community assembly processes following disturbance, as both the intensity and timing of a disturbance can create dynamic emergence responses among taxa following disturbances (Russell et al., 2015). Emergence from diapause in the egg bank and dispersal from surrounding lakes are likely the two major mechanisms of zooplankton recovery (Arnott & Yan, 2002;C aceres, 1998;Dalu et al., 2015). Dispersal from nearby lakes is unimportant in our study since none of our lakes are downstream from another lake and overland dispersal is probably negligible compared to the seed bank (Gray & Arnott, 2011). Cyclopoid copepods recovered the fastest and appeared to reach higher abundances following rotenone treatment than observed in the pretreatment monitoring phase (Figure 4c). Cyclopoids can diapause in a juvenile copepodid stage (Gyllström & Hansson, 2004), potentially allowing them to emerge at an advanced stage of development. However, some cyclopoid species may also diapause as adults (Elgmork, 1980) or even as fertilized adult females (Naess & Nilssen, 1991), giving these taxa a distinct advantage when conditions improve. M. edax is known to diapause as an adult female (Naess & Nilssen, 1991), which may help to explain why this species rebounded so quickly following rotenone (Figure 7). Our analysis of cyclopoid copepodid abundance (Appendix S3) further supports our hypothesis that dormancy strategies may be one mechanism for rapid cyclopoid recovery from rotenone. The exact mechanisms at play remain to be tested experimentally. Additionally, raptorial feeders, like cyclopoids (Barnett et al., 2007), are more mobile and may experience greater predation risk compared to less motile taxa; however, this pressure would be reduced following the removal of planktivorous fish with rotenone. Our results are consistent with similar  Anderson, 1970;Beal & Anderson, 1993). In our study, following a short-term dominance by early recovering species, later successional processes appeared to take hold. Cladocerans showed rapid population increases in spring, likely from a combination of overwintering individuals and emergence from diapause ( Figure 4b). Given that overwintering populations were either absent or very small, we hypothesize that emergence from the egg bank was the main source of the population-this remains to be tested. Cladoceran traits, such as fast growth rates, cyclical asexual reproduction, and efficiency in resource allocation, may be integral to their rapid population growth following disturbance (Haddad et al., 2008). Calanoids were the slowest group to recover in the treated lakes, with abundances not reaching pretreatment levels for the duration of the study (Figure 4a). Calanoids are thought to recover slowly due to emergence from diapause eggs, slow development time, reduced mate encounters at low densities, and selective feeding behaviors (Barnett et al., 2007). The slow recovery of calanoids following disturbance was observed in a similar study (Melaas et al., 2001). Notably, cyclopoid copepods also have relatively slow development rates, but it appears that this barrier to recovery may be overcome if taxa have the ability to diapause at advanced stages of development. Life history traits (e.g., development time and dormancy) may be key indicators of the likelihood and speed of recovery of different taxa.
Understanding Allee effects can also help explain differential recovery across taxa. Using the literature-derived critical density threshold for sexual reproduction (Gerritsen, 1980;Gray & Arnott, 2011), we observed that some species, like M. edax, never dipped below the threshold and thus maintained high densities (Figure 7). By contrast, other taxa, like D. pulicaria, fell below the threshold after rotenone treatment, but quickly exceeded the threshold in spring, likely due to asexual reproduction during favorable periods (Figure 7). Daphnia were able to maintain abundances above the threshold into fall. Daphnia typically use both abiotic (temperature and photoperiod) and biotic (crowding and food availability) cues to induce diapause (Gyllström & Hansson, 2004). Conversely, both B. longirostris and diaptomid copepods were largely unable to surpass the threshold and therefore failed to recover to pretreatment levels. The inability of B. longirostris to exceed the threshold is curious: It is cyclically parthenogenic and would be expected to grow rapidly in the spring and summer, much like Daphnia. However, despite having similar times to first reproduction (B. longirostris: 5.6 days; Daphnia pulex: 6.5 days) and similar egg development times (B. longirostris: 2.4 days at 20 C; D. pulex: 2.7 days at 20 C), B. longirostris produces far fewer eggs per clutch (maximum 15) compared to D. pulex (maximum 82) (Kuns & Sprules, 2000;Lynch, 1980). Additionally, the smaller B. longirostris swims more slowly than Daphnia and has a smaller encounter radius (i.e., the distance at which a conspecific is detected), increasing the critical threshold of B. longirostris compared to D. pulex (Gerritsen, 1980). Although we did not directly quantify Allee effects, it is a plausible mechanism to explain species-specific differences in recovery. Targeted experiments of mating success in the aftermath of rotenone could provide needed support for this hypothesis (Kramer et al., 2008). Thus, body size, life history traits, and reproductive mode may interact to play a critical role in shaping organismal recovery from disturbance.

Biotic controls on recovery
In general, biotic controls of zooplankton community recovery appeared to be weak. There was no signal of competition across taxonomic groups within lakes; indeed, we detected significant positive synchrony in 5 of 7 treatment lakes, as well as nearly all reference lakes (Table 7). This result concurs with that observed by Vasseur et al. (2014), where synchrony was the prevailing pattern for zooplankton, observed over a range of time frames and spatial scales, whereas asynchrony (i.e., compensatory dynamics in Vasseur et al., 2014) was rare. Thus, seasonal patterns, likely driven by both abiotic (e.g., temperature) and biotic (e.g., resources) factors, appear to outweigh the recovery dynamics for crustacean zooplankton in our study. There was some evidence for priority effects in the 2014 rotenone lakes, where, on average, cyclopoid abundance was 3.6Â greater in the year following treatment compared to the pretreatment period, while abundance in reference lakes was 2.5Â greater in the same time frame (Figure 4c). This large increase of early colonizing cyclopoids may have slowed the establishment of other F I G U R E 9 The relationship between fish stocking density (log 10 fish per hectare) and ΔShannon-Wiener zooplankton diversity across treatments, where a positive value indicates diversity increased following fish stocking (n = 17). Regression line across all treatments (AE95% confidence interval): ΔShannon-Wiener diversity = 0.19 Â log 10 stocking density À 0.85 taxonomic groups, though this effect was not observed in 2015 rotenone lakes. Furthermore, although we found a positive effect of fish stocking on zooplankton diversity, potentially demonstrating a release from priority effects, the effect was not different between reference and rotenone lakes (Figure 9). We conclude that the effects of fish stocking on zooplankton diversity appear to be independent of disturbance history, and overall, biotic interactions did not appear to play a significant role in recovery. Predation by other organisms can also influence zooplankton communities. Macroinvertebrates are thought to be moderately affected by rotenone, depending on their body size, habitat, and the presence or absence of gills (Vinson et al., 2010). The planktonic phantom midge larvae, Chaoborus, was present in 11 of our study lakes (5 rotenone and 6 reference). Though our sampling regime was not designed to target this taxon, it was frequently captured in daytime net tows. In two rotenone lakes, Chaoborus abundances declined by 51% and 94% immediately following treatment (abundances were already zero in October for the other three lakes), suggesting that it is sensitive to rotenone. Like the crustaceans, Chaoborus abundance returned quickly the following spring and peaked in the summer, often at much higher average abundances in the summer following treatment compared to the summer before treatment, anywhere from 1.1 to 10.7Â higher abundances (median = 5.3Â). This trend was not observed in reference lakes, where the same year-over-year comparison yielded, on average, a very slight increase in Chaoborus abundance (median = 1.1Â). We hypothesize that Chaoborus may be responding to both relaxed predation and competition from planktivorous fishes following rotenone treatment, in addition to high spring prey abundances ( Figure 3a). In particular, Chaoborus may favor small-bodied taxa, like Bosmina (von Ende & Dempsey, 1981), which may be limiting the recovery of Bosmina in rotenone lakes (Figure 7). We interpret these results cautiously because daytime sampling likely underestimates Chaoborus abundance, but these results lend some evidence to the complex nature of recovery from disturbance.

Abiotic controls on recovery
The direct role of abiotic factors in controlling zooplankton recovery was marginal. This result is a bit surprising, given that the positive synchrony we observed across taxonomic groups within lakes (Table 7) might imply a role for external abiotic forcing, as in Vasseur et al. (2014). However, we also observed that there was no significant spatial synchrony (i.e., synchronous changes across lakes) in treatment lakes (Table 6)-spatial synchrony across lakes within a region is suggestive of regional-scale, extrinsic drivers (Rusak et al., 1999). This latter result reinforces the hypothesized weak role of abiotic variables in our study. Water temperature and Secchi depth (a proxy for algal biomass) failed to explain changes in taxonomic group abundance in rotenone lakes compared to reference lakes in July and August, with one exception. Cladoceran abundance was positively related to average water column temperature, with warmer rotenone lakes having abundances similar to reference lakes (Figure 8). This is not surprising given that cladoceran emergence is more affected by temperature than the emergence of copepods (Jones & Gilbert, 2016) and higher temperatures increase developmental rates (Gillooly et al., 2002). Though not significant, it is worth noting that the two treatment lakes with the highest abundances of calanoid copepods relative to reference lakes also were the warmest and had the smallest Secchi depths, suggestive of more productive conditions (Figure 8a,b). These two lakes, McDowell and No Name, were relatively shallow (mean Z max = 7 m) compared to the other rotenone lakes (mean Z max = 18.4 m) ( Table 2). Additionally, No Name Lake would typically be isothermal by late summer, with anoxia in the bottom 1-2 m (Strecker, unpublished), suggestive of conditions that would be favorable for zooplankton, but less favorable for reintroduced trout, which prefer cooler water. Thus, although the direct role of abiotic features in influencing recovery may be limited, interactions between abiotic conditions (e.g., temperature) and traits (e.g., emergence from diapause and developmental rates) may be critical to understanding recovery from disturbance ( Table 1).

Characteristics of the disturbance
Characteristics of the disturbance itself may also influence recovery (Turner, 2010;White & Pickett, 1985). Though we cannot determine the importance of time since last disturbance or frequency of disturbance as a factor influencing recovery with our study design, we hypothesize that legacies of past disturbances may play an important role in shaping recovery. Cladoceran and total zooplankton abundance in No Name Lake recovered the most compared to all other rotenone lakes by the summer following treatment. No Name Lake was last treated with rotenone in $1949, 66 years prior to our study; all other rotenone lakes were last treated within the previous 13 years (Table 2). Additionally, several lakes have been treated multiple times in the past (WDFW, unpublished). Although it is speculative, it is possible that multiple rotenone applications over time may have depleted the egg bank of zooplankton, which would be especially detrimental to calanoid copepods, as they are thought to have shorter durations of egg viability in the sediment compared to cladocerans . The long history of fish stocking in Sierra Nevada lakes is thought to have depleted the egg bank of the copepod Hesperodiaptomus shoshone and restricted its recovery following cessation of stocking . Thus, historical rotenone treatments may have created legacy effects that persist, altering the ability of communities to recover from disturbance. These legacy effects are a function of both frequency and intensity of disturbance. For instance, the highest intensity disturbance had negative legacy effects on community biomass at different disturbance frequencies in experimental protist communities, suggestive of nonlinear and threshold dynamics potentially preventing recovery from occurring (Jacquet & Altermatt, 2020).
Lastly, the timing of disturbance may be a critical factor in understanding recovery. Rotenone was applied in October of 2014 and 2015 in our study. Other studies of the recovery of zooplankton from rotenone that observed similar recovery patterns were also treated in the fall (Hanson & Butler, 1994;Melaas et al., 2001). Crustacean zooplankton, especially copepods, can be cued to enter a period of dormancy as either diapausing eggs or other life stages, by seasonal changes in temperature and photoperiod, as well as food availability (Gyllström & Hansson, 2004) that can occur in late summer or fall (though timing can vary by species; C aceres, 1998). Rotenone may have a longer half-life in sediment compared to water (Vasquez et al., 2012).
Though there are few studies on the effects of rotenone on dormant life stages, one study suggests that calanoid resting eggs can tolerate rotenone at much higher concentrations compared to later life stages (Naess, 1991). However, resting egg hatching rates in rotenone treatments were reduced compared to pretreatment (Naess, 1991). Additionally, hatching of postdiapause embryos of Artemia franciscana (a stage equivalent to diapause embryos of copepods) was reduced by $80% at ecologically relevant rotenone concentrations (Covi et al., 2016). Critically, this inhibitory effect persisted even after the rotenone was removed, suggesting that this piscicide may have long-term effects on the egg bank of lakes, resulting in the prolonged recovery of certain taxa. Future studies on the effects of rotenone should focus on the dynamics of the egg bank to better understand barriers to species and community recovery.

Implications for disturbance ecology and management
The findings of our study emphasize the need to analyze multiple metrics of community structure to ensure recovery has occurred. Recovery of zooplankton community abundance did not coincide with recovery of community composition or diversity, which occurred much later or was only partially recovered, respectively. This may in part be due to high abundances of early recovering species, whereas species slower to recover still had not reached pretreatment abundances. Thus, the metrics chosen for analysis are vital in evaluating recovery (Frost et al., 2006). Recent studies have demonstrated the need to assess recovery trajectories using a variety of metrics, including multivariate composition data, to ensure that assessments of recovered communities are accurate (Goosem et al., 2016). Utilizing both univariate and multivariate metrics in analysis of recovering communities can be a more robust style of analysis and may help to avoid inadequate evaluations of ecological processes.
Although the results of our study are encouraging for the field of disturbance ecology, there are several limitations to be considered. Due to ice conditions, as well as inaccessibility due to a forest fire, some lakes were sampled more intensely than others following the disturbance. We do not believe that sampling frequency influenced our results because the significant negative effect of rotenone was observed immediately, without corresponding declines in reference lakes (Figure 3). The lack of reference samples in the months immediately following disturbance may have affected our interpretation of community recovery, but a difference in communities would still be likely due to a small but present overwintering zooplankton assemblage in reference lakes (Grosbois et al., 2017). The last limitation to consider is possible effects of rotenone on phytoplankton and rotifers, which were not assessed in our study. Duggan et al. (2015) observed a short-term shift in phytoplankton, and fluctuating rotifer abundances following rotenone treatment, but no long-term effects, indicating that food resource changes would likely not significantly alter the zooplankton community in the following years.
The recovery of communities to whole-ecosystem-scale disturbances is a dynamic process dependent on both biotic and abiotic characteristics. Species traits are a critical aspect of the recovery process, as both dormancy strategies and development time are important characteristics for resource managers to consider for recovering communities (Fraterrigo & Rusak, 2008). Priority effects can shape the recovery process but may depend on abiotic conditions and additional management actions following disturbance, such as fish restocking (Duggan et al., 2015;Tucker & Fukami, 2014). Further, though we did not quantify the long-term legacy effects of historical disturbance, these may have profound effects on community resilience and recovery, ultimately altering the processes that maintain ecosystem function.