Effects of Grazing and Nutrients on Phytoplankton Blooms and Microplankton Assemblage Structure in Four Temperate Lakes Spanning a Eutrophication Gradient

Phytoplankton assemblage dynamics are sensitive to biotic and abiotic factors, as well as anthropogenic stressors such as eutrophication, and thus are likely to vary between lakes of differing trophic state. We selected four lakes in Washington State, USA, ranging from oligo- to hypereutrophic, to study the separate and interactive effects of enhanced nutrient availability and zooplankton grazing on phytoplankton net growth rates and overall microplankton (phytoplankton and microzooplankton) assemblage structure. We collected water quality and plankton samples monthly in each lake from May to October 2014, and also conducted laboratory incubation experiments using ambient plankton assemblages from each lake with amendments of zooplankton grazers (5× ambient densities) and nutrients (Nitrogen + Phosphorus) in June, August, and October. In each set of monthly experiments, nested two-way ANOVAs were used to test the effects of enhanced grazers and nutrients on net chlorophyll a-based phytoplankton growth rates. Nested PERMANOVAs were used to test the effects of each factor on microplankton assemblage structure. Enhanced grazing reduced phytoplankton net growth in oligotrophic Cle Elum Lake and oligo-mesotrophic Lake Merwin in August (p < 0.001) and Merwin again in October (p < 0.05), while nutrient enhancement increased phytoplankton net growth in Lake Merwin in June (p < 0.01). Changes in microplankton assemblage composition were not detected as a result of either factor, but they were significantly different between sites (p < 0.001) during each month, and varied by month within each lake. Significant effects of both enhanced grazers and nutrients were detected in systems of low, but not high, trophic state, although this varied by season. We suggest that it is critical to consider trophic state when predicting the response of phytoplankton to bottom-up and top-down factors in lakes.


Introduction
Phytoplankton are a critical source of autochthonous energy for aquatic food webs, and phytoplankton assemblages naturally vary in biomass and species composition across a wide range of aquatic systems. However, eutrophication, an especially problematic anthropogenic stressor acting on freshwater lakes and rivers [1], is impacting the dynamics of phytoplankton assemblages in multiple and potentially ecologically damaging ways.
For example, increased availability of dissolved inorganic nitrogen (N) and phosphorus (P) may increase phytoplankton biomass from the "bottom up" by providing additional resources to support phytoplankton growth, and may also alter assemblage structure [2][3][4].
Indeed, multiple studies have demonstrated a positive relationship between phytoplankton biomass and nutrient availability in a range of natural lakes [5][6][7]. However, phytoplankton responses to enhanced nutrients often vary by species based on their specific functional traits [8,9] and may be further dependent upon co-occurring environmental conditions (e.g., temperature, light, and mixing) which alter photosynthetic and metabolic rates [10].
These differences in phytoplankton across systems, including species-specific traits interacting with abiotic factors, make predicting the response of phytoplankton biomass and assemblage structure to eutrophication somewhat challenging [11].
Eutrophication of lakes is also of particular concern since higher frequency and intensity of "harmful" phytoplankton blooms have been associated with increased nutrient concentrations occurring under conditions associated with climate change, e.g., higher temperature and increased periods of drought [12,13]. These conditions have been shown to especially stimulate the proliferation of cyanobacteria over other phytoplankton taxa [14][15][16], which is problematic due to their potential for rapid growth and the production of toxins that endanger humans and wildlife, as well as reduce water quality, dissolved oxygen concentrations, aquatic plant and epiphyte abundance, and biodiversity [3,17,18].
In addition to resource availability (i.e., "bottom-up" drivers), zooplankton grazing may exert control of phytoplankton assemblages from the "top down" [19], with differential feeding strategies and diet preferences shaping phytoplankton total biomass, bloom dynamics, and taxonomic composition [20,21]. For instance, daphnid cladocerans (who along with copepods are the dominant mesozooplankton (200-2000 µm) grazers in temperate lakes) are generally observed to limit phytoplankton biomass in lakes of varying size and depth due to their capacity to efficiently and non-selectively graze phytoplankton cells across a wide size range [22,23]. Copepods, on the other hand, primarily influence phytoplankton biomass and taxonomic composition through preferential consumption of microzooplankton (~20-200 µm, in particular ciliates) which releases certain phytoplankton taxa from grazing pressure through a planktonic trophic cascade [21,[24][25][26]. Of note, several species of Daphnia are known to consume and limit growth of specific cyanobacterial taxa [27,28], while a number of copepod taxa have been observed to contribute to cyanobacterial dominance by selectively ingesting eukaryotic algae or microzooplankton grazers of cyanobacteria [25,29,30].
Previous research has shown that both bottom-up [31,32] and top-down [21,24,33] factors are known to affect phytoplankton growth and assemblage structure, although the relative influence of these factors has been observed to vary across lakes of different size, depth, and latitude [5,[34][35][36]. The "trophic state" (i.e., combination of nutrient availability, algal biomass, and water clarity; [37][38][39]) of a lake may also influence this balance. In oligotrophic lakes, phytoplankton biomass is typically quite low, with the assemblage dominated by small cells [40], conditions under which Daphnia and other large cladocerans may have a relatively high grazing impact (e.g., [41][42][43]). In mesotrophic lakes, nutrient availability is often sufficient to support phytoplankton growth; however, mesozooplankton assemblages typically consist of both daphnid cladocerans and copepods, who together may exert grazing pressure that controls phytoplankton biomass and diversity [35,44,45]. Finally, in eutrophic and hypereutrophic lakes, high nutrient availability supports high phytoplankton growth rates and biomass (e.g., [6,46]), while mesozooplankton grazing may have a limited direct impact on phytoplankton biomass due to the dominance of small and/or inedible phytoplankton taxa [25,26,42].
Given these observed differences in phytoplankton response to nutrients and grazing in lakes of varying trophic state, along with the risk of more intense and/or harmful phytoplankton blooms occurring due to eutrophication, it is clear that a better understanding of nutrient-phytoplankton-grazer interactions is warranted in order to best model and manage these critical freshwater resources.
In response to this need, we experimentally tested the separate and interacting effects of enhanced grazing of metazoan zooplankton and enhanced nutrient availability on phytoplankton growth rates and the structure of the overall microplankton assemblage (here defined as phytoplankton plus unicellular heterotrophic protists) over the 2014 growing season (May-October) in four western Washington (USA) lakes that varied across a trophic gradient. We complemented our experiments with monthly field observations of dissolved inorganic nutrients (N and P), water temperature, Secchi depth, and microplankton and metazoan zooplankton abundance and composition. We expected that the effects of enhanced grazers and nutrients on phytoplankton growth and assemblage structure would vary in lakes across a gradient of trophic state, as illustrated in our simple conceptual model ( Figure 1). We categorized our four lake systems from oligotrophic to hypereutrophic according to patterns in the concentrations of P (soluble reactive phosphorus) and chlorophyll (chl) a, and lake depth. Based on the literature reviewed above, we hypothesized that the addition of replete nutrients (N and P) would have the greatest stimulating effect on the phytoplankton assemblage in lakes of low trophic state where nutrients are expected to be limiting, and that this effect would decrease in a linear fashion in lakes of higher trophic state. We further hypothesized that in oligotrophic lakes, where ambient mesozooplankton abundance was expected to be low to moderate, the effect of enhanced mesozooplankton grazing would be relatively high; but that additional mesozooplankton grazing pressure would have little to no effect on phytoplankton growth or assemblage structure in mesotrophic (since zooplankton grazing impact in these systems is already high) or eutrophic lakes (where phytoplankton are less efficiently grazed due to size and/or inedibility).
Water 2021, 13, x FOR PEER REVIEW 3 of 20 season (May-October) in four western Washington (USA) lakes that varied across a trophic gradient. We complemented our experiments with monthly field observations of dissolved inorganic nutrients (N and P), water temperature, Secchi depth, and microplankton and metazoan zooplankton abundance and composition. We expected that the effects of enhanced grazers and nutrients on phytoplankton growth and assemblage structure would vary in lakes across a gradient of trophic state, as illustrated in our simple conceptual model ( Figure 1). We categorized our four lake systems from oligotrophic to hypereutrophic according to patterns in the concentrations of P (soluble reactive phosphorus) and chlorophyll (chl) a, and lake depth. Based on the literature reviewed above, we hypothesized that the addition of replete nutrients (N and P) would have the greatest stimulating effect on the phytoplankton assemblage in lakes of low trophic state where nutrients are expected to be limiting, and that this effect would decrease in a linear fashion in lakes of higher trophic state. We further hypothesized that in oligotrophic lakes, where ambient mesozooplankton abundance was expected to be low to moderate, the effect of enhanced mesozooplankton grazing would be relatively high; but that additional mesozooplankton grazing pressure would have little to no effect on phytoplankton growth or assemblage structure in mesotrophic (since zooplankton grazing impact in these systems is already high) or eutrophic lakes (where phytoplankton are less efficiently grazed due to size and/or inedibility). level of effect on phytoplankton growth and assemblage structure predicted from the addition of nutrients and mesozooplankton grazers in four lakes of increasing "trophic state," as indicated by mean concentrations of dissolved PO4 and chlorophyll (chl) a, and total lake depth.

Study Sites
Cle Elum Lake, Lake Merwin, Lacamas Lake, and Vancouver Lake, all in Washington (WA) state, USA (Figure 2), were selected in this study because they span a gradient of eutrophy, including oligotrophic (Cle Elum), oligo-mesotrophic (Merwin), eutrophic (Lacamas), and hypereutrophic (Vancouver) lakes. level of effect on phytoplankton growth and assemblage structure predicted from the addition of nutrients and mesozooplankton grazers in four lakes of increasing "trophic state," as indicated by mean concentrations of dissolved PO 4 and chlorophyll (chl) a, and total lake depth.

Study Sites
Cle Elum Lake, Lake Merwin, Lacamas Lake, and Vancouver Lake, all in Washington (WA) state, USA (Figure 2), were selected in this study because they span a gradient of eutrophy, including oligotrophic (Cle Elum), oligo-mesotrophic (Merwin), eutrophic (Lacamas), and hypereutrophic (Vancouver) lakes. Cle Elum Lake (47.2796° N, 121.1041° W) is a large (19 km 2 ), deep (~90 m), dimictic reservoir of the Yakima River, a major tributary of the Columbia River (CR) in central WA, that was impounded in 1933 to meet the water needs of the Yakima Basin. Reservoir depth fluctuates annually as water is drawn down from the hypolimnion, primarily for irrigation [47]. The reservoir is located at high elevation (~678 m) in the Cascade Mountain Range, surrounded by largely forested land and wilderness.
Lake Merwin (45.9824° N, 122.4687° W) is a large (15.9 km 2 ), deep (mean = 31 m) reservoir of the North Fork of the Lewis River, a tributary of the CR in southwest WA [48]. Surface elevation of Lake Merwin is approximately 74 m. Merwin dam, built in 1931, is part of the Lewis River Hydroelectric project, operated primarily for power generation and fish passage. The surrounding watershed consists largely of land within the Gifford Pinchot National Forest and the site is a popular location for recreation and fishing.
Lacamas Lake (45.6184° N, 122.4276° W) is a small (1.3 km 2 ) eutrophic reservoir (mean depth = 7.8 m) located in southwest WA that was impounded in 1938 and has a surface elevation of approximately 57 m. The lake is monomictic and stratifies in summer months, during which the hypolimnion may become hypoxic [49]. Water level is substantially drawn down each fall, as the result of an epilimnetic spill event through near-surface gates at the dam [50]. Water from the lake flows from Lacamas Creek to the Washougal River, a tributary of the CR. Lacamas Lake has high summertime recreational use and receives regular public health advisories due to toxic cyanobacterial blooms in late summer.
Vancouver Lake (45.6747° N, 122.7279° W) is a broad (~9.3 km 2 ), shallow (mean depth ~1 m) urban floodplain lake located within the city of Vancouver, in southwest WA. The lake has a surface elevation of approximately 8 m and is connected to the CR by a gated channel on the southwest shore. The lake is a popular recreational site within the Vancouver, WA and Portland, OR metropolitan area, and has experienced recurrent harmful cyanobacterial blooms in late summer for several decades [26,46,51].

Field Sampling
Sampling was conducted monthly in each lake from May-October 2014. Surface water was collected in triplicate with a clean bucket, from which a 200-mL subsample was placed into amber bottles containing 5% Lugol's acid solution to preserve microplankton. For analysis of dissolved nutrient concentrations (PO4-P, NH4-N, NO3-N, NO2-N), Cle Elum Lake (47.2796 • N, 121.1041 • W) is a large (19 km 2 ), deep (~90 m), dimictic reservoir of the Yakima River, a major tributary of the Columbia River (CR) in central WA, that was impounded in 1933 to meet the water needs of the Yakima Basin. Reservoir depth fluctuates annually as water is drawn down from the hypolimnion, primarily for irrigation [47]. The reservoir is located at high elevation (~678 m) in the Cascade Mountain Range, surrounded by largely forested land and wilderness.
Lake Merwin (45.9824 • N, 122.4687 • W) is a large (15.9 km 2 ), deep (mean = 31 m) reservoir of the North Fork of the Lewis River, a tributary of the CR in southwest WA [48]. Surface elevation of Lake Merwin is approximately 74 m. Merwin dam, built in 1931, is part of the Lewis River Hydroelectric project, operated primarily for power generation and fish passage. The surrounding watershed consists largely of land within the Gifford Pinchot National Forest and the site is a popular location for recreation and fishing.
Lacamas Lake (45.6184 • N, 122.4276 • W) is a small (1.3 km 2 ) eutrophic reservoir (mean depth = 7.8 m) located in southwest WA that was impounded in 1938 and has a surface elevation of approximately 57 m. The lake is monomictic and stratifies in summer months, during which the hypolimnion may become hypoxic [49]. Water level is substantially drawn down each fall, as the result of an epilimnetic spill event through near-surface gates at the dam [50]. Water from the lake flows from Lacamas Creek to the Washougal River, a tributary of the CR. Lacamas Lake has high summertime recreational use and receives regular public health advisories due to toxic cyanobacterial blooms in late summer.
Vancouver Lake (45.6747 • N, 122.7279 • W) is a broad (~9.3 km 2 ), shallow (mean depth 1 m) urban floodplain lake located within the city of Vancouver, in southwest WA. The lake has a surface elevation of approximately 8 m and is connected to the CR by a gated channel on the southwest shore. The lake is a popular recreational site within the Vancouver, WA and Portland, OR metropolitan area, and has experienced recurrent harmful cyanobacterial blooms in late summer for several decades [26,46,51].

Field Sampling
Sampling was conducted monthly in each lake from May-October 2014. Surface water was collected in triplicate with a clean bucket, from which a 200-mL subsample was placed into amber bottles containing 5% Lugol's acid solution to preserve microplankton. For analysis of dissolved nutrient concentrations (PO 4 -P, NH 4 -N, NO 3 -N, NO 2 -N), subsamples from each bucket were filtered through 0.45-µm filters, frozen, and shipped to the University of Washington's Marine Chemistry Laboratory for measurement following the protocols of the WOCE Hydrographic Program using a Seal Analytical AA3 auto-analyzer. For measurement of chl a concentration, a 50-mL aliquot from each bucket was stored in an amber bottle in the field and filtered through Whatman GF/F filters upon return to the laboratory at Washington State University Vancouver. Filters were frozen for at least 24 h (but no longer than 7 days) before undergoing 24 h of extraction in 90% acetone. Chl a concentration was quantified using a Turner Model 10 AU fluorometer following the acidification method [52].
Metazoan zooplankton were collected on each sampling date via triplicate vertical tows from 0.5 m off the lake bottom to the surface, using a 73-µm mesh, 0.5-m diameter plankton net and preserved with 5% buffered formalin in the field. Water temperature (measured with a YSI95 probe [YSI/Xylem Inc., Yellow Springs, OH, USA]), lake depth, and water clarity (measured as Secchi disk depth) were also measured on each sampling date.

Grazing and Nutrient Enhancement Experiments
Two-factorial grazing and nutrient enhancement experiments were conducted three times using water and plankton from each lake (n = 12) over a summer bloom cycle, allowing for assessment of chl a-based phytoplankton growth rates in response to treatments during the pre-(June), mid-(August), and post-bloom (October) periods. For each experiment, twenty 500-mL polycarbonate incubation bottles were filled with unfiltered lake water, obtained using an acid-washed carboy. Four replicate initial control bottles containing only lake water with the natural plankton assemblage were immediately subsampled for measurement of both chl a concentration and microplankton composition and abundance (as described above). One set of quadruplicate final control bottles and three sets of quadruplicate treatment bottles were prepared as follows: "Control" bottles contained only the natural assemblage of plankton; "Grazer" bottles contained the natural assemblage of plankton plus an amendment of the numerically-dominant zooplankton grazer taxon (described in detail below); "Nutrient" bottles contained the natural assemblage of plankton with an amendment of N and P (described in detail below); and "Both" bottles contained the natural assemblage of plankton plus amendments of both grazers and nutrients. Control and treatment bottles were sealed with Parafilm and a lid to prevent air space and turbulence-causing bubbles.
Nutrient amendments consisted of the addition of NaNO 3 (454 µg/L) and NaH 2 PO 4 ·H 2 O (108 µg/L) to treatment bottles, based on the highest mean concentrations of both constituents occurring in Vancouver Lake, the most eutrophic system included in our study [46,53]. All bottles were incubated for 24 h under white light, at the ambient temperature and light:dark cycle, on a slowly rotating (0.5-1 rpm) plankton wheel in a temperature-controlled room. At the end of the incubation period all final bottles were subsampled to measure chl a concentration and the microplankton assemblage as described above. We chose a 24-h incubation period to provide sufficient time for a phytoplankton growth response to occur, as has been observed previously in the most eutrophic lake [51], but not so long as to introduce potential error and confounding effects associated with increased zooplankton respiration and excretion, and interactions among grazers.
Approximately 5-8 h prior to each experiment, live zooplankton samples were collected from a given lake using a 73-µm mesh, 0.5-m diameter plankton net with a live cod end. Upon return to the laboratory, a subsample of the net tow was preserved and then examined microscopically to determine the numerically-dominant taxon in the 200-2000 µm size range and to estimate ambient zooplankton density. We utilized > 200 µm-sized grazer taxa for amendments due to their ability to ingest a broad array of microplankton taxa [33]. Live individuals of the numerically-dominant taxon were then selectively transferred from the net tow sample into holding beakers containing filtered lake water and held from 4 to 6 h to allow for acclimation and gut emptying [54]. Adult zooplankton grazers were amended at 5× ambient lake density in each "Grazer" incubation bottle.

Taxonomic Enumeration and Identification
For identification of microplankton, Lugol's preserved samples were settled for 12-48 h in Utermöhl chambers before examination with a Leica DMI 4000B (Leica Microsystems, Inc., Buffalo Grove, IL, USA) inverted microscope. At least 300 organisms were counted per sample from fields along transects, as described in [55]. Microplankton cells were enumerated, sized, and identified to genus or species, where possible, using [56] and [57]. Biomass (µg C/L) was estimated using formulas derived in [58] and [59]. To ensure that larger, and potentially rare, taxa were not excluded from counts, organisms > 30 µm in size were counted at 200× and organisms ≤ 30 µm in size were counted at 400×, with at least 50 organisms for each sample counted at 200×. In order to keep units consistent for comparison among all microplankton taxa (especially eukaryotic algae and cyanobacteria), the counting unit was individual cells, rather than colonies.

Rate Calculations and Statistical Analyses
Net phytoplankton growth rates (r) for all experimental incubations were calculated based on changes in chl a concentration using the formula r = ln(C f /C i )/t, where C f is concentration at the end of the incubation, C i is concentration at the start of the incubation, and t is the incubation period (1 d). Nested two-way ANOVA was used to analyze the effect of treatments on growth rates, where net phytoplankton growth was the response variable and the treatment (fixed) factors (i.e., grazers and nutrients and their interaction) were nested within lake (random factor). Three ANOVA tests were conducted (one each for June, August, and October) to reduce disparity of variance and to separate spatially and temporally autocorrelated data. In August and October, data were cube root-transformed to meet test assumptions of normality and homogeneity of variance. ANOVA analysis was conducted using linear mixed effects models, using the R package 'lmerTest' [63], because of its soundness in analyzing unbalanced data, crossed factors, nested designs, and inclusion of random effects, along with type III sums of squares, recommended for unbalanced designs and when the importance of factor order is not known. Where significance of a main effect was found by ANOVA, post-hoc comparisons were made by means of t-tests within an experiment for each lake.
Nested factorial PERMANOVA [64] was used to test the effects of grazers, nutrients, and their interaction on the full microplankton assemblage from each lake in experiments conducted in June, August, and October. As was done for analysis of growth rates, three tests were conducted (one for each month). PERMANOVA was selected because of its robustness in handling crossed and nested designs as well as non-normal data typical of ecological communities [64]. The analysis of microplankton assemblage structure was based on taxonomic counts (log 10 +1 transformed abundance) and PERMANOVA was conducted using 1000 permutations and a dissimilarity matrix calculated with Bray-Curtis distance. Homogeneity of variance was tested using an examination of beta dispersion in each PERMANOVA model. Nonmetric Multidimensional Scaling (NMDS), also conducted using log 10 +1 transformed abundance and Bray-Curtis distance, was used to visually explore patterns in microplankton taxonomic composition between all experimental incubations (3 control and 9 treatment replicates per experiment) and separately for each experimental month as a visual accompaniment to PERMANOVA. Multivariate analysis was conducted using the package 'vegan' [65]. All data were analyzed using R software, version 3.6.1 (available at https://www.r-project.org/ (accessed on 8 April 2021)) [66].

Field Observations
A gradient of biological, chemical, and physical factors was observed across our four study systems from May to October 2014 ( Figure 3). Monthly nutrient concentrations varied over the study period in each system, but were generally lower in Merwin, higher in Cle Elum, and highest in Lacamas and Vancouver lakes (Figure 3). While nutrient concentrations were comparable between systems, two notable differences were apparent: nitrate plus nitrite (NO 3 +NO 2 -N) was dramatically higher in Lacamas from May to August (580.3-105.1 µg /L), and phosphate (PO 4 -P) was markedly highest in Vancouver from August to October (69.4-11.9 µg/L). Similarly, chl a concentration was lowest in Cle Elum (0.4-1.3 µg/L) and Merwin lakes (0.9-5.9 µg/L) throughout the study period, higher and similar in Vancouver was conducted using the package 'vegan' [65]. All data were analyzed using R software, version 3.6.1 (available at https://www.r-project.org/ (accessed on 8 April 2021)) [66].

Field Observations
A gradient of biological, chemical, and physical factors was observed across our four study systems from May to October 2014 ( Figure 3). Monthly nutrient concentrations varied over the study period in each system, but were generally lower in Merwin, higher in Cle Elum, and highest in Lacamas and Vancouver lakes (Figure 3). While nutrient concentrations were comparable between systems, two notable differences were apparent: nitrate plus nitrite (NO3+NO2-N) was dramatically higher in Lacamas from May to August (580.3-105.1 µg /L), and phosphate (PO4-P) was markedly highest in Vancouver from August to October (69.4-11.9 µg/L). Similarly, chl a concentration was lowest in Cle Elum (0.4-1.3 µg/L) and Merwin lakes (0.9-5.9 µg/L) throughout the study period, higher and similar in Vancouver

Succession of Microplankton and Zooplankton Assemblage Structure
Similar to the seasonal pattern of chl a concentration in all four lakes, microplankton biomass was lowest in Cle Elum and increased from Merwin to Lacamas and Vancouver lakes, with variation in the timing of peak biomass between systems (Figure 4). Diatoms commonly constituted the majority of microplankton biomass in all four lakes, with Cle Elum and Merwin also containing high proportions of flagellate and ciliate biomass, and Lacamas and Vancouver lakes containing a high proportion of cyanobacterial biomass in late summer. Cyanobacteria were more common in Merwin than Cle Elum in August through October, but they were not a major contributor of microplankton biomass in either lake. lakes, with variation in the timing of peak biomass between systems (Figure 4). Diatoms commonly constituted the majority of microplankton biomass in all four lakes, with Cle Elum and Merwin also containing high proportions of flagellate and ciliate biomass, and Lacamas and Vancouver lakes containing a high proportion of cyanobacterial biomass in late summer. Cyanobacteria were more common in Merwin than Cle Elum in August through October, but they were not a major contributor of microplankton biomass in either lake.
In Cle Elum, peak microplankton biomass occurred in May and June, declined in July, and experienced a second smaller peak in September and October. Diatoms, flagellates, and dinoflagellates typically comprised the greatest autotrophic portion of microplankton biomass, while ciliates constituted a major portion of assemblage biomass from May to July and again in October (Figure 4). In Merwin, peak microplankton biomass occurred in May and steadily declined through the summer and into the fall. Diatoms and flagellates dominated the community biomass in Merwin from May through June, after which the proportion of diatoms decreased, and the proportion of ciliates and dinoflagellates increased (Figure 4). In Lacamas microplankton biomass was greatest in May and declined to a low in August before experiencing a second smaller peak in September and October. Diatoms comprised the greatest portion of the assemblage biomass in May and June but then decreased, while cyanobacteria In Cle Elum, peak microplankton biomass occurred in May and June, declined in July, and experienced a second smaller peak in September and October. Diatoms, flagellates, and dinoflagellates typically comprised the greatest autotrophic portion of microplankton biomass, while ciliates constituted a major portion of assemblage biomass from May to July and again in October (Figure 4).
In Merwin, peak microplankton biomass occurred in May and steadily declined through the summer and into the fall. Diatoms and flagellates dominated the community biomass in Merwin from May through June, after which the proportion of diatoms decreased, and the proportion of ciliates and dinoflagellates increased (Figure 4). In Lacamas microplankton biomass was greatest in May and declined to a low in August before experiencing a second smaller peak in September and October. Diatoms comprised the greatest portion of the assemblage biomass in May and June but then decreased, while cyanobacteria concurrently increased in biomass from July to October. Cyanobacteria dominated the microplankton assemblage biomass in Lacamas from August through October (Figure 4).
In Vancouver, microplankton biomass increased through the summer and was greatest in August and September. The microplankton assemblage was more mixed in Vancouver than the other lakes; diatoms made up a large portion of the assemblage biomass through-out the study period, with a notable proportion of green algae and dinoflagellates present in May and June, and cyanobacteria were dominant in July and September (Figure 4).
Zooplankton abundance was lowest in Cle Elum, moderately greater in Merwin and Lacamas, and much higher in Vancouver from May through October 2014 ( Figure 5). Cle Elum typically contained high proportions of copepods and daphnid cladocerans throughout the study period, with copepods dominant and rotifers subdominant in May and an increasing proportion of daphnid cladocerans present throughout the summer and into fall. In Merwin, daphnids were frequently dominant from May through October, with a high proportion of rotifers present in May and bosminid cladocerans in August through October. Zooplankton abundance was similar in Vancouver and Lacamas in May and June, but abundance steadily declined in Lacamas and increased in Vancouver throughout the summer. In Lacamas, daphnid cladocerans were dominant in May and June but varied in dominance with copepods from July to October. Vancouver had a more mixed zooplankton assemblage compared with the other lakes. While copepods were more consistently present than daphnid cladocerans, high relative abundances of smaller taxa, including rotifers (in June and August) and bosminids (in July and October), were more common in Vancouver than in any of the other systems.
concurrently increased in biomass from July to October. Cyanobacteria dominated the microplankton assemblage biomass in Lacamas from August through October (Figure 4).
In Vancouver, microplankton biomass increased through the summer and was greatest in August and September. The microplankton assemblage was more mixed in Vancouver than the other lakes; diatoms made up a large portion of the assemblage biomass throughout the study period, with a notable proportion of green algae and dinoflagellates present in May and June, and cyanobacteria were dominant in July and September ( Figure  4).
Zooplankton abundance was lowest in Cle Elum, moderately greater in Merwin and Lacamas, and much higher in Vancouver from May through October 2014 ( Figure 5). Cle Elum typically contained high proportions of copepods and daphnid cladocerans throughout the study period, with copepods dominant and rotifers subdominant in May and an increasing proportion of daphnid cladocerans present throughout the summer and into fall. In Merwin, daphnids were frequently dominant from May through October, with a high proportion of rotifers present in May and bosminid cladocerans in August through October. Zooplankton abundance was similar in Vancouver and Lacamas in May and June, but abundance steadily declined in Lacamas and increased in Vancouver throughout the summer. In Lacamas, daphnid cladocerans were dominant in May and June but varied in dominance with copepods from July to October. Vancouver had a more mixed zooplankton assemblage compared with the other lakes. While copepods were more consistently present than daphnid cladocerans, high relative abundances of smaller taxa, including rotifers (in June and August) and bosminids (in July and October), were more common in Vancouver than in any of the other systems.

Responses to Grazer and Nutrient Enhancement
We conducted 12 experiments to assess the effects of amended nutrients and amended zooplankton grazers on algal growth rates in each lake three times between June and October 2014. In each experiment, incubations were conducted under ambient temperature and light cycles, and the grazer taxon amended in each experiment aligned with the taxon that was most abundant in that lake at that time (Table 1). Table 1. Dates, locations, temperatures, dark:light cycles, and mesozooplankton grazer taxa amended in the incubation experiments conducted from June to October 2014. * = Bosmina spp. were most abundant in Vancouver Lake during the October experiment, but could not be amended at 5× ambient levels; Acanthocyclops robustus was second most abundant and could be amended at 5× ambient levels. Phytoplankton chl a-based net growth rates (d −1 ) varied widely by month and lake throughout our experiments ( Figure 6). In June (expected pre-bloom conditions), the results of nested two-way ANOVA measured for Cle Elum, Lacamas, and Vancouver lakes revealed no significant effect of enhanced grazers, enhanced nutrients, or the interaction of grazers and nutrients on phytoplankton net growth rates ( Table 2). A one-way ANOVA on phytoplankton growth rates from Merwin in June, conducted separately due to the lack of data for the interaction treatment, showed that nutrients significantly increased phytoplankton growth (F = 12.4, p < 0.01) compared to the grazer-enhanced and control treatments. In August (expected bloom conditions), a nested two-way ANOVA revealed that enhanced grazers had a significant effect on phytoplankton growth (F = 16.04, p < 0.001), and post-hoc comparisons found that phytoplankton growth rates were reduced by enhanced grazers in Merwin and Cle Elum (both at p < 0.001). In October (expected post-bloom conditions), enhanced grazers significantly reduced phytoplankton growth (F = 5.67, p < 0.05) in Merwin (p < 0.05). Overall, phytoplankton net growth rates were significantly different between lakes in June, August, and October, as the random nesting factor (lake) was found to be highly significant (all at p < 0.001) in the two-way ANOVA conducted on experimental data from each month ( Table 2).

Location
The NMDS results revealed well-defined clusters associated with each lake and month, indicating strong differences in species composition between lakes and between early summer (June) and late summer/fall (August and October) assemblages (Figure 7). Nested PERMANOVA analyses, conducted separately for each experimental month, found no significant effect of either enhanced nutrients, enhanced grazers, or their interaction on the microplankton composition of each lake in each month (Figure 8, Table 3). However, PERMANOVA revealed that microplankton composition was highly significantly different between lakes in all three months (all at p < 0.001). 021, 13, x FOR PEER REVIEW 11 of 20   Table 2. Results of nested two-way ANOVAs testing the effect of enhanced grazers and enhanced nutrients on chl a-based phytoplankton net growth rates (per day). ANOVAs were conducted separately for each month and utilized linear mixed effects methods to include nesting within lake (site) as a random factor. The ANOVA in June did not include Lake Merwin, but all lakes were included in analyses in August and October. * = p < 0.05, *** = p < 0.001. <0.001 *** month, indicating strong differences in species composition between lakes and between early summer (June) and late summer/fall (August and October) assemblages (Figure 7). Nested PERMANOVA analyses, conducted separately for each experimental month, found no significant effect of either enhanced nutrients, enhanced grazers, or their interaction on the microplankton composition of each lake in each month ( Figure 8, Table 3). However, PERMANOVA revealed that microplankton composition was highly significantly different between lakes in all three months (all at p < 0.001).

Discussion
Overall, our results show that the enhancement of zooplankton grazers and/or nutrients (N + P) affected phytoplankton growth rates only in lakes of a lower trophic state  Table 3. Results of nested two-way PERMANOVAs conducted separately in June, August, and October, testing the effects of enhanced grazers and enhanced nutrients, nested within lake, on microplankton assemblage structure.

Discussion
Overall, our results show that the enhancement of zooplankton grazers and/or nutrients (N + P) affected phytoplankton growth rates only in lakes of a lower trophic state (Merwin and Cle Elum). Enhanced zooplankton grazers significantly reduced phytoplankton net growth in Cle Elum and Merwin lakes in August, and Merwin again in October, and enhanced nutrients only significantly increased phytoplankton net growth in Merwin in June. The total microplankton assemblage structure was not altered in any lake by either enhanced grazing or enhanced nutrients, but they showed very distinct differences between lakes and by month.

Microplankton Assemblage Structure across Lakes of Varying Trophic State
According to the original Plankton Ecology Group (PEG) model of seasonality and succession of plankton assemblages in temperate lakes [67], phytoplankton biomass typically reaches peak biomass (i.e., a "bloom") in the spring, followed by a mid-summer low, and then a second bloom in late summer and early fall. The PEG model predicts that late summer blooms in eutrophic lakes are as large or larger in magnitude than the spring bloom, and are ultimately dominated by large and/or inedible phytoplankton taxa; in oligotrophic lakes the magnitude of the late summer bloom is much lower than in spring and is dominated by small phytoplankton taxa [67].
Typically, the onset of the spring bloom in temperate lakes is strongly linked to abiotic factors, including changes in light availability and water column mixing that makes subsurface nutrient stores available to phytoplankton in surface waters; whereas, the magnitude and decline of the late summer bloom is driven more by biotic (e.g., grazing and competition) factors [32,67,68]. These patterns were more comprehensively described in the "revised" PEG model [23], which predicted that both grazing and nutrient availability are dominant drivers of the late summer bloom in eutrophic lakes, but nutrient availability is typically the strongest factor influencing the late summer bloom in oligotrophic lakes.
We observed patterns of biomass and assemblage composition in our study lakes that generally align with the revised PEG model. In oligotrophic Cle Elum and Merwin lakes, phytoplankton biomass was consistently low with relatively muted peaks of diatom and flagellate biomass in late spring (May-June). There were two phytoplankton blooms in the eutrophic Lacamas Lake-a relatively small bloom during May and June dominated by diatoms and flagellates, and a second, larger bloom of mostly cyanobacteria in September. In the typically hypereutrophic Vancouver Lake, phytoplankton biomass was consistently higher than 50 µg chl a/L over the study period, with one very large (~150 µg chl a/L) bloom in September consisting of both cyanobacteria and diatoms. In both Lacamas and Vancouver lakes, the cyanobacterial blooms consisted of harmful filamentous genera (i.e., Dolichospermum spp. and Aphanizomenon flos-aquae). Notably, the magnitude of the phytoplankton bloom in Vancouver Lake during our study in 2014 was considerably lower than the average bloom size observed in this lake between 2007 and 2020 (~500 µg chl a/L; [26] and Rollwagen-Bollens unpublished data).

Effects of Enhanced Grazing and Enhanced Nutrients on Microplankton Assemblages
Neither of our treatments (i.e., enhanced grazers or enhanced nutrients) had a significant effect on microplankton assemblage structure within experimental incubations. These results suggest that grazer and nutrient enhancements, at the densities, concentrations, or ratios amended, did not alter selective pressures on microplankton taxa beyond those already present within each lake at the time of the experiments. However, microplankton assemblages were significantly different between lakes in each set of our experiments (June, August, October), indicating that there are strong differences in species composition across trophic state. Moreover, within each lake the plankton assemblages also differed over time, particularly between early summer (June) and late summer/fall (August and October).
Contrasts in physical factors such as water clarity, lake depth, elevation, and seasonal temperature ranges between our four study lakes may help explain why microplankton assemblage patterns appear to be strongly structured by differences associated with site (lake and, by extension, trophic state) and month. For instance, Lacamas is highly stratified in summer, which prevents mixing and allows the formation of bottom water hypoxia [49,50], while Vancouver is very shallow and well mixed, but highly turbid due to the resuspension of unconsolidated sediment at the lake bottom [46]. Both Cle Elum and Merwin are deep, and factors influencing euphotic zone mixing, such as winds and water column stability (not measured by us), are likely affecting phytoplankton dynamics in these lakes [31,32]. A study conducted in Sicilian lakes experiencing similar climate conditions noted that total phytoplankton biomass was connected to nutrient concentrations, but assembly was structured by euphotic depth, thermocline, and mixing depth [40]. Because site (lake) and month were strongly related to our observed assemblage patterns, the influence of light penetration and mixing dynamics on phytoplankton assemblage and bloom dynamics should be explored further in these systems, and regularly included in assemblage studies more broadly.

Effects of Enhanced Grazing ("Top Down") and Nutrients ("Bottom Up") on Phytoplankton Growth
Zooplankton grazer enhancement significantly reduced phytoplankton net growth in oligotrophic Cle Elum and oligo-mesotrophic Merwin in August, and again in Merwin in October. These results are in line with our predictions (Figure 1) and provide evidence that top-down control by daphnid cladocerans likely played a role in phytoplankton bloom decline in mid to late summer in Cle Elum and Merwin, lakes of relatively low trophic state. During these months the large-bodied cladoceran genus Daphnia was the most abundant mesozooplankton taxon in both Cle Elum and Merwin and was therefore the grazer taxon amended in those experiments. Daphnia are filter feeders, typically consuming organisms within a size range of~1 to 30 µm with minimal selectivity, unless prey taxa exhibit successful avoidance strategies or consumption is inhibited by prey size or shape [21,69]. Daphnia are known to consume and thrive on high-quality diets of diatoms and chrysophytes [69], which comprised a major portion of the chl-containing microplankton assemblage of both Cle Elum and Merwin in August and October. These results align with studies conducted in oligotrophic lakes in the Rocky Mountain region of Idaho (USA) [70] as well as oligotrophic Lake Tahoe (CA-NV, USA) and oligo-mesotrophic Castle Lake in the Siskiyou Mountains of northern CA (USA) [42], where the addition of mesozooplankton grazers had highly significant negative effects on phytoplankton chl a biomass and net growth rates.
In contrast, yet in keeping with our expectations, we found no evidence of a grazing effect by enhanced densities of copepods or daphnid cladocerans on phytoplankton growth in the eutrophic and hypereutrophic lakes (Lacamas and Vancouver, respectively). This lack of effect of amended grazing pressure was likely because in both lakes the dominant phytoplankton taxa were relatively inaccessible to the dominant mesozooplankton graz-ers, due to size range or potentially toxicity. For instance, in June, Lacamas experienced high abundance and biomass of the colonial diatom Fragillaria crotonensis, the large size of which may have reduced grazing effectiveness by Daphnia [71], potentially explaining why phytoplankton growth was not affected in our grazer-enhanced incubations from this lake. Similarly, as the summer progressed in both eutrophic Lacamas and Vancouver lakes, the high abundance and biomass of large filamentous cyanobacteria, Dolichospermum and Aphanizomenon, may have continued to limit grazer filtration rates, or caused selective consumption of non-filamentous and less abundant taxa, such as diatoms or chlorophytes [25,27]. This is supported by several previous studies in Vancouver Lake that show copepod grazing to be highly selective, with a preference for heterotrophic and mixotrophic protists and an avoidance of cyanobacteria [25,51]. In fact, copepods may be coexisting with, and even facilitating, dense filamentous cyanobacterial blooms that occur in Vancouver and other eutrophic lakes [26,51,72,73].
With respect to nutrient (bottom-up) factors, our results mostly aligned with our expectation that enhanced nutrient availability would stimulate higher phytoplankton growth in lakes on the low end of the trophic gradient, although we only observed a significant effect of added nutrients in oligotrophic Lake Merwin during June and none of the lakes during any other month. It is possible that we did not see a strong effect of added nutrients on phytoplankton net growth in most of our experimental incubations because nutrient amendments were not high enough to stimulate growth or that nutrients were not generally limiting in these lakes. However, it is likely that our nutrient (N and P) amendments were adequate, as treatment concentrations were frequently much higher than ambient (~1.5-100× P, and 4-450× N), with the exception of high N concentrations in Lacamas in June. In the Merwin experiment conducted in June, when nutrients did have a significant effect on phytoplankton growth, ambient concentrations of N and P were lower than in the other lakes in that month, including Cle Elum. Thus, June may have been a somewhat rare instance of nutrient limitation in Merwin (and possibly a rare case of higher nutrient availability in Cle Elum at that time). However, ambient nutrient availability did not appear to restrict phytoplankton growth during the late summer or fall in Merwin, nor in any of our other three lakes, such that control of phytoplankton bloom dynamics must be limited by other factors within these systems.
For instance, the taxonomic composition of the microplankton assemblage may have influenced the phytoplankton growth response to our experimental nutrient enrichments. Flagellates consistently comprised a major portion of the microplankton assemblage biomass in Cle Elum and Merwin compared to Lacamas and Vancouver. Chrysophyte flagellates are known to be common in clear-water lakes of lower trophic state, with many of these taxa exhibiting a preference for waters with low alkalinity [57,74]. Many chrysophyte genera are mixotrophic, having dual strategies for nutrient acquisition, and can consume bacteria [74,75] at significant rates, making them functionally similar to dinoflagellates and cryptomonads [76]. Mixotrophs can meet nutritional needs through phagotrophy and therefore are less likely to be nutrient-limited than fully autotrophic phytoplankton [75]. The positive phytoplankton growth response to our experimental nutrient enhancement in Merwin occurred in June, when autotrophic taxa (diatoms) comprised~70% of the total microplankton assemblage; however, in both Cle Elum and Merwin it was more common for mixotrophic taxa (flagellates and dinoflagellates) to comprise a greater proportion of the microplankton assemblage than autotrophic taxa. Thus, the high proportion of mixotrophic phytoplankton, which are less likely to be nutrient-limited than autotrophs, in Cle Elum and Merwin may in part explain why we saw a limited growth response to nutrient enhancement in experimental incubations in these two lakes of lower trophic state.

Conclusions
In summary, our results show that enhanced zooplankton grazers and enhanced nutrients affected phytoplankton net growth rates in study lakes that were more oligotrophic, but not in lakes that were eutrophic or hypereutrophic. Our results demonstrate that the effects of enhanced grazing and nutrients on phytoplankton bloom dynamics vary temporally within a given lake and across lakes spanning a gradient of eutrophication. Thus, we found that top-down (grazing) and bottom-up (nutrient) controls on phytoplankton growth are both important, but that the relative strength of these controls varies between systems and over time. These results further imply that the effects of bottom-up and top-down factors on phytoplankton dynamics vary by lake trophic state, and thus we suggest that future studies and predictive models of phytoplankton dynamics include a characterization of lake trophic state as a dependent ecosystem parameter.
Somewhat unexpectedly, enhanced grazer and nutrient treatments did not alter microplankton assemblage structure in our experiments, but the strong clustering patterns associated with month and site (lake) highlight the important influence that seasonal and other in-lake factors (not measured by us) can have in structuring microplankton assemblages. While nutrients are generally considered predictive of total phytoplankton biomass in a given system, other physical factors, such as light, mixing, and temperature, may more strongly regulate microplankton assemblage structure. Further, nutrients may not be the only or even the strongest factor limiting phytoplankton growth or structuring assembly and bloom dynamics in oligotrophic systems; grazing also needs to be considered.
We conclude that it is critically important to consider trophic state as well as the relative abundance of the dominant mesozooplankton taxa present when modeling or predicting the future response to the nutrient enrichment of lakes, particularly where assemblages contain a high proportion of mixotrophic taxa, or when considering increasing zooplankton abundance (i.e., biomanipulation) as a mechanism to reduce phytoplankton blooms.