Mesozooplankton Community Dynamics and Grazing Potential Across Algal Bloom Cycles in a Subtropical Estuary

Mesozooplankton, as abundant grazers of microalgae in coastal systems, have the potential to prevent or mitigate harmful algal blooms (HABs) and their effects. The Indian River Lagoon (IRL) is a subtropical estuary in eastern Florida (United States) where repeated blooms, dominated by the toxic dinoflagellate Pyrodinium bahamense, the brown tide species Aureoumbra lagunensis, pico/nano planktonic cyanobacteria and other nano-eukaryotes, have highlighted the need to better understand fluctuations in the grazing potential of mesozooplankton populations across bloom cycles. Mesozooplankton and abiotic environmental data were collected at five sites in the northern IRL system at 6-week intervals from November 2013 through June 2016. A total of 107 taxa from 14 phyla were detected. Communities varied across sites, dates and between bloom and non-bloom periods, with densities up to 338 individuals L–1. Eight taxa comprising 85–94% of the total population at each site were identified as primary potential grazers, including barnacle nauplii, cladocerans, adult copepods, gastropod veligers, larvaceans, and polychaete metatrochophores. Although abundant, the estimated grazing potential of the primary taxa, calculated from their measured densities and previously published grazing rates, suggest that mesozooplankton lack the capacity to suppress phytoplankton once they reach bloom levels. These findings illustrate the utility of monitoring data and underscore the importance of systematically evaluating algal bloom controls with a consideration for the dynamic conditions of each unique ecosystem.


INTRODUCTION
Mesozooplankton are important grazers of phytoplankton in a wide variety of marine systems. Comprised of holoplankton such as copepods and meroplanktonic larvae ranging from 200 µm to 2 mm, mesozooplanktonic grazers ingest an estimated 12-18% of global oceanic primary production annually (Calbet, 2001;Calbet and Saiz, 2005). Ingestion rates are typically greatest in highly productive coastal systems (Calbet, 2001) such as estuaries, where eutrophication and other anthropogenic alterations support the formation of phytoplankton blooms (Sunda et al., 2006;Anderson et al., 2008;Heisler et al., 2008;Raven et al., 2020).
The Indian River Lagoon (IRL) estuary is a biodiverse barrier island lagoon stretching from warm temperate to subtropical climate zones across 251 km of Florida's east coast (United States). The estuary system consists of three connected regions: the Mosquito Lagoon, the Indian River Lagoon and the Banana River Lagoon. Phytoplankton blooms have long been observed in the IRL, but they experienced a dramatic increase in intensity and change in composition in 2011 (Phlips et al., 2021). Prior to 2011, relatively large dinoflagellates and diatoms were the dominant bloom species (Phlips et al., 2015). However, widespread seagrass losses in 2009-2010, and extraordinarily cold winter water temperatures in 2010, were followed by blooms of small-celled species (<5 µm in diameter). In 2011, picocyanobacteria and Pedinophyceae species bloomed throughout the northern IRL, followed by repeated intense blooms of the brown tide pelagophyte Aureoumbra lagunensis, and other small-celled taxa (Phlips et al., 2021). None of these taxa were documented at such levels in the estuary prior to 2011 but have occurred repeatedly in the years since (Phlips et al., 2015(Phlips et al., , 2020Schaefer et al., 2019). The blooms have decreased light penetration and dissolved oxygen concentrations, leading to large-scale losses of habitat and changes in biological communities (Phlips et al., 2015;Lapointe et al., 2020;Lewis et al., 2020Lewis et al., , 2021Lazensky et al., 2021) that have defined the IRL for decades.
Harmful algal blooms (HABs) are not isolated to the IRL or similar systems. In the United States alone, such blooms occur in all 50 states, caused by over 100 taxa (Anderson et al., 2021) that impact human and ecosystem health, tourism, and commercial fishing and aquaculture industries with a combined value of $7 billion in 2018/2019 (National Marine Fisheries Service, 2021). Climate change and other anthropogenicallydriven environmental effects are expanding ranges and increasing bloom frequencies of many phytoplankton taxa (Gobler, 2020). It is therefore crucial to better understand bloom impacts on ecosystems, and to evaluate potential top-down controls such as grazing by a variety of organisms.
Microzooplankton (20-200 µm) are recognized for their potential to graze blooms of small-celled phytoplankton (e.g., Sautour et al., 2000;Sherr and Sherr, 2007;Calbet, 2008;Schmoker et al., 2013) like those occurring in the IRL system. In contrast, mesozooplankton have received less attention as top-down bloom controls because they are typically deemed inefficient consumers of small cells (Calbet and Landry, 1999;Calbet et al., 2003;Davis et al., 2012;Morison et al., 2020). However, some studies have documented the ability of mesozooplankton to consume particles <5 µm (Calbet et al., 2000;Davis et al., 2012). Furthermore, these larger zooplankton are major components of the IRL system (Youngbluth et al., 1976;Rey et al., 1991;Dix and Hanisak, 2015), are a critical link between phytoplankton and higher trophic levels (Sommer and Stibor, 2002;Lobry et al., 2008;David et al., 2016), and should be considered in the big picture as bloom impacts and potential control mechanisms are explored. This study investigated mesozooplankton in the IRL system with the following objectives: (1) Evaluate the spatial and temporal variability of mesozooplankton communities at bloom hotspots. (2) Estimate the grazing potential of mesozooplankton populations on natural levels of phytoplankton biomass. (3) Identify key mesozooplankton grazers as potential candidates for bloom control.

Field Procedures
Mesozooplankton were sampled at five sites across the northern Indian River Lagoon (IRL) system (Figure 1) that were selected because of previous phytoplankton bloom activity (Phlips et al., 2015):  Table 1). Daytime plankton samples (n = 4) were collected with a ring net (153-µm mesh, 30-cm diameter) every 6 weeks from November 17, 2013 through June 27, 2016 for all sites except COA, which began on June 25, 2014. The net was towed horizontally within 1 m of the surface, effectively representing the entire shallow and generally well-mixed water column, which was supported by the presence of several benthic and demersal taxa (e.g., Caprella and Cerapus amphipods, Sphaeroma isopods, Oxyurostylis cumaceans, and several species of benthic foraminiferans and harpacticoid copepods) in the samples. Tow duration was 1 min and was reduced as needed to prevent clogging, which occurred when phytoplankton and zooplankton densities were at their highest and when byproducts (e.g., mucus and molts) were abundant. Tow volumes (0.02-35.63 m 3 , mean = 7.98 ± 0.41 m 3 ) were calculated from tow speeds and confirmed by in-net (Model 2030R, General Oceanics, Inc., Miami, FL, United States) and external (Model FP211, Global Water Instrumentation, Inc., Gold River, CA, United States) flowmeters. All resulting samples were immediately preserved in 5% formalin with sodium tetraborate buffer.
During each sampling event, total depth (m), water temperature ( • C) and salinity (ppt) were measured with a CastAway R CTD (SonTek, San Diego, CA, United States). Average water column temperatures and salinities were calculated using data from the profiles at standardized depths (0.15, 0.45, 0.75, and 1.05 m) within the range of a typical net tow. Surface dissolved oxygen (DO, mg L −1 ) and pH were measured with a YSI Model 85 handheld multiparameter meter (YSI, Inc., Yellow Springs, OH, United States). Water clarity was determined via Secchi disk depth (m) using a weighted, limnological disk (20-cm diameter).

Laboratory Procedures
After 1-2 weeks of fixation, samples were re-sieved (153 µm), transferred to 1-µm filtered seawater, and divided zero to eight times with a Motoda plankton splitter (Motoda, 1959) to obtain aliquot densities at which all individuals could be easily differentiated in a 10 cm × 10 cm gridded tray. Using stereomicroscopy (35× magnification), zooplankton from each grid were enumerated and identified to the lowest possible taxonomic level using several sources (e.g., Davis, 1949;Grice, 1960;Smith and Johnson, 1996;Johnson and Allen, 2012) until ≥300 individuals were recorded (Davidson and Clem, 2003;BS EN 14407:2004. Counts were normalized by volume for comparison across samples. Means are ± 1 SE.

Phytoplankton and Grazing Potential
Algal biomass data (A. lagunensis, cyanobacteria, diatoms, dinoflagellates, others) for the study period were obtained from weekly to biweekly monitoring by the Phlips Lab at the University of Florida. At each site, daytime phytoplankton samples (n = 5) were collected with an integrated water sampler (Sch 40 PVC pipe, 3.2 cm diameter) lowered vertically from the surface to within 0.1 m of the bottom. The five samples were pooled to reduce the effect of patchiness, and duplicate aliquots were preserved, one with Lugol's iodine solution and one with buffered glutaraldehyde (Phlips et al., 2010(Phlips et al., , 2020. Phytoplankton composition was determined using the Utermöhl method (Utermöhl, 1958). Preserved samples were settled in 19-mm diameter cylindrical chambers. Phytoplankton cells were identified and counted at 400 and 100× magnification using inverted phase contrast microscopy (Leica Microsystems, Wetzlar, Germany). At 400×, 30-100 grids were examined until a minimum of 100 cells of a single taxon were counted. At 100×, a total bottom count was obtained for taxa >30 µm in size. Fluorescence microscopy was used on glutaraldehyde-preserved aliquots to enumerate picoplanktonic cyanobacteria (i.e., picocyanobacteria, which is mainly spherical Synechococcus spp. in the IRL) at 1,000× magnification (Phlips et al., 1999). Light microscopy was used to identify and enumerate A. lagunensis based on morphological features evident from previous work (Phlips and Badylak, 2013), which initially characterized A. lagunensis via light, scanning electron and transmission electron microscopy that was later supported by immunological assays and 18S rRNA gene sequencing conducted by colleagues at the Florida Fish and Wildlife Research Institute and Stony Brook University. To estimate cell biovolumes, subsamples of seawater were filtered onto 0.2 µm Whatman R Nuclepore TM filters (GE Healthcare, Chicago, IL, United States) and mounted with immersion oil between a microscope slide and cover slip. Biovolumes were estimated by assigning combinations of geometric shapes to fit the cell characteristics of individual taxa (Smayda, 1978). Specific phytoplankton dimensions were measured for ≥100 randomly selected cells. Species of variable size (e.g., many diatom taxa) were placed into size categories. Phytoplankton carbon values (µg C mL −1 ) were estimated by multiplying biovolume estimates (expressed as 10 6 µm 3 mL −1 ) by conversion factors for different taxonomic groups: 0.061 for diatoms, 0.22 for cyanobacteria and small picoplanktonic eukaryotes (including A. lagunensis), and 0.16 for dinoflagellates and other taxa (Strathmann, 1967;Ahlgren, 1983;Sicko-Goad et al., 1984;Work et al., 2005).
The phytoplankton data were used to characterize sites, define bloom periods and compare with estimated grazing rates. Bloom periods at each site were defined as dates with an algal biomass >2 µg C mL −1 (Phlips et al., 2015(Phlips et al., , 2020 for at least one of the five phytoplankton categories. The phytoplankton monitoring schedule was more frequent and sometimes offset from that of the mesozooplankton monitoring. Therefore, a zooplankton sampling event was considered to occur during a phytoplankton bloom if bloom conditions were observed the same day that the zooplankton were collected, or both before and after (typically 1-7 days) they were collected. The same strategy was employed to define non-bloom periods. Relative chlorophyll (RFU) data from nearby continuous water quality monitoring stations (St. Johns River Water Management District Aquarius WebPortal 1 ) were used to support these designations. If the bloom status on a particular zooplankton sampling date remained unresolved, that date was defined as 'unknown' and was excluded from bloom/non-bloom comparisons.
Key potential grazers were identified as primarily herbivorous or omnivorous taxa reaching densities >10 ind. L −1 during the study period. Data on these taxa were used to demonstrate whether mesozooplankton may have the capacity to suppress algal blooms in this estuary through grazing. Ingestion rates (IR) were obtained from studies of the key grazers, or closely related taxa known to occur in the region, at the same or similar life stages, and examined at water temperatures within the range of this study (10.6-31.8 • C). The minimum and maximum estimated grazing potentials (EGP min and EGP max , µg C mL −1 d −1 ) for each taxon during each sampling event were obtained by multiplying their in situ densities (ind. mL −1 ) by their IR min and IR max values (µg C ind. −1 d −1 ), respectively. The EGP min and EGP max values were summed across the suite of eight key grazers to calculate the minimum and maximum Estimated Metagrazing Potential (EMP min and EMP max , µg C mL −1 d −1 ) for the entire mesozooplankton community at each site over time. This EMP range provides an extremely rough but relatively inclusive estimate of mesozooplankton grazing potential at a site -'best' and 'worst-case' scenarios that consider a variety of conditions and prey types.

Data Analysis
Univariate statistical analyses were conducted in The Real Statistics add-in package for Microsoft Excel 2 . Following Levene's tests for homoscedasticity, Welch's analysis of variance (ANOVA, α = 0.05) and Games-Howell post hoc tests were used to detect differences among sites for abiotic parameters, species richness (S) and zooplankton density (ind. L −1 ). Welch's t-tests were used to identify differences in key grazer densities from bloom and non-bloom periods within each site.
Multivariate community analyses were performed on square root transformed count data using the PRIMER v5 software package (PRIMER-E Ltd., Plymouth, United Kingdom). A twoway crossed analysis of similarities (ANOSIM, α = 0.05, Bray-Curtis similarities) was used to find differences in zooplankton community composition among sites and dates. Within each site, differences in communities between bloom and non-bloom periods were detected with a one-way ANOSIM. Groups were considered different only when R > 0.250 and p < 0.05 (Clarke and Gorley, 2001). Taxa driving the differences were identified with similarity percentage (SIMPER) analyses.

RESULTS
Across all sites, a total of 439 samples were collected and 107 taxa from 14 phyla were identified. Per sample mesozooplankton richness and density ranged from 5 to 25 taxa and 0.1 to 338.1 ind. L −1 , respectively. 2 real-statistics.com

Site Characteristics
Key physical and biological characteristics varied among the five sites during the study period (Table 1). Pairwise comparisons revealed that MLA was generally the most saline site (compared with TIV, BRC, and COA: p ≤ 0.026) with the lowest water clarity (as Secchi depth, compared with TIV and COA: p ≤ 0.002). No differences were detected among sites with regard to water column temperature, surface DO or surface pH (p ≥ 0.313). MLA and TIV exhibited higher species richness than BRN and COA (p ≤ 0.038), with BRC situated between the two groups. In contrast, MLA and TIV, along with BRN, hosted lower zooplankton densities than one or both of the other sites (p ≤ 0.012).
Zooplankton communities were dominated by arthropods (primarily copepods) at all sites, with MLA also containing a large percentage of molluscs ( Table 1). A. lagunensis was the prevailing alga overall (43.2-56.1%), with percentages of other groups varying considerably by site. Bloom conditions occurred across all sites at a wide range of water temperatures during both wet and dry seasons (Figure 2), with blooms being detected from 5 to 26 days at BRN and MLA, respectively. Prior to June 2015, phytoplankton communities at all sites were assemblages of cyanobacteria, diatoms, dinoflagellates and others, which formed occasional blooms up to 12 µg C L −1 . During this period, A. lagunensis was a minor contributor to the system (Figure 2). In late 2015, total algal biomass dominated by A. lagunensis and 'other' taxa, rose at all sites and remained elevated for several months, peaking at 29.89 µg C L −1 (85% A. lagunensis) at BRN on March 07, 2016. At the height of this bloom, the presence of other phytoplankton taxa was minimal.

Grazing Potential
Eight taxa comprising 85-94% of the total population at each site were identified as primary potential grazers: Amphibalanus spp. barnacle nauplii, Evadne sp. cladocerans, gastropod veligers (one species), Oikopleura sp. larvaceans, polychaete metatrochophores (one species), and the copepods A. tonsa, O. colcarva, and P. crassirostris. Of 40 density comparisons of these taxa between bloom and non-bloom periods, 23 were lower during bloom events (p ≤ 0.049, Table 2). The greatest impacts were seen on FIGURE 2 | For each of the five mesozooplankton sampling sites, relative percentages (top) and biomasses (µg C mL −1 , bottom) of all five phytoplankton categories monitored during the study. The gray dashed lines denote the 2 µg C mL −1 bloom threshold.
populations at BRN and BRC, where all eight and seven out of eight key grazers declined during blooms, respectively.
A review of the relevant literature on these taxa and their close relatives uncovered per-species ingestion rates (IRs) ranging from 0.001 to 13.1 µg (1-13,100 ng) C ind .−1 d −1 for Oithona and Oikopleura, respectively ( Table 3). The largest community-wide Estimated Metagrazing Potential (EMP max ) of 1.22 µg C mL −1 d −1 was documented at TIV on October 1, 2015, resulting from high densities of all three copepod species. On several occasions, the EMP max exceeded the minimum phytoplankton biomass documented at a particular site (Figure 3), signaling that mesozooplankton grazing can be high enough to impact phytoplankton populations. However, EMP max never exceeded the 2 µg C mL −1 bloom threshold, suggesting that mesozooplankton lack the grazing capacity needed to suppress phytoplankton in this system once they reach bloom levels.

DISCUSSION
This investigation documents the spatial and temporal variation in mesozooplankton communities from a subtropical estuary affected by recurring phytoplankton blooms and provides a rough estimate of their grazing potential. It is important to acknowledge some limitations of this study. First, although the zooplankton tows sampled the majority of the well-mixed water column, only collecting during the day could have poorly estimated the abundances of some grazers due to diel vertical migration. Secondly, although grazing potential ranges were based on rates measured in a variety of conditions, they were still calculated indirectly from other studies. Grazing within a system is affected by a combination of factors, including those not examined by this study, such as phytoplankton quality and the abundance of alternative microzooplankton prey for omnivorous grazers. However, this study presents a method of utilizing datasets often readily available to scientists and resource managers working in bloom-impacted estuaries, providing an approximation of grazing potential to further monitoring and management efforts when direct grazing measurements are unavailable. Results suggest that, while mesozooplankton grazing in the Indian River Lagoon system may exert some control on phytoplankton when abundances are low, it was not sufficient during the study period to suppress phytoplankton populations once they reached bloom levels. The largest Estimated Metagrazing Potential (EMP max ) calculated for the mesozooplankton community was 1.22 µg C mL −1 d −1 , lower than the 2 µg C mL −1 threshold used herein to define bloom events. The brown tide pelagophyte A. lagunensis exhibited the largest algal biomass of all the phytoplankton groups surveyed during the study period. The prevalence of this alga and the composition of the zooplankton communities may help explain the apparent inability of mesozooplankton grazing to control blooms in the northern IRL system. Two major factors affecting zooplankton grazing efficacy are prey cell size and the trophic modes of potential grazers. Many of the most abundant bloom-forming phytoplankton in the IRL and elsewhere are <5 µm, including A. lagunensis, pedinophytes, Synechococcus and other picocyanobacteria. Several studies have documented a tendency for mesozooplankton to graze on larger cells. Such findings are consistent with this study, which recorded the peak EMPs prior to the height of the A. lagunensis bloom, when algal biomasses were dominated by larger diatoms and dinoflagellates.
In addition, omnivorous mesozooplankton may also prey on microzooplankton that are often the more effective grazers of small autotrophs (Sautour et al., 2000;Calbet and Landry, 2004;Liu et al., 2005;Calbet, 2008;Davis et al., 2012). In this study, the copepods A. tonsa and O. colcarva were the most abundant mesozooplankton. Both genera are omnivorous and consume microzooplankton as a major and sometimes dominant part of their diets (Lonsdale et al., 1979(Lonsdale et al., , 2000Castellani et al., 2005). Numerous studies have documented feeding in Acartia and Oithona under a variety of conditions. Jonsson and Tiselius (1990) observed the feeding behavior of A. tonsa on several ciliate taxa and on the alga Cryptomonas baltica. The authors noted that the copepods would devote time to suspension feeding on algal cells only when they were particularly abundant. At lower algal densities, A. tonsa would switch to raptorial feeding on ciliates. Rollwagen-Bollens and Penry (2003) investigated the feeding dynamics of Acartia spp. at two locations in San Francisco 3 | Ingestion rates (IR, µg C ind −1 d −1 ) from studies of key grazers, or closely related taxa known to occur in the region, at the same or similar life stages, and examined at water temperatures within the range of this study.  *Unspecified water temperature within 1-2 • C of ambient in ship channel, Port Aransas, TX, United States; † senescent culture; ‡ growing culture; § mixed culture with ciliates; } stage II nauplii; }} stage III nauplii; }}} stage VI nauplii; Ⴋ autotrophic nanoflagellates < 10 µm; í low silica; ɦ high silica.

IR
Bay (CA, United States) that hosted different assemblages of potential prey. The findings revealed that the copepods consumed different prey at the two sites but tended to target motile ciliates and flagellates >10 µm. Using in situ mesocosms containing natural densities of A. lagunensis and Synechococcus sp. from the upper Laguna Madre lagoon (Texas, United States), Buskey et al. (2003) tested whether A. tonsa indirectly affected phytoplankton densities by consuming herbivorous ciliates. Ciliate numbers increased in mesocosms where the copepods were removed, signaling that A. tonsa preyed on the microzooplankton.
However, the rise in ciliate grazers had different impacts on the two autotrophs, reducing Synechococcus sp. but leaving A. lagunensis populations unaffected. A laboratory study by Saiz et al. (2014) on Oithona davisae feeding revealed that adult females preyed most heavily on organisms >4 µm, including ciliates, heterotrophic flagellates and nauplii from other copepod species. The authors hypothesized that these preferences were the result of ambush feeding behavior that requires prey to be large and motile enough to trigger hydromechanical detection. Finally, in a study of the variability of mesozooplankton grazing FIGURE 3 | Estimated Metagrazing Potential maxima (EMP max , dark orange lines) and minima (EMP min , pale orange lines) for the suite of key mesozooplankton grazers across all five monitoring sites. To aid comparisons between grazing potential and phytoplankton biomass, the phytoplankton bloom threshold (2 µg C mL −1 , gray dashed lines) and the minimum phytoplankton biomass recorded at each site (blue solid lines), are denoted.
on phytoplankton at two contrasting sites in Hong Kong, Chen et al. (2017) documented a high degree of carnivory in mesozooplankton populations comprised largely of Acartia and Oithona. The preference for microzooplankton prey reduced mean grazing impacts on phytoplankton to <4% at both sites.
Another factor likely contributing to the persistence of A. lagunensis in the northern IRL system is a set of characteristics and behaviors exhibited by the species that reduces grazer fitness: (1) potential low nutritive value and toxicity, (2) allelopathy toward competitive phytoplankton, and (3) the secretion of protective exopolymers. During the brown tide-rich blooms in late 2015 and 2016, mesozooplankton populations declined at all sites, presumably driven by deteriorating conditions from the rise in A. lagunensis biomass. After zooplankton declined in the Laguna Madre lagoon (Texas, United States) following an outbreak of A. lagunensis, Buskey and Hyatt (1995) investigated the condition of A. tonsa fed on cultures of the brown tide alga. Adult egg production rates and naupliar survival both decreased. The findings suggested that A. lagunensis was mortally toxic to nauplii and a poor food for adults, producing egg release rates similar to those of starved copepods. In a study that assessed the development of A. tonsa nauplii grazing on another brown tide organism, Aureococcus anophagefferens, Smith et al. (2008) found that naupliar development was delayed when larvae consumed the alga. Interestingly, the nauplii selectively fed on more palatable algae only when the cells were larger than A. anophagefferens.
Recent work has revealed that brown tides can also inhibit the growth of competing phytoplankton, potentially lowering the availability of other food sources for grazers. In an investigation of the allelopathic effects of A. lagunensis and A. anophagefferens on a wide variety of phytoplankton, Kang and Gobler (2018) documented reductions of up to 96% in the abundance of competing autotrophs. Filtrate of Aureoumbra cultures also reduced the photosynthetic efficiency of Synechococcus and the diatom Chaetoceros calcitrans by 53 and 19%, respectively. A particularly concerning discovery was that strains of Aureoumbra isolated from the IRL in 2012 and 2013 had a higher allelopathic potency compared to those obtained from Laguna Madre, Texas 20 years ago. The 2015-2016 blooms in the present study were usually dominated by A. lagunensis biomass, but frequently contained elevated levels of other phytoplankton as well. However, diatoms and dinoflagellates, two well-documented primary food sources for mesozooplankton, never exceeded 0.64 µg C mL −1 during this period. These results may reflect the use of allelopathic chemicals by A. lagunensis to reduce competitive phytoplankton populations.
The utilization of extracellular polymeric substances (EPS) may also give A. lagunensis a competitive edge. In laboratory experiments using cultures isolated from Texas, Liu and Buskey (2000a) found that the alga formed a mucus layer of secreted EPS. As populations progressed from growth to decline, per cell EPS production increased threefold, which was maximized under hypersaline conditions. The authors continued their investigation by feeding high and low EPS cultures of A. lagunensis to three protozoan species common in the bloom area (Liu and Buskey, 2000b). All taxa showed reduced growth rates in the high EPS treatment. A subsequent experiment with one of those taxa, the ciliate Aspidisca sp., revealed that the slower growth resulted from reduced grazing, possibly from mucus clogging its feeding apparatus or from impaired grazing caused by the altered swimming behavior observed in the high EPS treatment.
Two of the most promising candidates for grazing down brown tide and other phytoplankton during this study were the copepod P. crassirostris and the larvacean Oikopleura sp. Both taxa were major contributors to the EMP peaks documented at COA and BRC, and while A. tonsa exhibited the highest densities overall, P. crassirostris was the most abundant primarily herbivorous grazer. However, densities of the copepod dropped by an average 81% during bloom periods. Researchers have also documented reduced feeding and increased mortality of P. crassirostris during grazing experiments on A. lagunensis and two other common IRL bloom taxa (Ma et al., 2021a,b and unpublished data), suggesting it is not a viable mitigator of many algal blooms in this estuary. In contrast, Oikopleura sp. was significantly more abundant during bloom periods at COA, although this was driven largely by its high densities during diatom and dinoflagellate-rich blooms in 2014 and summer 2015. However, the larvacean was still one of the most abundant taxa, with densities >3 ind. L −1 , during the start of the 2015-2016 A. lagunensis blooms at MLA and TIV. Although we are unaware of any studies investigating the effects of A. lagunensis on Oikopleura, Badylak and Phlips (2008) tracked the abundances of Oikopleura dioica and the copepods A. tonsa and O. colcarva during a bloom of the toxic dinoflagellate P. bahamense in Tampa Bay (FL). While populations of both copepods declined, O. dioica abundances mimicked those of P. bahamense, suggesting it was more tolerant to the toxic alga. This resilience during some blooms, coupled with the ability to ingest copious small phytoplankton caught in mucous feeding nets (Acuña and Kiefer, 2000;Scheinberg et al., 2005;Troedsson et al., 2007), positions Oikopleura as the most likely mesozooplankter to suppress blooms of phytoplankton <5 µm in the northern IRL system.
However, another layer of complexity exists for grazers able to consume A. lagunensis, in particular. The same EPS coating that deters some herbivores may keep the algal cells from being digested or even damaged when they are consumed. In a laboratory study testing the viability of A. lagunensis consumed by A. tonsa, Bersano et al. (2002) found undamaged cells that proceeded to inoculate new algal cultures after being sequestered inside copepod fecal pellets for 3 days. The authors postulated that the EPS layer helped protect cells from digestion, allowing them to re-enter the water column and persist as primary producers. This process has been documented for other algae as well. Porter (1976) detailed the grazing of the green alga Sphaerocystis schroeteri by freshwater cladocerans from the genus Daphnia. Cells of this autotroph form colonies that are surrounded by a protective gelatinous sheath. Observations revealed that S. schroeteri was fed on at a lower rate than unsheathed algal taxa and was resistant to digestion when consumed. In fact, >90% of ingested cells passed through the gut undamaged. The author noted that grazing did break up cell clusters and remove the protective sheath of many colonies, but this removal allowed for the uptake of nutrients by S. schroeteri from within the guts of Daphnia, boosting algal growth that could compensate for the minor losses due to grazing.
These studies provide ample evidence that A. lagunensis can reduce grazing opportunities and foster malnourishment, starvation and population declines of mesozooplankton grazers. In the present study, dramatic reductions of primary grazers at all sites during the most intense blooms may have been caused by one or more of these 'anti-grazer' mechanisms, allowing the blooms to persist for several months. The IRL system is a complex and diverse estuary where phytoplankton communities that were dominated by diatoms and dinoflagellates have been largely displaced over the past decade by smaller taxa that appear to evade grazing by mesozooplankton. The frequency of smallcelled phytoplankton blooms is likely to increase as warming ocean temperatures favor these taxa over larger autotrophs (Flombaum et al., 2013;Phlips et al., 2015;Cloern, 2018). Recent research has revealed additional competitive advantages for some small species, such as the discovery of a protective resting stage for A. lagunensis that may facilitate future blooms and geographic expansion (Kang et al., 2017). This projected spread forewarns that more coastal systems will grapple with similar threats to economic, public and ecosystem health stemming from blooms of shared taxa. However, each estuary facing recurring phytoplankton blooms is unique with changing biological communities that provide different probabilities for top-down control. Therefore, it is important to systematically evaluate these potential controls with a consideration for the dynamic conditions of each ecosystem.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KJ, LS, and HA conceived and designed the mesozooplankton monitoring plan. LS and HA collected and curated the mesozooplankton and abiotic environmental data. EP supplied the phytoplankton data. LS analyzed the data and wrote the manuscript with support from KJ, HA, and EP. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by St. Johns River Water Management District, award #27786. Funding for the publication of this manuscript was provided by the Smithsonian National Museum of Natural History and by the Indian River Lagoon National Estuary Program (IRLNEP) through an EPA Region 4 supplemental funding grant. This is Smithsonian Marine Station contribution #1168.