Whole-Ecosystem Experiments Reveal Varying Responses of Phytoplankton Functional Groups to Epilimnetic Mixing in a Eutrophic Reservoir

Water column mixing can influence community composition of pelagic phytoplankton in lakes and reservoirs. Previous studies suggest that low mixing favors cyanobacteria, while increased mixing favors green algae and diatoms. However, this shift in community dominance is not consistently achieved when epilimnetic mixers are activated at the whole-ecosystem scale, possibly because phytoplankton community responses are mediated by mixing effects on other ecosystem processes. We conducted two epilimnetic mixing experiments in a small drinking water reservoir using a bubble-plume diffuser system. We measured physical, chemical, and biological variables before, during, and after mixing and compared the results to an unmixed reference reservoir. We observed significant increases in the biomass of cyanobacteria (from 0.8 ± 0.2 to 2.4 ± 1.1 μg L−1, p = 0.008), cryptophytes (from 0.7 ± 0.1 to 1.9 ± 0.6 μg L−1, p = 0.003), and green algae (from 3.8 to 4.4 μg L−1, p = 0.15) after our first mixing event, likely due to increased total phosphorus from entrainment of upstream sediments. After the second mixing event, phytoplankton biomass did not change but phytoplankton community composition shifted from taxa with filamentous morphology to smaller, rounder taxa. Our results suggest that whole-ecosystem dynamics and phytoplankton morphological traits should be considered when predicting phytoplankton community responses to epilimnetic mixing.


Introduction
Water column mixing can have a substantial impact on phytoplankton community composition in lakes and reservoirs [1][2][3].Specifically, previous research suggests that waterbodies with low mixing in the epilimnion are more likely to be dominated by cyanobacteria, while waterbodies with high mixing are often dominated by green algae and diatoms.This is because many cyanobacterial taxa are buoyant in the water column and can position themselves to take advantage of nutrients, thereby forming blooms which are disrupted and entrained by turbulent mixing [4].In contrast, green algae and diatoms are more dense and quick-sinking than cyanobacteria and benefit from mixing because it re-suspends them higher in the photic zone [1,5].Notably, these predicted responses of different phytoplankton groups to epilimnetic mixing are based primarily on laboratory, mesocosm, and modeling studies, rather than ecosystem experiments [6].
In keeping with these findings from smaller-scale studies, managers sometimes install bubble-plume epilimnetic mixing systems in lakes and reservoirs with the goal of shifting the dominance of the phytoplankton community from cyanobacteria, the phytoplankton taxon most responsible for blooms and toxin production in freshwaters, to less harmful taxa such as green algae and diatoms [3].However, at the whole-ecosystem scale, the efficacy of mixing systems in achieving this shift in community dominance is equivocal [3].Heo and Kim (2004) and Becker et al. (2006) reported that mixing systems caused a shift in community dominance from cyanobacteria to green algae and diatoms in Lake Dalbang, South Korea and Bleiloch Reservoir, Germany, respectively, although total phytoplankton biomass did not decrease in either case [7,8].However, the operation of mixing systems to reduce either cyanobacteria or total phytoplankton biomass has been unsuccessful in other cases.For example, mixers were unsuccessful in reducing cyanobacterial dominance or phytoplankton biomass in Chaffey Reservoir, Australia [9].
One possible explanation for inconsistent responses of phytoplankton to mixing at the whole-ecosystem scale may be that mixers are operated for different durations and at different intensities.Six of the 14 case studies examined in a recent review of artificial mixing for cyanobacterial control reported mixing systems that were operated intermittently [3], potentially due to the high cost of continuous operation, system maintenance, or other logistical constraints.This intermittent operation may decrease mixer efficacy, as reported for Lake Nieuwe Meer, The Netherlands [10].In addition, the depth and intensity of mixing may lead to unexpected, ecosystem-specific responses of phytoplankton to mixing.For example, some bubble-plume mixing systems are used to completely destratify the water column [9], while others are operated with the goal of maintaining a clearly defined epilimnion and hypolimnion [8].
The variable effects of whole-ecosystem mixing on phytoplankton community composition may also be attributable to the challenge of scaling up laboratory, mesocosm, and modeling approaches to whole ecosystems.At the whole-ecosystem scale, mixing has broader effects on lake and reservoir dynamics than simply the entrainment or resuspension of various phytoplankton taxa.Mixing could decrease the transparency of the water column by stirring up sediments and increasing turbidity in the water column (particularly in shallow waterbodies), leading to light limitation of photosynthesis and subsequent decrease in phytoplankton biomass [11,12].Alternatively, stirring of nutrient-rich sediments into the epilimnion could also provide a nutrient subsidy to phytoplankton growth [13,14].Mixing can also affect thermal stratification by deepening the thermocline and decreasing stratification strength, which could affect the optimal depth for phytoplankton growth [6,15].Furthermore, while our current understanding of mixing effects on phytoplankton is largely based on the relative buoyancy and density of aggregated phytoplankton taxonomic groups such as cyanobacteria and green algae, there may be other important phytoplankton traits that should be considered when predicting phytoplankton responses to mixing.Phytoplankton morphology-including the size, shape, and surface area to volume (SA:V) ratio-and the presence or absence of flagella have important implications for sinking rate and motility in the water column and are therefore highly relevant when considering phytoplankton responses to mixing [16][17][18].Within the green algae taxonomic group, taxa with round morphology, low SA:V ratios, and fast sinking rates would likely benefit from resuspension, whereas taxa with elongate morphologies and high SA:V ratios may not.While the use of functional traits to track community responses to disturbance is increasingly common in ecology [19,20], its application for management purposes has lagged, despite its ease and utility compared to traditional taxonomic methods.Considering a broad array of morphological traits when assessing phytoplankton responses to mixing may reveal changes that are not apparent based on aggregate measures of phytoplankton biomass [21].In sum, the current approach to understanding mixing effects on freshwater phytoplankton, which was largely developed using lab, mesocosm, and modeling studies and is based on aggregated phytoplankton groups, does not account for phytoplankton functional traits or mixing effects on ecosystem dynamics and how these may feed back to alter phytoplankton responses.To improve our understanding of phytoplankton responses to epilimnetic mixing at the whole-ecosystem scale, we conducted two pulsed mixing experiments in a small, eutrophic drinking water reservoir.We measured physical, chemical, and biological variables in both our experimental reservoir and a nearby reference reservoir before, during, and after mixing events to answer the following questions: (Q1) How does epilimnetic mixing affect the concentration of total phytoplankton and phytoplankton functional groups at the whole-ecosystem scale?and (Q2) How do the effects of epilimnetic mixing on physical and chemical variables (e.g., thermocline depth, light availability, and water column nutrient concentrations) mediate phytoplankton responses to mixing?
We found that epilimnetic mixing caused variable responses in phytoplankton biomass.Biomass across a suite of functional groups increased after our first mixing event but not after our second mixing event.This differential response was likely due to entrainment of upstream turbidity by the first mixing event, leading to a nutrient subsidy for phytoplankton in the lacustrine zone.After the second mixing event, despite no change in total biomass, functional group dominance shifted from taxa with large, colonial, filamentous morphology to taxa with medium/large, flagellated morphology.Phytoplankton with large, filamentous morphology decreased while small, siliceous taxa increased in the 24 h immediately after both mixing events, supporting previous hypotheses that mixing decreases cyanobacterial filaments while causing increases of dense, siliceous taxa such as diatoms.

Study Sites
We conducted epilimnetic mixing experiments in Falling Creek Reservoir (FCR; Figure 1), a small, eutrophic drinking water reservoir located in southwestern VA, USA (37 • 18 12" N, 79 • 50 14" W).FCR is owned and operated by the Western Virginia Water Authority (WVWA) and has a surface area of 0.119 km 2 and a maximum depth of 9.3 m at full pond.FCR is dimictic, with a summer thermal stratification period that generally lasts from April-October.The naturally-occurring summer thermocline depth is at 4-5 m [22].Mean water residence time during the study period in FCR was 68 ± 16 days (1 S.D.) [23].
FCR has an epilimnetic mixer (EM) that was installed by the WVWA to reduce the formation of cyanobacterial blooms.The EM is located at a depth of 5 m and extends 230 m up the reservoir from the deepest site near the dam (Figure 1).The EM consists of an onshore air compressor connected to a single diffuser line with five discrete sections of porous hose that release bubbles of atmospheric air into the water column at any rate from 0-25 SCFM (standard cubic feet per minute; see Text S1 for an in-depth description of the EM design).When activated, the EM causes turbulent mixing in the surface water of FCR from 0-5.5 m depth in the water column but does not destratify the entire water column, unlike many other mixing systems (Figure 2, Figure S1) [24].We focused on this surface layer of turbulent mixing, hereafter referred to as the "treatment zone," when assessing the effects of our epilimnetic mixing experiments.
Beaverdam Reservoir (BVR; Figure 1) serves as a reference ecosystem for FCR and is located 3 km upstream of FCR (37 • 18'59" N, 79 • 49'8" W).BVR is also owned by the WVWA and serves as a secondary drinking water source but does not have an EM system.BVR drains into FCR via a 1.7 km-long stream, which serves as the primary inflow to FCR [22].At full pond, BVR has a surface area of 0.39 km 2 and a maximum depth of 13 m, and is generally stratified from April-November with a typical thermocline depth of ~5-6 m [25].Throughout the study, BVR was maintained at a maximum depth of 11.5 m.Both reservoirs have primarily forested watersheds since abandonment of agricultural land in the 1930s [22].

Experimental Mixing Manipulations
In 2016, the EM in FCR was activated twice with different durations and mixing intensities during our monitoring period (12 May-26 August).Because these mixing experiments were conducted using a newly-installed mixing system in a drinking water reservoir from which managers were actively withdrawing water at the time of our experiments, we were limited to operating the mixing system for relatively short trial periods due to management concerns about the effects of mixing on water quality.In addition, because of the substantial planning, personnel, and logistical effort involved in conducting intensive sampling at two drinking water reservoirs over a short period, we planned the dates of our mixing experiments a priori for the time of year when cyanobacterial blooms in FCR had historically occurred, rather than adaptively activating the mixer during blooms that were developing.The first mixing event (hereafter, EM1) occurred on 30 May, when the mixer was activated from 10:00 to 16:00 at 25 SCFM (at maximum capacity); the second mixing event (EM2) occurred on 27-28 June, when the mixer was activated for 24 h (12:00 to 12:00) at 15 SCFM.conducted using a newly-installed mixing system in a drinking water reservoir from which managers were actively withdrawing water at the time of our experiments, we were limited to operating the mixing system for relatively short trial periods due to management concerns about the effects of mixing on water quality.In addition, because of the substantial planning, personnel, and logistical effort involved in conducting intensive sampling at two drinking water reservoirs over a short period, we planned the dates of our mixing experiments a priori for the time of year when cyanobacterial blooms in FCR had historically occurred, rather than adaptively activating the mixer during blooms that were developing.The first mixing event (hereafter, EM1) occurred on 30 May, when the mixer was activated from 10:00 to 16:00 at 25 SCFM (at maximum capacity); the second mixing event (EM2) occurred on 27-28 June, when the mixer was activated for 24 h (12:00 to 12:00) at 15 SCFM.The duration and flow rate of the two epilimnetic mixing events were chosen to preserve thermal stratification, as required by the reservoir managers.The two events were designed to be complementary, using different flow rates and durations to achieve the same change in metalimnetic boundary depth (a decrease of ~1 to 1.3 m) during both experiments (Text S1) [26].In addition, mixing was conducted with the goal of homogenizing phytoplankton biomass concentrations throughout the treatment zone (0-5.5 m depth) at the deepest site of our experimental reservoir, which we determined by real-time monitoring of phytoplankton biomass depth distributions during our mixing experiments.The duration and flow rate of the two epilimnetic mixing events were chosen to preserve thermal stratification, as required by the reservoir managers.The two events were designed to be complementary, using different flow rates and durations to achieve the same change in metalimnetic boundary depth (a decrease of ~1 to 1.3 m) during both experiments (Text S1) [26].In addition, mixing was conducted with the goal of homogenizing phytoplankton biomass concentrations throughout the treatment zone (0-5.5 m depth) at the deepest site of our experimental reservoir, which we determined by real-time monitoring of phytoplankton biomass depth distributions during our mixing experiments.

Field Sampling
Throughout the monitoring period, we collected depth profiles of physical, chemical, and biological variables one to three times per week at the deepest sites in FCR and BVR (Figure 1).Profiles were obtained immediately before mixing (before 10:00 on the day the mixer was activated), immediately after mixing (at 18:00 on 30 May for EM1 and at 13:00 on 28 Jun for EM2), 24 h post-mixing, and 48-72 h post-mixing in FCR.In BVR, profiles were obtained 24 h before mixing in FCR, 24 h after mixing in FCR, and 48-72 h post-mixing in FCR.We used a 4-Hertz Conductivity, Temperature, and Depth (CTD) profiler (SeaBird Electronics, Bellevue, WA, USA) equipped with a turbidity sensor (WET Labs ECO fluorometer, SeaBird Scientific, Bellevue, WA, USA) to collect high-resolution (~1 cm) profiles of temperature and turbidity [27].Profiles at 1-m resolution of photosynthetically active radiation (PAR) were collected with an LI-250 underwater light meter (LI-COR, Lincoln, NV, USA) [28].We also collected depth profiles of total nitrogen (TN), total phosphorus (TP), nitrate (NO 3 − ), ammonium (NH 4 + ), and soluble reactive phosphorus (SRP) concentrations using a 4-L van Dorn sampler (Wildco, Yulee, FL, USA) at the depths of the outflow valves of FCR (0.1, 1.6, 3.8, 5, 6.2, 8, and 9 m) and evenly spaced on the vertical profile of BVR (0.1, 3, 6, 9, and 11 m) [29].Dissolved samples for NO 3 − , NH 4 + , and SRP analysis were filtered through Whatman GF/F filters immediately after collection, and both dissolved and total samples were frozen until analysis.Finally, we collected CTD and FluoroProbe profiles once a week at four sites in FCR upstream of the deepest site in the riverine zone of the reservoir (Figure 1).
To assess phytoplankton responses to mixing, we used a FluoroProbe (bbe moldaenke, Schwentinental, Germany) to collect high-resolution profiles (~10 cm) of phytoplankton biomass at each site [30,31].The FluoroProbe is a submersible in situ fluorometer that captures pigment signatures of four spectral groups that approximately correlate to the biomass of aggregated diatoms, dinoflagellates, and chrysophytes (hereafter, brown algae); green algae; cyanobacteria; and cryptophytes [32].FluoroProbe sensors are often used to track the concentration and vertical distribution of these aggregated phytoplankton groups in lakes and reservoirs [30,33].The FluoroProbe was calibrated for site-specific use at FCR and BVR via a yellow substances offset that corrects for any spectral changes due to dissolved organic matter in the reservoir [32].We also collected water samples at the same depths as our nutrient samples in the treatment zones of both FCR and BVR that were filtered through Whatman GF/C filters for chlorophyll-a analysis to verify total phytoplankton biomass as reported by the FluoroProbe.
To assess functional group responses within the phytoplankton assemblage, 250 mL of water were collected with a van Dorn sampler at the depth of the fluorescence maximum as determined by the FluoroProbe.These samples were collected from FCR on every sampling day from two weeks before EM1 until two weeks after EM2 (16 May-11 July) and from BVR two weeks before and after EM2 (23 June-7 July), and immediately preserved with Lugol's iodine solution for phytoplankton microscopy.
We also monitored wind and precipitation using a research-grade meteorological station installed on the dam at FCR (Campbell Scientific, Logan, UT, USA).

Laboratory Analyses
In the laboratory, all water chemistry samples were processed following standard methods.TN and TP samples were digested with alkaline persulfate [34] and then analyzed colorimetrically using flow injection analysis (Lachat ASX 520 Series, Lachat Instruments, Loveland, CO, USA).SRP, NO 3 − , and NH 4 + were also analyzed colorimetrically using flow injection analysis on the Lachat [35].
Chlorophyll-a was analyzed using Standard Method 10200H [35] on a UV-Vis Spectrophotometer (UV-Vis 1601, Shimadzu Scientific Instruments, Durham, NC, USA).Pelagic phytoplankton samples from the depth of fluorescence maximum of FCR were identified and counted using a Sedgwick-Rafter cell (Wildco, Yulee, FL, USA) at 400× resolution on a Nikon Eclipse TS100 inverted microscope (Nikon, Tokyo, Japan).At least 300 cells per sample were counted and identified to the genus level.The biovolumes of the first 10 individuals in each genus were calculated by approximating the cells to known geometric shapes using measured linear dimensions [36].Picophytoplankton, or phytoplankton with a maximum linear dimension ≤2 µm, were not counted.All microscopy was conducted by the same taxonomist (M.E.L.).

Calculation of Thermal Stability Metrics and Light Attenuation
To assess thermal stratification in FCR and BVR, we calculated several metrics of lake thermal structure (buoyancy frequency, Schmidt stability, lake number, parent thermocline depth, and lower metalimnetic boundary depth) at the deepest site of each reservoir for each sampling day using the LakeAnalyzer program [37] in MATLAB software (v.R2016a 9.0.0.341360,MathWorks, Natick, MA, USA).LakeAnalyzer defines parent thermocline depth as either the same as the thermocline depth (depth of greatest water density gradient) or as the depth of a secondary local density gradient maximum if (1) it is deeper than the thermocline depth and (2) its magnitude is at least 20% of the thermocline depth [38].The attenuation coefficient (K d ) for photosynthetically active radiation (PAR) was calculated as the slope of the best fit line for the natural logarithm of PAR (measured in µmol m −2 s −1 ) plotted against depth (m) for each sample day [39].

Statistical Analyses
To assess experimental mixing manipulation effects, we conducted a before-after-control-impact (BACI) analysis on the physical, chemical, and biological variables measured in our experimental and reference reservoirs [40].To focus our analysis on the layer directly affected by the EM in FCR, we restricted our analysis to the treatment zone in FCR (0-5.5 m) and an equivalent layer depth in BVR (0-6 m).Two treatment zone response variables were analyzed for each of the four aggregate phytoplankton groups measured by the FluoroProbe on every sampling day: mean biomass concentration in the treatment zone at the deepest site of each reservoir (using the hypsographic data for the reservoir, Figure 1) and depth of maximum fluorescence within the treatment zone.We examined the depth of maximum fluorescence (F max ) because many water bodies exhibit deep chlorophyll maxima, and changes in F max depth could indicate a phytoplankton response to changes in thermocline depth, light availability, or nutrient availability within the water column as a result of mixing [15,41,42].
Our BACI analysis was a Welch's t-test on the average difference between measurements at FCR and BVR summed over the two-week period before mixing versus the two-week period after mixing [43].We chose a two-week window because phytoplankton spectral groups were largely able to recover their midday pre-mixing depth distributions within two weeks.Volume-weighted nutrient concentrations across the treatment zone before and after mixing and the lake physics parameters were also analyzed using a BACI analysis during the two weeks before versus two weeks after mixing.For volume-weighting, we divided the treatment zone of FCR and BVR into layers with each discrete-depth nutrient sample occurring in the middle of a layer.In cases where nutrient samples were unevenly distributed on the treatment zone depth profile, we divided the depth between nutrient samples equally between the shallower sample and the deeper sample (e.g., since nutrient samples were taken at 0.1, 1.6, 3.8, and 5 m in FCR and the treatment zone was 0-5.5 m, our volume-weighting layers in the treatment zone were 0.1-0.8m, 0.9-2.7 m, 2.8-4.4 m, and 4.5-5.5 m).Our nutrient sample concentrations were then weighted based on the relative volume of the water layer in which they occurred using a detailed bathymetric profile of each reservoir (0.1 m resolution) and used to calculate the average nutrient concentration throughout the treatment zone of each reservoir on each sampling day.All analyses were conducted in the R statistical environment (R version 3.4.3,R Core Development Team, 2017).

Relative Abundance of Phytoplankton Functional Groups
To assess the relative abundance of phytoplankton functional groups, which provides complementary information to the FluoroProbe profiles on phytoplankton responses to mixing, we followed the morphology-based functional group (MBFG) classification developed by Kruk et al. (2010) [16].In this classification system, phytoplankton are categorized by three quantitative variables (maximum linear dimension, cell biovolume, and SA:V ratio) and four presence/absence variables (presence or absence of flagella, mucus secretion, intracellular gas vesicles, and siliceous exoskeletons) into seven MBFGs [16,17].To assign our identified genera at each study site to MBFGs, we took the mean value of the quantitative variables measured for each genus over our study period and assigned presence/absence traits based on whether the majority of observed natural units (cells or colonies) of that genus possessed the trait.We used an open-access R script provided in the supplemental material of Kruk and Segura (2012) to assign MBFG classification for each genus at each study site based on observed traits (Table S1) [17].This classification was then used to calculate the relative abundance of each MBFG at each study site over time.

Thermal Structure
The goal of our study was to conduct two pulsed mixing experiments to determine the response of phytoplankton to epilimnetic mixing at the whole-ecosystem scale.Here, we report on the thermal structure of our experimental reservoir relative to the reference reservoir to determine whether our mixing manipulations caused changes in the thermal structure of FCR.
The parent thermocline in FCR ranged from 3.3-7.0m due to mixing manipulations (Figure 2).The first and second mixing experiments decreased the thermocline depth at FCR's deepest site by 1.0 m (p = 0.05) and 0.3 m (p = 0.008), respectively.The lower metalimnetic boundary in FCR decreased by 0.8 m after the first mixing event (EM1) and by 0.4 m after the second mixing event (EM2; Figure 2), though these changes were not significant compared to BVR (both p > 0.11).Mixing also affected thermal structure in the region beyond the diffuser in FCR, suggesting that the patterns observed at the deepest site were representative of much of the reservoir (Figure S1).(maximum linear dimension, cell biovolume, and SA:V ratio) and four presence/absence variables (presence or absence of flagella, mucus secretion, intracellular gas vesicles, and siliceous exoskeletons) into seven MBFGs [16,17].To assign our identified genera at each study site to MBFGs, we took the mean value of the quantitative variables measured for each genus over our study period and assigned presence/absence traits based on whether the majority of observed natural units (cells or colonies) of that genus possessed the trait.We used an open-access R script provided in the supplemental material of Kruk and Segura (2012) to assign MBFG classification for each genus at each study site based on observed traits (Table S1) [17].This classification was then used to calculate the relative abundance of each MBFG at each study site over time.

Thermal Structure
The goal of our study was to conduct two pulsed mixing experiments to determine the response of phytoplankton to epilimnetic mixing at the whole-ecosystem scale.Here, we report on the thermal structure of our experimental reservoir relative to the reference reservoir to determine whether our mixing manipulations caused changes in the thermal structure of FCR.
The parent thermocline in FCR ranged from 3.3-7.0m due to mixing manipulations (Figure 2).The first and second mixing experiments decreased the thermocline depth at FCR's deepest site by 1.0 m (p = 0.05) and 0.3 m (p = 0.008), respectively.The lower metalimnetic boundary in FCR decreased by 0.8 m after the first mixing event (EM1) and by 0.4 m after the second mixing event (EM2; Figure 2), though these changes were not significant compared to BVR (both p > 0.11).Mixing also affected thermal structure in the region beyond the diffuser in FCR, suggesting that the patterns observed at the deepest site were representative of much of the reservoir (Figure S1).Schmidt stability in FCR exhibited a statistically significant decrease after EM2 relative to BVR (p = 0.003; Figure 3, Table 1).While Schmidt stability did not decrease in FCR after EM1, BVR Schmidt stability in FCR exhibited a statistically significant decrease after EM2 relative to BVR (p = 0.003; Figure 3, Table 1).While Schmidt stability did not decrease in FCR after EM1, BVR exhibited a sharp increase in Schmidt stability in the two weeks after EM1-as air temperature warmed-that was not observed in FCR (Figure 3).Thus, both mixing events altered the thermal structure of FCR relative to BVR (Table 1).There was no statistically significant difference in mean treatment zone temperature at the deepest site in FCR relative to BVR after either mixing event (Table 1; both p ≥ 0.76).
The meteorological station at FCR did not record any large rain or wind storms during the monitoring period, indicating that changes in thermal stability were due to the mixing manipulations rather than weather events (Figure S2).Furthermore, because changes in parent thermocline depth post-mixing persisted for weeks after each mixing manipulation, changes in parent thermocline observed immediately after mixing events were unlikely to be solely the result of convective mixing due to evaporation and cooling at night.exhibited a sharp increase in Schmidt stability in the two weeks after EM1-as air temperature warmed-that was not observed in FCR (Figure 3).Thus, both mixing events altered the thermal structure of FCR relative to BVR (Table 1).There was no statistically significant difference in mean treatment zone temperature at the deepest site in FCR relative to BVR after either mixing event (Table 1; both p ≥ 0.76).
The meteorological station at FCR did not record any large rain or wind storms during the monitoring period, indicating that changes in thermal stability were due to the mixing manipulations rather than weather events (Figure S2).Furthermore, because changes in parent thermocline depth post-mixing persisted for weeks after each mixing manipulation, changes in parent thermocline observed immediately after mixing events were unlikely to be solely the result of convective mixing due to evaporation and cooling at night.

Aggregate Phytoplankton Groups
To assess phytoplankton response to epilimnetic mixing, we measured phytoplankton using fluorescence-based depth profiles before, during, and after mixing in both reservoirs.Measures of total phytoplankton biomass using the FluoroProbe were well-correlated with chlorophyll-a concentrations measured using standard methods in the epilimnia of FCR and BVR throughout the study period (Figure S3).Prior to the EM1 mixing, aggregate groups exhibited broadly similar patterns of vertical distribution and concentration across the treatment zone of the deepest sites of FCR and BVR, supporting the use of BVR as a reference system for FCR.These data allowed us to assess changes in the vertical distribution and total biomass of aggregate phytoplankton groups in the two weeks after each mixing event.
Aggregate phytoplankton groups (green algae, cyanobacteria, brown algae, and cryptophytes) exhibited substantial changes in their vertical distribution in FCR during the two weeks post-mixing (Figure 4, Figures S4 and S5).One hour after each mixing event, phytoplankton concentrations were homogenized throughout the water column in FCR, but over the course of the two weeks postmixing, all aggregate groups re-formed fluorescence maxima in the treatment zone (Figure 4, Figures S4 and S5).Only green algae exhibited a statistically significant change in Fmax depth in the two weeks after EM1, rising from 5.4 ± 0.2 m (1 S.D.) pre-mixing to 3.4 ± 1.4 m post-mixing.No aggregate groups exhibited significant Fmax depth changes after EM2 (Table 1).

Aggregate Phytoplankton Groups
To assess phytoplankton response to epilimnetic mixing, we measured phytoplankton using fluorescence-based depth profiles before, during, and after mixing in both reservoirs.Measures of total phytoplankton biomass using the FluoroProbe were well-correlated with chlorophyll-a concentrations measured using standard methods in the epilimnia of FCR and BVR throughout the study period (Figure S3).Prior to the EM1 mixing, aggregate groups exhibited broadly similar patterns of vertical distribution and concentration across the treatment zone of the deepest sites of FCR and BVR, supporting the use of BVR as a reference system for FCR.These data allowed us to assess changes in the vertical distribution and total biomass of aggregate phytoplankton groups in the two weeks after each mixing event.
Aggregate phytoplankton groups (green algae, cyanobacteria, brown algae, and cryptophytes) exhibited substantial changes in their vertical distribution in FCR during the two weeks post-mixing (Figure 4, Figures S4 and S5).One hour after each mixing event, phytoplankton concentrations were homogenized throughout the water column in FCR, but over the course of the two weeks post-mixing, all aggregate groups re-formed fluorescence maxima in the treatment zone (Figure 4, Figures S4 and S5).Only green algae exhibited a statistically significant change in F max depth in the two weeks after EM1, rising from 5.4 ± 0.2 m (1 S.D.) pre-mixing to 3.4 ± 1.4 m post-mixing.No aggregate groups exhibited significant F max depth changes after EM2 (Table 1).Mean biomass concentrations of aggregated phytoplankton groups responded differently to EM1 and EM2 (Figure 5).Specifically, EM1 stimulated cyanobacteria and cryptophytes (both p ≤ 0.008; Table 1, Figure 5C,G) and also led to an increase in green algae (p = 0.15, Figure 5A).Increases in these aggregate groups summed to a 36% increase in mean phytoplankton biomass in the treatment zone of FCR after mixing (p = 0.04, Figure 5I).In contrast, green algae in FCR exhibited a significant decrease after EM2 relative to BVR (p = 0.0002; Figure 5B), while cyanobacterial and cryptophyte biomass did not change (both p ≥ 0.69), resulting in a 12% decrease in mean phytoplankton biomass in FCR.While the BACI analysis for mean total biomass across the treatment zone after EM2 suggests that this decrease is statistically significant (p = 0.02), there was also a 37% increase in phytoplankton biomass in BVR shortly after EM2.Consequently, the statistically significant response of total treatment zone biomass after EM2 is likely driven by both decreased biomass in the experimental reservoir and increased biomass in the reference reservoir.Brown algae did not exhibit changes in volume-weighted treatment zone biomass after either mixing event (both p ≥ 0.23, Figure 5E,F).
Upstream sites in FCR also experienced an increase in mean phytoplankton biomass in the surface waters after EM1 (Figure 6), suggesting that the phytoplankton response to mixing at the deepest site was representative of the rest of the reservoir.In addition, the first mixing event entrained biomass from the furthest upstream site along the bottom of the reservoir and downstream to Sites 2 and 3 (Figure 1), resulting in an increase in phytoplankton biomass in the metalimnion of the transitional zone of the reservoir within 24 h after mixing (Figure 6).This entrainment pattern was repeated after EM2.

Morphology-Based Phytoplankton Groups
To complement our fluorescence-based profiles of phytoplankton biomass, we collected samples from the depth of maximum fluorescence throughout the monitoring period in FCR to assess changes in morphology-based phytoplankton functional groups after mixing.These data allow us to address the hypothesis that epilimnetic mixing should cause decreases in colonial, filamentous taxa while stimulating dense, round, quick-sinking taxa.Total phytoplankton biovolume at the depth of F max calculated from microscopy exhibited generally the same pattern as the total phytoplankton concentration reported by the FluoroProbe with some variability, as expected, because of the inherent differences in the methodologies (Figure S6).In particular, all of the FluoroProbe data we report here are averaged across the treatment zone (0-5.5 m) in FCR, while the functional group samples are point samples taken at the depth of the fluorescence maximum only.
Before both mixing events, the fluorescence maximum in FCR was dominated by Dolichospermum spp.(formerly known as Anabaena spp.), a large, colonial, filamentous cyanobacterium (Figure 7B).After each mixing event, phytoplankton samples analyzed microscopically revealed substantial variability in morphology-based functional group (MBFG) responses (Figure 7).After EM1, five of the seven morphology-based functional groups experienced an immediate decline followed by a steep increase in biovolume at the Fmax depth in FCR.Overall, this led to a substantial increase in total biovolume at the Fmax depth one week after EM1 (Figure 7A).This pattern of an immediate post-mixing decline followed by a substantial increase was most noticeable in large, high SA:V filamentous taxa (MBFG 3), which corresponds to large, colonial, filamentous cyanobacteria (Dolichospermum) in our reservoirs (Table S1).MBFG 3 declined by four orders of magnitude immediately after mixing but then rebounded and continued to increase for a net gain in biovolume one week after mixing.The other four MBFGs which exhibited this pattern included small, high SA:V phytoplankton (MBFG 1), medium-sized phytoplankton with no special traits (MBFG 4), medium/large flagellated phytoplankton (MBFG 5), and large, mucilaginous, low SA:V colonies (MBFG 7).These four MBFGs correspond to a variety of green algae and golden algae taxa, ranging from taxa with elongated morphologies and no flagella (such as Ankistrodesmus) to flagellated taxa with round morphologies (such as Nephroselmis) to colonial taxa with irregular shapes (such as Dinobryon).After EM1, five of the seven morphology-based functional groups experienced an immediate decline followed by a steep increase in biovolume at the F max depth in FCR.Overall, this led to a substantial increase in total biovolume at the F max depth one week after EM1 (Figure 7A).This pattern of an immediate post-mixing decline followed by a substantial increase was most noticeable in large, high SA:V filamentous taxa (MBFG 3), which corresponds to large, colonial, filamentous cyanobacteria (Dolichospermum) in our reservoirs (Table S1).MBFG 3 declined by four orders of magnitude immediately after mixing but then rebounded and continued to increase for a net gain in biovolume one week after mixing.The other four MBFGs which exhibited this pattern included small, high SA:V phytoplankton (MBFG 1), medium-sized phytoplankton with no special traits (MBFG 4), medium/large flagellated phytoplankton (MBFG 5), and large, mucilaginous, low SA:V colonies (MBFG 7).These four MBFGs correspond to a variety of green algae and golden algae taxa, ranging from taxa with elongated morphologies and no flagella (such as Ankistrodesmus) to flagellated taxa with round morphologies (such as Nephroselmis) to colonial taxa with irregular shapes (such as Dinobryon).Two of the morphology-based functional groups increased immediately after EM1; these were small siliceous flagellates (MBFG 2), which correspond most closely to dinoflagellates and other taxa within the brown algae aggregate group, and non-flagellated, siliceous phytoplankton (MBFG 6), which correspond most closely to diatoms (Figure 7A, Table S1).Of all the MBFGs, only small, siliceous flagellates did not experience an overall increase in biovolume in the week after EM1; though they increased by five orders of magnitude immediately after the mixing event, they subsequently declined to pre-mixing levels (Figure 7A).These changes in biovolume resulted in a temporary shift in community dominance at the Fmax from large, filamentous cyanobacteria (MBFG 3) to medium-sized phytoplankton with no special traits and medium/large flagellates (MBFGs 4 and 5).However, MBFG 3 regained community dominance one week after EM1 (Figure 7B).
Four of the five MBFGs that declined immediately after EM1 also declined after EM2, but unlike EM1, most of these groups did not undergo a subsequent rebound in biomass at the Fmax depth (Figure 7A).Groups that declined immediately after both mixing events included small, high SA:V phytoplankton (MBFG 1), large, high SA:V filamentous taxa (MBFG 3), medium-sized phytoplankton with no special traits (MBFG 4), and medium/large flagellated phytoplankton (MBFG 5).Nonflagellated, siliceous phytoplankton (MBFG 6), which increased after the first mixing event, declined after the second mixing event, while large, mucilaginous, low SA:V colonies (MBFG 7), which had decreased after the first mixing event, increased after EM2.Small, siliceous flagellates (MBFG 2) were the only group to increase after both mixing events.In sum, changes in biovolume among the MBFGs resulted in a slight decrease in total biovolume at the depth of fluorescence maximum after EM2, unlike the increase observed after EM1 (Figure 7A).As in EM1, these changes in biovolume resulted in a community dominance shift from large, filamentous cyanobacteria (MBFG 3) to flagellates (MBFG 5) immediately after EM2; however, unlike EM1, MBFG 3 never regained community dominance in the two weeks after EM2 (Figure 7B).Two of the morphology-based functional groups increased immediately after EM1; these were small siliceous flagellates (MBFG 2), which correspond most closely to dinoflagellates and other taxa within the brown algae aggregate group, and non-flagellated, siliceous phytoplankton (MBFG 6), which correspond most closely to diatoms (Figure 7A, Table S1).Of all the MBFGs, only small, siliceous flagellates did not experience an overall increase in biovolume in the week after EM1; though they increased by five orders of magnitude immediately after the mixing event, they subsequently declined to pre-mixing levels (Figure 7A).These changes in biovolume resulted in a temporary shift in community dominance at the F max from large, filamentous cyanobacteria (MBFG 3) to medium-sized phytoplankton with no special traits and medium/large flagellates (MBFGs 4 and 5).However, MBFG 3 regained community dominance one week after EM1 (Figure 7B).
Four of the five MBFGs that declined immediately after EM1 also declined after EM2, but unlike EM1, most of these groups did not undergo a subsequent rebound in biomass at the F max depth (Figure 7A).Groups that declined immediately after both mixing events included small, high SA:V phytoplankton (MBFG 1), large, high SA:V filamentous taxa (MBFG 3), medium-sized phytoplankton with no special traits (MBFG 4), and medium/large flagellated phytoplankton (MBFG 5).Non-flagellated, siliceous phytoplankton (MBFG 6), which increased after the first mixing event, declined after the second mixing event, while large, mucilaginous, low SA:V colonies (MBFG 7), which had decreased after the first mixing event, increased after EM2.Small, siliceous flagellates (MBFG 2) were the only group to increase after both mixing events.In sum, changes in biovolume among the MBFGs resulted in a slight decrease in total biovolume at the depth of fluorescence maximum after EM2, unlike the increase observed after EM1 (Figure 7A).As in EM1, these changes in biovolume resulted in a community dominance shift from large, filamentous cyanobacteria (MBFG 3) to flagellates (MBFG 5) immediately after EM2; however, unlike EM1, MBFG 3 never regained community dominance in the two weeks after EM2 (Figure 7B).Functional group concentrations in the Fmax of BVR, which experienced no mixing events, were less variable throughout the EM2 monitoring period of 23 June to 7 July than in FCR (Figure S7).In contrast to FCR, which experienced shifts in the dominant functional groups in the Fmax after mixing events, medium/large flagellated phytoplankton (MBFG 5) was the dominant functional group in the Functional group concentrations in the F max of BVR, which experienced no mixing events, were less variable throughout the EM2 monitoring period of 23 June to 7 July than in FCR (Figure S7).In contrast to FCR, which experienced shifts in the dominant functional groups in the F max after mixing events, medium/large flagellated phytoplankton (MBFG 5) was the dominant functional group in the F max of BVR throughout the EM2 monitoring period, with proportional biovolumes ranging from 0.57-0.92(Figure S7).

Sediment and Nutrient Dynamics
We collected a suite of environmental variables in both FCR and BVR before, during, and after our mixing events to determine whether the effects of epilimnetic mixing on whole-ecosystem dynamics might feed back to affect the phytoplankton response to mixing.We hypothesized that inconsistent responses of phytoplankton responses to mixing at the whole-ecosystem scale might be due to unforeseen effects of mixing on variables such as nutrient concentrations, turbidity, and light.
There were no significant changes in nutrient concentrations in FCR after EM2 (all p > 0.1; Table 1).Median turbidity in the treatment zone at the deepest site in FCR increased after EM1 by 0.4 NTU (p = 0.001, Figure 3), but not after EM2 (p = 0.99).Median turbidity at upstream sites 1 and 2 in FCR was 0.6-0.7 NTU higher than at the deepest site during the EM1 monitoring period of 16 May-6 June (Figure 9).Some of this higher-turbidity water was subsequently entrained downstream during the week following EM1 (Figure 9).Despite a small increase in water column turbidity at the deepest site after EM1, the photosynthetically active radiation (PAR) attenuation coefficient (K d ) did not significantly change in FCR after either EM1 or EM2 (both p ≥ 0.28, Table 1).Fmax of BVR throughout the EM2 monitoring period, with proportional biovolumes ranging from 0.57-0.92(Figure S7).

Sediment and Nutrient Dynamics
We collected a suite of environmental variables in both FCR and BVR before, during, and after our mixing events to determine whether the effects of epilimnetic mixing on whole-ecosystem dynamics might feed back to affect the phytoplankton response to mixing.We hypothesized that inconsistent responses of phytoplankton responses to mixing at the whole-ecosystem scale might be due to unforeseen effects of mixing on variables such as nutrient concentrations, turbidity, and light.
Median turbidity in the treatment zone at the deepest site in FCR increased after EM1 by 0.4 NTU (p = 0.001, Figure 3), but not after EM2 (p = 0.99).Median turbidity at upstream sites 1 and 2 in FCR was 0.6-0.7 NTU higher than at the deepest site during the EM1 monitoring period of 16 May-6 June (Figure 9).Some of this higher-turbidity water was subsequently entrained downstream during the week following EM1 (Figure 9).Despite a small increase in water column turbidity at the deepest site after EM1, the photosynthetically active radiation (PAR) attenuation coefficient (Kd) did not significantly change in FCR after either EM1 or EM2 (both p ≥ 0.28, Table 1).

Discussion
Contrary to theoretical expectations, yet in keeping with previous reports of inconsistent responses of phytoplankton to epilimnetic mixing and thermocline deepening, total phytoplankton biomass and cyanobacteria increased after the first whole-reservoir mixing experiment (EM1), while there was no change in total biomass after the second mixing experiment (EM2).Below, we explore the responses of various aggregated phytoplankton spectral groups and morphology-based

Discussion
Contrary to theoretical expectations, yet in keeping with previous reports of inconsistent responses of phytoplankton to epilimnetic mixing and thermocline deepening, total phytoplankton biomass and cyanobacteria increased after the first whole-reservoir mixing experiment (EM1), while there was no change in total biomass after the second mixing experiment (EM2).Below, we explore the responses of various aggregated phytoplankton spectral groups and morphology-based phytoplankton functional groups, as well as the possible ecosystem drivers of the different responses between EM1 and EM2.
Like all organisms, phytoplankton interact with their environment via their traits [19].We chose to assess phytoplankton mixing response using morphology-based functional groups because the traits in this classification system (size, shape, buoyancy, density, and motility) are highly relevant to how phytoplankton experience water column mixing and are also commonly used in other studies of phytoplankton functional traits [18].
For both of our mixing experiments, large, high SA:V, filamentous taxa (MBFG 3) were the dominant functional group at the depth of chlorophyll maximum in our experimental reservoir immediately prior to mixing.This group, along with small siliceous flagellated taxa (MBFG 2), were the most responsive to mixing events in the 24 h after mixing, although they displayed opposite patterns: large filamentous taxa decreased immediately after mixing, while small siliceous flagellated taxa increased (Figure 7).This finding agrees with previous work suggesting that bloom-forming, filamentous cyanobacteria would lose their buoyancy and have their filaments broken apart by turbulence during mixing, while denser siliceous taxa would be stimulated by resuspension in the photic zone [1,7,15].However, after EM1, the immediate decrease in filamentous taxa and increase in small siliceous taxa at the deepest site was quickly subsumed by an increase in biovolume of five of the seven MBFGs in the week after mixing.The MBFGs that experienced a net increase in biovolume in the week after EM1 include taxa with a wide range of fluorescence signatures, supporting the increases in average treatment zone biomass of green algae, cyanobacteria, and cryptophytes reported by the FluoroProbe after EM1 (Figures 5 and 7A, Table S1).Interestingly, although both MBFGs corresponding to brown algae (MBFGs 2 and 6) increased in biovolume after EM1, our FluoroProbe BACI test reported no change in biomass of brown algae during the EM1 monitoring period (Figure 5).This discrepancy may be because small siliceous flagellates (MBFG 2) decreased during the EM1 monitoring period, while non-flagellated, siliceous phytoplankton (MBFG 6) increased during the EM1 monitoring period (Figure 7A), resulting in a compensation effect between MBFG 2 and MBFG 6 [21].
Whole-ecosystem experiments are inherently fraught with difficulty, as they are challenging to replicate.However, we argue that we can still derive important insights from whole-ecosystem manipulations that cannot be deduced from laboratory or mesocosm experiments, provided that we proceed cautiously with our interpretation of results [44][45][46].Keeping this caveat in mind, we posit that the increase in biomass and biovolume across multiple spectral groups and MBFGs in FCR after EM1, which was not observed after EM2, may be due to entrainment of turbidity from upstream regions of FCR to the deepest site after EM1.
The two most salient differences in our environmental variables between EM1 and EM2 are an increase in total phosphorus and turbidity at the deepest site in FCR after EM1 that was not observed in our reference reservoir after EM1 and did not occur in either reservoir after EM2.As shown in Figure 9, the increase in turbidity at the deepest site in FCR after EM1 likely originated in shallow, turbid upstream regions of the reservoir and was subsequently entrained downstream over the course of the week post-EM1.This possibility is supported by 3-D hydrodynamic modeling of the same reservoir by Chen et al. (2018).Using field turbidity data, Chen et al. (2018) found that the intensity of EM operation in FCR was an important driver of short-term particle transport, and that mixing facilitates the transport of particulate matter from the shallow upstream region to the deep lacustrine region in the reservoir [26].The increased turbidity at the deepest site in FCR after the first mixing event coincided with an increase in volume-weighted TP throughout the treatment zone, as well as an increase in volume-weighted green algae, cyanobacteria, and cryptophyte biomass and total phytoplankton biomass in the treatment zone of FCR.None of these responses were observed in our reference reservoir throughout the monitoring period or in FCR after EM2.Taken together, these post-EM1 responses suggest that the increased total phytoplankton biomass in the treatment zone after the first mixing event may have been due to entrainment of suspended sediment from the upstream riverine zone of FCR into the deepest site, providing a nutrient subsidy that enabled an increase of green algae, cyanobacteria, cryptophyte, and overall total phytoplankton biomass.The possibility of a turbidity-induced nutrient subsidy is supported by three pieces of evidence: changes in nutrient concentrations post-EM1 at the depth of chlorophyll maximum in FCR, a back-of-the-envelope calculation of nutrient uptake rates at the depth of chlorophyll maximum post-EM1, and the difference in seasonal timing and intensity between EM1 and EM2, as described below.
First, in addition to an increase in total phosphorus across the treatment zone in FCR after EM1, there was also a decrease in SRP and the ratio of SRP:TP at 3.8 m.This depth (3.8 m) corresponds closely to the Fmax depth of phytoplankton biomass in the two weeks after EM1 (4.6 ± 0.5 m), and especially to the Fmax depth of green algae (3.4 ± 1.4 m).Consequently, the decrease of both SRP and SRP:TP at 3.8 m, combined with an increase in TP throughout the epilimnion, suggest uptake of SRP by phytoplankton and subsequent incorporation into biomass, resulting in the phytoplankton Taken together, these post-EM1 responses suggest that the increased total phytoplankton biomass in the treatment zone after the first mixing event may have been due to entrainment of suspended sediment from the upstream riverine zone of FCR into the deepest site, providing a nutrient subsidy that enabled an increase of green algae, cyanobacteria, cryptophyte, and overall total phytoplankton biomass.The possibility of a turbidity-induced nutrient subsidy is supported by three pieces of evidence: changes in nutrient concentrations post-EM1 at the depth of chlorophyll maximum in FCR, a back-of-the-envelope calculation of nutrient uptake rates at the depth of chlorophyll maximum post-EM1, and the difference in seasonal timing and intensity between EM1 and EM2, as described below.
First, in addition to an increase in total phosphorus across the treatment zone in FCR after EM1, there was also a decrease in SRP and the ratio of SRP:TP at 3.8 m.This depth (3.8 m) corresponds closely to the F max depth of phytoplankton biomass in the two weeks after EM1 (4.6 ± 0.5 m), and especially to the F max depth of green algae (3.4 ± 1.4 m).Consequently, the decrease of both SRP and SRP:TP at 3.8 m, combined with an increase in TP throughout the epilimnion, suggest uptake of SRP by phytoplankton and subsequent incorporation into biomass, resulting in the phytoplankton concentration increase.
Second, a back-of-the-envelope calculation of the uptake rate of SRP at 3.8 m during the week immediately after the first mixing event supports the possibility of a P subsidy leading to increased phytoplankton growth.From 30 May-6 June, SRP decreased by ~50%, from 12.3 to 6.6 µg L −1 at 3.8 m.Using the bathymetry of the reservoir to calculate SRP mass in the 2-m layer encompassing that depth and the microscopy-based total phytoplankton density of 1.7 × 10 8 cells L −1 at 3.8 m one week after mixing (6 June), we calculate an uptake rate of 5.9 × 10 −22 mol P cell −1 s −1 for the week, which is within the range of reported uptake rates of a diverse range of freshwater phytoplankton [47] (pp.151-161).This is a conservative estimate of uptake as phytoplankton biomass was lower in the week leading up to 6 June, which would imply a higher rate of P uptake per cell earlier in the week.However, even a much higher uptake rate would likely be less than the maximum P uptake rates of 10 −17 to 10 −14 mol P unit −1 s −1 exhibited by most phytoplankton [47] (pp.151-161).Furthermore, the ratio of the increase in total phytoplankton biomass to the increase in TP from the two weeks before the first mixing to the two weeks afterwards was 0.36, which corresponds closely with chlorophyll-a:TP ratios for water bodies similar to FCR [48].
Finally, because EM1 was conducted in late spring and at a high flow rate intensity of 25 SCFM, it is not surprising to see entrainment of turbid water after EM1 that was not replicated after EM2, which was conducted in the summertime and at a lower mixing intensity.Median upstream turbidity in FCR in the two weeks prior to EM1 (16 May-6 June) was 0.3 NTU higher than during the EM2 monitoring period of 13 June-27 June at both Site 1 and Site 2 (Figure S8), likely due to sediment in elevated spring flows, which either settled out or was redistributed and mineralized by the time of EM2.In addition, a more intense flow rate of 25 SCFM during EM1 may have caused more sediment stirring and quicker entrainment of suspended sediment than the lower flow rate of 15 SCFM during EM2.
The differences between EM1 and EM2 suggest that consideration of seasonal timing of mixing may be important in predicting ecosystem responses to epilimnetic mixing.The seasonal timing of our experiments was chosen based on when cyanobacterial blooms were historically highest in this reservoir (WVWA, unpublished data), and indeed, bloom-forming cyanobacteria were the dominant taxa at the F max depth prior to both mixing events (Figure 7, Table S1).However, the net increase in both cyanobacteria and total phytoplankton biomass after our first mixing experiment suggest that other environmental factors, such as the amount of recent precipitation in the watershed and the strength of thermal stratification, should also be considered when activating epilimnetic mixers.
Our second mixing experiment resulted in a change in community composition at the F max depth from a community dominated by filamentous cyanobacteria (specifically, Dolichospermum; MBFG 3) to a community of medium/large flagellates representing a diverse suite of phytoplankton taxa (MBFG 5).MBFG 5 is a spectrally diverse group, including green algae, brown algae, and cryptophytes (Table S1).This shift is in keeping with theoretical expectations from previous studies that mixing should disrupt large, colonial cyanobacteria while promoting taxa with rounder, denser morphologies.No change in community dominance occurred in our reference reservoir, which was dominated by medium to large flagellates throughout the EM2 monitoring period.
Taking a morphology-based approach, we found that while total phytoplankton biomass did not change after EM2, there was a substantial shift in the phytoplankton assemblage from large, high SA:V filamentous taxa to medium/large flagellated taxa at the depth of fluorescence maximum, and that large, high SA:V filamentous taxa consistently decreased while small, siliceous taxa consistently increased in the 24 h after both mixing events.This result has two implications for our understanding of phytoplankton responses to changes in mixing.First, we found it helpful to assess phytoplankton using morphological traits which are relevant to how phytoplankton experience mixing in the water column.Our use of morphological traits was focused on tracking changes in phytoplankton community dominance in response to mixing.However, it is possible that this approach could also be useful for managing cyanobacterial blooms; as cyanobacteria exhibit a wide variety of morphologies, using morphological traits to assess phytoplankton responses to management strategies could be an effective tool for water resource managers.An added benefit of the MBFG approach is that it is often easier to assess morphological traits than to conduct taxonomic identification of phytoplankton species, which requires extensive training and can be inconsistent among taxonomists.Second, if mixing is being used as a management strategy to decrease certain phytoplankton taxa but the assemblage is dominated by morphologies that tend to increase in biovolume or community dominance post-mixing, epilimnetic mixing will likely be ineffective in controlling phytoplankton biomass.
Our results also reiterate the importance of assessing phytoplankton assemblages using multiple measures across multiple sites and depths.While fluorescence-based sensors can give inaccurate biomass readings due to quenching or adaptive changes in pigment composition or pigment to biomass ratios of phytoplankton populations [49] and their aggregate data mask individual taxa responses, they are still useful for rapidly obtaining data on aggregated phytoplankton groups at finely resolved temporal and spatial scales.In our study, FluoroProbe biomass estimates were closely correlated with biomass estimates obtained using standard methods (absorbance of chlorophyll-a; Figure S3) and our functional trait data provided complementary information to the FluoroProbe for determining phytoplankton mixing responses.While our study only measured functional groups at the depth of the fluorescence maximum at the deepest site, we were able to measure aggregated phytoplankton groups across depth and along the reservoir gradient with the FluoroProbe (Figure 6).These upstream FluoroProbe casts provided valuable information showing increases in phytoplankton biomass across all sites after our first mixing event and not after our second mixing event, suggesting that phytoplankton responses along the reservoir gradient may have been similar to those observed at the deepest site.
We suggest that future whole-ecosystem mixing studies collect additional phytoplankton data at more resolved temporal and spatial scales as well as zooplankton data.Because phytoplankton can respond very rapidly (sub-hourly) to perturbation by adjusting their depth in the water column, we may have missed short-term responses to mixing overnight when sampling was not conducted.In future work, it would be helpful to take integrated depth samples of phytoplankton morphological groups at a higher time frequency before, during, and after epilimnetic mixing throughout the reservoir.Finally, while we measured a suite of physical and chemical variables that affect phytoplankton, we did not measure zooplankton, which can exert top-down grazing pressure.Moving forward, it would be useful to consider how zooplankton dynamics alter phytoplankton responses to epilimnetic mixing.
Consideration of both whole-ecosystem dynamics and functional traits with regards to phytoplankton mixing responses suggests possible future lines of inquiry.If epilimnetic mixing has the potential to entrain turbidity into new regions of a lake or reservoir, as demonstrated by Chen et al. (2018), which provides nutrient subsidies to phytoplankton, assessment of other, non-morphological functional traits, such nutrient uptake rates, could also be useful in predicting phytoplankton mixing responses.Furthermore, while we assessed the effects of short, relatively intense mixing pulses on phytoplankton in a small, shallow reservoir, epilimnetic mixers are operated over a variety of durations and intensities in lakes and reservoirs of varying sizes and bathymetries worldwide [3].Assessment of a range of mixing regimes in other lakes and reservoirs are needed to determine which ecosystem drivers and phytoplankton functional traits are most important for explaining phytoplankton responses to mixing and changes in mixed layer depth.

Conclusions
We found that epilimnetic mixing has variable effects on the concentration of total phytoplankton and the relative abundance of functional groups at the whole-ecosystem scale (Q1).Specifically, total phytoplankton biomass across a variety of morphology-based functional groups increased in the week

Figure 1 .
Figure 1.Bathymetric map of Falling Creek Reservoir (FCR), the focal experimental reservoir, and the reference Beaverdam Reservoir (BVR), Vinton, Virginia, USA.The inset is Falling Creek Reservoir, and EM denotes the epilimnetic mixer.

Figure 1 .
Figure 1.Bathymetric map of Falling Creek Reservoir (FCR), the focal experimental reservoir, and the reference Beaverdam Reservoir (BVR), Vinton, Virginia, USA.The inset is Falling Creek Reservoir, and EM denotes the epilimnetic mixer.

Water 2019 ,
11, x FOR PEER REVIEW 7 of 23

Figure 2 .
Figure 2. Water column temperature of (A) FCR and (B) BVR during summer 2016.The solid white line is the parent thermocline depth and the dashed white line is the lower metalimnetic boundary.Sampling dates are denoted by black points and the values in between these points are interpolated.The two mixing events in FCR (EM1 and EM2) are denoted by vertical black lines.

Figure 2 .
Figure 2. Water column temperature of (A) FCR and (B) BVR during summer 2016.The solid white line is the parent thermocline depth and the dashed white line is the lower metalimnetic boundary.Sampling dates are denoted by black points and the values in between these points are interpolated.The two mixing events in FCR (EM1 and EM2) are denoted by vertical black lines.

Figure 3 .
Figure 3. (A,C) Schmidt stability (left y-axis) and thermocline depth (right y-axis) in (A) Falling Creek Reservoir (FCR) and (C) Beaverdam Reservoir (BVR) during the two mixing experiments; (B,D) Turbidity in (B) FCR and (D) BVR during the mixing experiments.The two mixing events in FCR (EM1 and EM2) are denoted by vertical gray lines in (A and B).

Figure 3 .
Figure 3. (A,C) Schmidt stability (left y-axis) and thermocline depth (right y-axis) in (A) Falling Creek Reservoir (FCR) and (C) Beaverdam Reservoir (BVR) during the two mixing experiments; (B,D) Turbidity in (B) FCR and (D) BVR during the mixing experiments.The two mixing events in FCR (EM1 and EM2) are denoted by vertical gray lines in (A and B).

Figure 4 .
Figure 4. Epilimnetic concentrations of (A) green algae and (B) brown algae in Falling Creek Reservoir (FCR) before and after EM2.The response in depth profile patterns shown here exemplify patterns of recovery by phytoplankton spectral groups during both EM1 and EM2 (Figures S4 and S5).Pre-mix refers to 3-4 h prior to the experimental mixing, and post-mix refers to <1 h after mixing.

Figure 4 .
Figure 4. Epilimnetic concentrations of (A) green algae and (B) brown algae in Falling Creek Reservoir (FCR) before and after EM2.The response in depth profile patterns shown here exemplify patterns of recovery by phytoplankton spectral groups during both EM1 and EM2 (FiguresS4 and S5).Pre-mix refers to 3-4 h prior to the experimental mixing, and post-mix refers to <1 h after mixing.

Figure 5 .
Figure 5. (A,C,E,G,I) Mean treatment zone phytoplankton concentrations at the deepest site in both reference Beaverdam (Control, orange lines) and experimental Falling Creek (Impact, blue lines) Reservoirs during the two weeks before and after EM1 and (B,D,F,H,J) EM2.p-values are from the Before-After-Control-Impact (BACI) analyses (seeTable 1 for all statistics).Black points indicate linearly interpolated values.Mixing events are denoted by vertical black lines.Note y-axis scale differs by plot.

Figure 5 .
Figure 5. (A,C,E,G,I) Mean treatment zone phytoplankton concentrations at the deepest site in both reference Beaverdam (Control, orange lines) and experimental Falling Creek (Impact, blue lines) Reservoirs during the two weeks before and after EM1 and (B,D,F,H,J) EM2.p-values are from the Before-After-Control-Impact (BACI) analyses (seeTable 1 for all statistics).Black points indicate linearly interpolated values.Mixing events are denoted by vertical black lines.Note y-axis scale differs by plot.

Figure 6 .
Figure 6.(A-E) Longitudinal profiles of FCR (A) immediately before EM1, (B) immediately after EM1, (C) 24 h after EM1, (D) immediately before EM2, and (E) immediately after EM2.Vertical dashed lines represent where FluoroProbe profiles were taken.Location of Sites 1-5 can be found on the inset map of Figure 1.(F) Mean phytoplankton biomass across the top 5.5 m of the water column (or the full water column if shallower than 5.5 m) at all sampling sites in FCR during the monitoring period.Solid vertical lines denote mixing events.

Figure 6 .
Figure 6.(A-E) Longitudinal profiles of FCR (A) immediately before EM1, (B) immediately after EM1, (C) 24 h after EM1, (D) immediately before EM2, and (E) immediately after EM2.Vertical dashed lines represent where FluoroProbe profiles were taken.Location of Sites 1-5 can be found on the inset map of Figure 1.(F) Mean phytoplankton biomass across the top 5.5 m of the water column (or the full water column if shallower than 5.5 m) at all sampling sites in FCR during the monitoring period.Solid vertical lines denote mixing events.

Water 2019 , 23 Figure 7 .
Figure 7. (A) Biovolume of morphology-based functional groups (MBFGs) 1-7 and total biovolume of phytoplankton at the Fmax depth in FCR during the monitoring period.Solid vertical lines denote mixing events.(B) Proportional biovolume of the seven MBFGs in Falling Creek Reservoir at the depth of chlorophyll maximum, highlighting shifts in community dominance.The seven MBFGs are: 1-small, high SA:V cells; 2-small siliceous flagellates; 3-large, high SA:V filaments; 4-mediumsized, unspecialized cells; 5-medium/large flagellates; 6-non-flagellated siliceous cells; and 7large, mucilaginous, low SA:V colonies.Mixing events are denoted by vertical black lines.

Figure 7 .
Figure 7. (A) Biovolume of morphology-based functional groups (MBFGs) 1-7 and total biovolume of phytoplankton at the F max depth in FCR during the monitoring period.Solid vertical lines denote mixing events.(B) Proportional biovolume of the seven MBFGs in Falling Creek Reservoir at the depth of chlorophyll maximum, highlighting shifts in community dominance.The seven MBFGs are: 1-small, high SA:V cells; 2-small siliceous flagellates; 3-large, high SA:V filaments; 4-medium-sized, unspecialized cells; 5-medium/large flagellates; 6-non-flagellated siliceous cells; and 7-large, mucilaginous, low SA:V colonies.Mixing events are denoted by vertical black lines.

Figure 8 .
Figure 8. (A) Volume-weighted total phosphorus (TP) in the epilimnion of both Falling Creek Reservoir (FCR, blue) and Beaverdam Reservoir (BVR, orange) during the 2 weeks before and after EM1; (B) soluble reactive phosphorus (SRP) at the depth of fluorescence maximum (Fmax) in FCR in the 2 weeks before and after EM1; and (C) SRP:TP at the depth of Fmax in FCR in the 2 weeks before and after EM1.

Figure 8 .
Figure 8. (A) Volume-weighted total phosphorus (TP) in the epilimnion of both Falling Creek Reservoir (FCR, blue) and Beaverdam Reservoir (BVR, orange) during the 2 weeks before and after EM1; (B) soluble reactive phosphorus (SRP) at the depth of fluorescence maximum (F max ) in FCR in the 2 weeks before and after EM1; and (C) SRP:TP at the depth of F max in FCR in the 2 weeks before and after EM1.

Water 2019 , 23 Figure 9 .
Figure 9. (A-D) Turbidity across a longitudinal cross-section in Falling Creek Reservoir (FCR) from (A) immediately before EM1 (30 May), (B) 24 h after EM1 (31 May), (C) 4 days after EM1 (03 June), and (D) one week after EM1 (6 June).Contours are interpolated from CTD casts taken at the upstream and deepest sites (denoted by the vertical lines; upstream sampling sites are shown on the inset in Figure 1).

Figure 9 .
Figure 9. (A-D) Turbidity across a longitudinal cross-section in Falling Creek Reservoir (FCR) from (A) immediately before EM1 (30 May), (B) 24 h after EM1 (31 May), (C) 4 days after EM1 (03 June), and (D) one week after EM1 (6 June).Contours are interpolated from CTD casts taken at the upstream and deepest sites (denoted by the vertical lines; upstream sampling sites are shown on the inset in Figure 1).

Table 1 .
Results of Before-After-Control-Impact (BACI) analysis for physical, chemical, and biological parameters.Pre-EM column values display the mean ± 1 standard deviation in the two weeks before mixing, and post-EM column values display the mean ± 1 standard deviation in the two weeks after mixing.b Chemical variables are volume-weighted and summed over the treatment zone in FCR and equivalent zone in BVR.c Phytoplankton variables are averaged over the treatment zone in FCR and equivalent zone in BVR.d Significant values are highlighted in bold.e Asterisk denotes tests for which the sample size was too small to perform Welch's t-test.f F max depth = depth of fluorescence maximum. a