Temporal and spatial variability of phytoplankton and mixotrophs in a temperate estuary

: A significant proportion of phototrophic species are known to be mixotrophs: cells that obtain nutrients through a combination of photosynthesis and prey ingestion. Current methods to estimate mixotroph abundance in situ are known to be limited in their ability to help identify conditions that favor mixotrophs over strict autotrophs. For the first time, we combine microscopic analysis of phototrophic taxa with immunoprecipitated bromodeoxyuridine (BrdU)-labeled DNA amplicon sequencing to identify and quantify active and putative mixotrophs at 2 locations in a microtidal temperate estuary. We analyze these data to examine spatial and temporal variability of phytoplankton and mixotrophs. Microscopy-based phototrophic diversity and abundances reveal expected seasonal patterns for our 2 stations, with the start of growth in winter and highest abundances in summer. Diatoms tend to dominate at the site with less stratification, while dinoflagellates and euglenids are usually more prominent at the stratified station. The BrdU-based mixotroph identifications are translated to the microscopy identification and abundances to estimate the proportion of mixotrophs (cells >10 μm in size) at both sites. The average proportion of potential mixotrophs is higher at the station with higher stratification (51%) compared to the station with lower stratification (30%), and potential mixotrophs tend to be higher in summer, although we did not conduct BrdU experiments in any of the other seasons. Combining the identification of active mixotrophs through the uptake of BrdU-labeled bacteria with robust abundance measurements can expand our understanding of mixotrophs across systems.


INTRODUCTION
Phototrophic plankton form the base of most aquatic food webs, producing organic carbon that will be utilized by higher trophic levels. Phototrophic biomass is commonly estimated by measuring chl a concentrations in water samples because the chl a pigment is present in almost all phototrophs and has been shown to be a direct proxy for their biomass (Steele 1962, Strickland 1965, Cullen 1982, albeit an imperfect one (Buchanan 2020). As a result of the emphasis on the collection of chl a concentration data, the potential importance of mixotrophs in per-forming some combination of both photosynthesis and ingestion goes largely unknown (Millette et al. 2018). It is necessary to understand the nutritional modes being used by the phototrophic assemblage because these factors contribute to variability in biogeochemical cycling (Lomas et al. 2014, Brierley 2017, water quality (Shashi Shekhar et al. 2008, Valiela et al. 2016, and productivity and stability (Ptacnik et al. 2008, Vallina et al. 2014, Tian et al. 2017. For example, recent modeling studies that examined the potential contribution of mixotrophy to carbon cycling propose that mixotrophs could substantially increase the amount of primary production and the transfer of carbon to higher trophic levels (Mitra et al. 2016, Ward & Follows 2016. A major impediment to incorporating mixotrophy into environmental studies is the inability to reasonably estimate the abundance of mixotrophs in a system. Fluorescently labeled bacteria and fluorescent microspheres have long been used to estimate mixotroph abundance and ingestion rates, but research shows that these methods underestimate true abundance (Anderson et al. 2017, Li et al. 2021. More recently, a few studies have identified potential mixotrophs in field samples via taxonomy, referring to phototrophs that have been found to consume prey in previous experiments (Haraguchi et al. 2018, Leles et al. 2019, Cesar-Ribeiro et al. 2020, Schneider et al. 2020. For example, Schneider et al. (2020) used microscopy-based taxonomic information from a long-term North Sea dataset to link the occurrence and abundance of potential mixotrophs with abiotic conditions (offshore versus estuarine). While this method of potential mixotroph identification provides the opportunity to reanalyze available microscopy samples and estimate the proportion of mixotrophs, it likely overestimates their abundance since species presence does not necessarily translate to mixotroph behavior.
Here, we apply a relatively new molecular method to constrain the identification of mixotrophs in microscopy data to known active mixotrophs in the same environment. By feeding live bacteria labeled with bromodeoxyuridine (BrdU) to a plankton sample, it is possible to ascertain cells that ingest the bacteria and incorporate the BrdU into their DNA through an immunoprecipitation process (Fay et al. 2013). From there, the immunoprecipitated samples are amplified and sequenced to identify the mixotrophs taxonomically as amplicon sequence variants (ASVs) that are associated with protists known to contain chloroplasts (Gast et al. 2018). Mixotrophs identified from the BrdU experiments can then be matched with taxa identified in microscopy samples to estimate the abundance and proportion of mixotrophs greater than 10 μm in size. The BrdU method improves upon previous research by curating a list of active mixotrophs specific to our study site. Improved estimation of the abundance and proportion of mixotrophs can advance our understanding of how mixotrophs vary in response to local environmental conditions.
The goal of this study is to investigate the spatial and temporal variability of the phytoplankton (strict phototrophs) and mixotrophs present in a temperate estuary (Waquoit Bay, MA, USA) and identify possible environmental conditions driving their variabil-ity. Due to intensive long-term monitoring, there is good temporal and spatial coverage on chl a concentration and other environmental data in Waquoit Bay, but data on the phototrophic taxa present are extremely limited. This research project has 3 specific objectives: (1) identify the phototrophic taxa present in Waquoit Bay, (2) determine which of these taxa are capable of phagotrophy, and (3) examine what environmental factors are associated with the spatial and temporal variability in phyto-and mixotrophic abundance and proportion.
Data were collected on the abundance of phototrophic taxa and environmental conditions over a 15 mo period (2018−2019) at 2 of the long-term monitoring locations in the Waquoit Bay Estuary. Laboratory experiments were conducted over a 3 mo period in the summer of 2019 to determine which phototrophic taxa were capable of ingesting bacteria. This study provides some of the first data on the proportion of potential mixotrophs within the phototrophic assemblage throughout the year in a temperate estuary.

Waquoit Bay
The Waquoit Bay Estuary (Fig. 1) is a shallow microtidal system that is part of the National Estuarine Research Reserve System (NERRS). At its deepest, the estuary is approximately 3 m deep at low tide, but on average the estuary is only 1.5 m deep at mean low tide. The average tidal range for the estuary is roughly 1 m. Salinity ranges from 27 to 32 at the mouth of the estuary, while the tidally influenced tributaries experience a much greater salinity amplitude of 5 to 30.
Each reserve within the NERRS is required to maintain at least 4 System-Wide Monitoring Program (SWMP) stations where a series of water quality data are collected (Centralized Data Management Office− NERRS SWMP manual https://cdmo.baruch.sc.edu/ request-manuals/). In addition to the SWMP, the Waquoit Bay National Estuarine Research Reserve (WBNERR) coordinates a citizen science water quality monitoring program, the Waquoit BayWatchers (WBW). The WBW program involves 12 to 16 volunteers who monitor 10 locations throughout the Waquoit Bay Estuary. At each site, the volunteers measure standard water quality parameters and collect water samples which are analyzed for a suite of nutrient components and chl a concentration.
The WBNERR monitoring stations at Childs River (CR) and Menauhant (MH) were sampled between May 2018 and August 2019 in coordination with the WBW program (Fig. 1). These 2 dockside stations were sampled every 2 wk in June through September and sampled once a month from October through May, for a total of 23 collection dates. These stations were selected because of the regular sample frequency and the difference in chl a concentration; CR experiences relatively high chl a concentrations during the summer months, while MH generally has low chl a concentrations overall (Howes et al. 2013).
The MH site has an average depth of approximately 2 m at low tide and is located at one of the estuary's southern inlets where the marine water mass from Vineyard Sound exchanges with the estuarine outflow. The average water depth at CR is approximately 1 m at low tide, and the site receives freshwater flow from surface water and groundwater discharge upstream while also experiencing mesohaline tidal influence from Vineyard Sound (Tomasky-Holmes et al. 2013). Due to dense residential developments within the sub-watershed, excess nitrogen from leaching septic systems has resulted in an overabundance of macroalgae, excess phytoplankton production, and thus frequent summer hypoxia events at the CR site (Valiela et al. 1992, D'Avanzo & Kremer 1994, Foster & Fulweiler 2019.

Environmental data
At each station, the WBW volunteers used a YSI Pro2030 sonde to collect water temperature (°C), sa linity, and dissolved oxygen concentrations (mg l −1 ) at the surface (0.25 m below water surface) and bottom (within 2−5 cm of benthic sediments) of the water column. Sampling and measurements occur red before 09:00 h on an outgoing tide, usually within 3 h of slack low tide. For the nutrient and chl a analysis, the volunteers collected 2 l of water in 2 pre-rinsed 1 l acid-washed amber Nalgene bottles with a custom-designed pole sampler. The water samples were filtered onto GF/Fs within 2 h of sampling and frozen at −20 °C until analyzed for nutrient (μM) and chl a (μg l −1 ) concentrations at the Center for Coastal Studies Water Quality Laboratory in Provincetown, MA, following a protocol developed by WBNERR based on standard operating procedures outlined by the EPA (US EPA 2006).
Using the surface and bottom values for water temperature and salinity, a stratification index was calculated according to Eq. (4) in Cohen (1985). The index is between 0 and 1, with 0 being the highest stratification and 1 being no stratification.

Microscopy
In tandem with the volunteers, additional water samples were collected using a 5 l Niskin bottle (General Oceanic) from just below the surface for microscopy analysis of phototrophic taxa at both stations. Triplicate 500 ml water samples were immediately preserved with 5% Lugol's solution in glass amber bottles and sealed with electrical tape. Each 500 ml Lugol's sample was concentrated in the laboratory by settling for 24 h in a 1 l container and then gently pipetting liquid off the top to reduce the total volume. Periodically, this liquid was checked to confirm that organisms were not being removed from the concentrated sample. The final concentrated volume for samples ranged from 9 to 85 ml. Samples collected at CR for 6 sampling dates (13 and 27 Jul 2018, 10 and 25 Aug 2018, and 2 and 16 Aug 2019) were not concentrated because cell abundances were very high. To identify taxa to the lowest classification level and estimate abundances (cells ml −1 ), samples were analyzed with a Zeiss Axiovert S 100 microscope at 400× magnification on a Sedgewick rafter slide (Sherr & Sherr 1993). This magnification allowed us to count all microplankton and larger nanoplankton species (~10−200 μm length), but cells that were too small to be identified were excluded from the counts. A minimum of 300 cells were counted per sample. The focus was on phototrophic taxa, but a single heterotrophic dinoflagellate (Protoperidinium sp.) was included in the counts and analysis.

BrdU-labeled bacterial ingestion experiments
From June 2019 to August 2019, BrdU incorporation experiments were conducted using water collected via Niskin bottles at both sites. BrdU, a thymidine nucleotide analog, is incorporated into prey genomic DNA. When labeled prey are eaten, BrdU is transferred to the grazer genomic DNA via digestion, assimilation, and replication. This material is then recovered by immunoprecipitation and followed by eukaryotic ribosomal amplicon sequencing.
Seawater was taken back to the laboratory and prefiltered through 100 μm Nitex mesh, and triplicate incubations were set up for +BrdU and −BrdU bacterial additions (see Text S1 in the Supplement at www. int-res. com/ articles/ suppl/ m677 p017 _ supp. pdf for details on labeling bacteria with BrdU). The 100 μm seawater (250 ml) was placed into a 500 ml Whirl-Pak bag, and bacteria were added to a final concentration of 10 5 to 10 6 cells ml −1 . Bags were incubated for 48 h at 15°C with a 14 h light:10 h dark cycle. From each incubation bag, 100 ml was collected onto 47 mm, 3 μm Isopore filters, which reduced the amount of bacterial material in the final sample. DNA was isolated following the hot detergent method reported by Gast et al. (2004) using 400 μl of lysis buffer.

Immunoprecipitation of +BrdU DNA and sequence analysis
DNA from +BrdU incubations was immunoprecipitated using an anti-BrdU monoclonal antibody fol-lowing the process reported by Fay et al. (2013) and Gast et al. (2018), with modifications (see Text S1 for details). Amplicons for sequencing were generated from −BrdU DNA and +BrdU immunoprecipitated DNA by PCR amplification of the V4 region of the 18S ribosomal gene using primers 574V4F (5'-[TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG]CGG TAA YTC CAG CTC YV-3') and 1132V4R (5'-[GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G]CCG TCA ATT HCT TYA ART-3'), described in Hugerth et al. (2014) and modified to include 5' adapter sequences (in square brackets) for Illumina MiSeq. Each sample was amplified in triplicate using up to 3 μl template DNA, 1.25 units AmpliTaq Gold 360 DNA polymerase, 2 mM MgCl 2 , 2 μl 2.5 μM dNTPs, and 2.5 μl 10× reaction buffer (25 μl total volume). Amplification conditions were 95°C for 8 min; 40 cycles of 95°C for 30 s, 58°C for 30 s, 72°C for 90 s; 72°C for 5 min; 4°C hold. Replicate reactions were pooled and purified using Qiagen Min-Elute. Purified amplicons were sent to the University of Rhode Island Genomics and Sequencing Center for library preparation and Illumina MiSeq sequencing (300 bp paired end; 600 cycle kit V2).
Only forward reads were analyzed because the fragment was too large (580 bp) to permit efficient assembly with the reverse reads. QIIME2 was used to demultiplex, denoise, remove chimeras, and quality control the forward read Illumina data (see Text S1 for details). The size of the quality-controlled forward read fragment was 200 bp and covered the 5' half of the V4 fragment routinely used for amplicon diversity (Stoeck et al. 2010, Tragin et al. 2018. Al though these fragments were short, the sequence quality was high and contained useful information for taxonomic identification. ASVs were grouped at 100% identity and taxonomy assigned using the Silva 132 database. ASVs identified as bacteria, meta zoa, fungi, and macroalgae were removed from the dataset, as were those that occurred only once (singletons). Raw sequence reads have been de posited at the NCBI Sequence Read Archive (PRJNA690935).
For each experiment, taxa were identified as putative bacterivores based on comparison of tag sequences between +BrdU and −BrdU samples. ASV abundances were converted to a percentage of the total tags for each sample, and the average of each −BrdU ASV was subtracted from the average of the corresponding +BrdU ASV. The ASV was considered a putative bacterivore if the subtracted value was positive and > 0.1% of the total average amplicon abundance and if it was present in at least 2 of the +BrdU samples. If only 1 of the +BrdU samples had the ASV, its abundance needed to be >1% of the total average amplicon abundance to be considered. Putative bacterivores identified as taxa containing chloroplasts were then considered as putative mixotrophs. This approach was based on prior work (Fay et al. 2013), and the use of 0.1 and 1% as cutoffs was to represent the more abundant amplicons in the datasets. This balanced the effect of variability between incubation replicates and reduced the influence of non-specific recovery of extremely abundant DNA (e.g. from diatom taxa). The putative mixotroph ASV taxonomic identifications were compared to the microscopy taxa, and the qualitative identification of a microscopy taxon as a potential mixotroph was made if there were matches at the genus level.

Analysis
ANOSIM (Primer v7; Clarke & Gorley 2015) was applied to the microscopy data to examine whether the data groupings based on factors of site, season, and month were significantly different. Microscopy cell counts were log(x + 1) transformed to balance the contributions of dominant and rare variants. Bray-Curtis resemblance matrices were generated. Oneway unordered ANOSIM was run with 999 permutations for each factor.
We utilized the +BrdU ASV taxonomic information to identify mixotroph taxa in the microscopy samples and their corresponding abundances in 2 ways. First, for only the 2019 summer microscopy samples collected in tandem with water for +BrdU experiments, taxa detected by each +BrdU experiment were identi fied in the corresponding microscopy samples. These are referred to as active mixotrophs because only the taxa that were identified to be grazing when the sample was collected were considered. Second, we applied any taxa detected across all +BrdU experiments to all the microscopy samples. Mixotrophic abundance estimations using the second approach are referred to as potential mixotrophs because although data on whether these taxa were consuming bacteria when the samples were collected were not available, they have demonstrated the ability to consume bacteria in other Waquoit Bay samples. Due to the limitation of microscopy (greater than ~10 μm), estimates of mixotroph abundance were restricted to larger cells.
Poisson generalized linear models (GLMs) with an offset analysis were applied to examine potential environmental drivers of different taxonomic group abundances at each station using microscopy data. The environmental data (independent variables) we considered were water temperature, salinity, turbidity, ammonium, nitrate + nitrite (NO x ), silicate, ortho-phosphate, dissolved inorganic nitrogen (DIN) and total nitrogen (TN) concentration, DIN: dissolved inorganic ortho-phosphate (DIP) ratio, and stratification index (0−1). The offset term was the proportion of the microscope slide that was counted for each sample. The independent variables for each model run were tested for collinearity using the car package in R with the variance inflation factor (VIF) function (Fox & Weisberg 2016). The VIF function tests for variation−inflation in the GLM analysis. Any environmental variable that had a collinearity VIF value of 3 or greater (Zuur et al. 2010) was removed so that analyses were run only with variables that had < 3 VIF. Any variables with p-values greater than 0.05 were also removed until the selected GLM for each taxonomic group at a station had only sig nificant variables with no collinearity. All GLM runs were conducted in R (version 3.6.3) using the built-in linear regression (glm) function.

Environmental data
Most of the environmental data collected by the WBW from May 2018 to August 2019 were significantly different between MH and CR (Table 1). Average salinity, dissolved bottom O 2 concentrations, and stratification index were significantly high er (meaning lower stratification) at MH (Table 1). Average chl a, ammonium, NO x , silicate, TN concentrations, and DIN:DIP ratio were significantly higher at CR (Table 1). Average temperature, turbidity, and phosphate concentrations were not significantly different between MH and CR (Table 1).

Microscopy data
The average (± SE) phototroph abundance at MH was 1189 ± 221 cells ml −1 compared to 7128 ± 2070 cells ml −1 at CR (Fig. 2). Typically, the highest abundances at CR occurred in July and August in both years, with noticeable peaks in October and December 2018 (Fig. 2). At MH, cellular abundances were highest in summer 2018, declined into fall, and then started to increase again in winter 2019 (Fig. 2). For both stations, summer abundances were higher in 2018 compared to 2019 (Fig. 2). Forty taxa were identified over the course of 15 mo; all 40 occurred at MH, and 37 occurred at CR (Table 2). Twenty-five diatoms, 12 dinoflagellates, 1 cryptophyte, 1 chlo rophyte, and 1 euglenid were identified (Table 2). Of the taxa listed in Table 2, all contain chloroplasts except for Protoperidinium. Given the higher total cellular concentration at CR, in dividual taxon abundance was norm alized to total abundance to make comparisons between the 2 stations (Table 2).
Based on the ANOSIM analysis, there was a significant difference in the phototrophic assemblage between MH and CR ( Table 3). The proportion of diatoms was significantly higher at MH, and the proportions of dinoflagellates and euglenids were significantly higher at CR, while the cryptophyte genus, Teleaulax, was consistently prominent at both stations (Table 2, Fig. 3). The phototrophic assemblage was also significantly different be tween seasons and months (Table 3). At CR, there was a strong seasonal variability in the proportion of different taxonomic groups, with cryptophytes and euglenids more prominent in the summer and diatoms dominating in winter and spring (Fig. 3a). There were exceptions to this pattern, with a dinoflagellate bloom in December 2018 and diatoms dominating in late July and early August 2019 (Fig. 3a). At MH, diatoms were consistently more dominant throughout the year, with an occasional increase in cryptophytes (Fig. 3b). The overall proportion of dinoflagellates was higher at CR than at MH, but their proportions remained relatively consistent throughout the year at each station (Fig. 3).
At CR, the 3 most dominant taxa were Teleaulax sp., Leptocylindrus spp., and Eutreptiella gymnastica (Table 2). At MH, the 3 most dominant taxa were Leptocylindrus spp., Teleaulax sp., and Skeletonema spp. (Table 2). Nine diatom taxa composed a significantly higher proportion of the phototrophic assemblage at MH compared to CR, while 4 dinoflagellate and 1 diatom taxa composed a significantly higher proportion of the phototrophic assemblage at CR (Table 2). At CR, the abundances of diatoms, dinoflagellates, and cryptophytes were negatively related to some form of inorganic nitrogen, and all but cryptophytes were positively related to phosphate (Table 4). At MH, the abundances of all taxonomic groups were negatively related to phosphate and positively related to temperature (Table 5). Diatoms were positively related to the stratification index (less stratified), while dinofla gellates and cryptophytes were ne gatively related to the index (Table 5).

Mixotrophy data
Between CR and MH, 53 unique ASVs identified as containing plastids were associated with BrdUlabeled bacterial ingestion over 6 sampling dates (Table 6). More than half of the mixotroph ASVs were shared at both sites, but 9 occurred only at CR and 18 at MH (Table 6). Almost half of the taxa identified as mixotrophs by this process were dinoflagellates, but other common taxa included Microglena uva-maris and Ostreococcus spp. Seven of the ASVs could be matched with 6 genera identified in microscopy samples that were collected simul taneously:  Table 6). The BrdU method could be biased towards the detection of mixotrophic dinoflagellates relative to other mixotrophic taxa due to their high ribosomal gene copy number, but it was possible to detect chlorophyte, cryptophyte, and ochrophyte ASVs, and some were as abundant as the most prevalent dinoflagellate ASVs (e.g. Ostreococcus at >1%). Only a few of the dinoflagellate ASVs occurred at more than 1% in the mixotrophic predictions for each sample. Finally, 98 heterotrophic taxa ASVs were identified by the method (Table S1) : phosphate concentration; TN: total nitrogen; DIN:DIP: (NH 4 + + NO x )/PO 4 3− ratio; S: stratification index. *Significant difference between the 2 stations (paired, 2-tailed t-test, p < 0.05) some were much more abundant than the dinoflagellates (e.g. a jakobid ASV at > 8% and a cercozoan ASV at > 3%). While there may be bias associated with the response of individual organisms to the incubation conditions, this approach allowed the identification of many actively grazing taxa.
Forty-six ASVs that were identified from the BrdU experiments were not associated with taxa from the microscopy samples. These consisted of organisms that were too small to be accurately identified via microscopy (e.g. Ostreococcus spp., Micromonas sp., M. uva-maris), organisms whose identification was too generalized (e.g. Suessiaceae uncultured, Dinophyceae uncultured, Eustigmatales uncultured), or organisms that were not identified in microscopic samples, despite being large enough (e.g. Go ny anulx spinifera, Pelagodinium sp., Paragymnodinium sp.).
We used ASVs from the BrdU experiments to estimate the abundance and proportion of active mixotrophs and potential mixotrophs in the microscopy samples. Active mixo trophs are taxa that were identified to be grazing on bacteria for a specific sample date (Table 6) and then associated with microscopy samples at the sample location and date; this could only be applied to samples collected between June and August 2019 (Fig. 4). At CR, there was an average of 1623 ± 146 active mixotrophs ml −1 (Fig. 4a), which accounted for 34% of the average total num-ber of cells counted (Fig. 4b). At MH, there was an average of 149 ± 33 active mixo trophs ml −1 (Fig. 4a), which account ed for 26% of the average proportion of total cells counted (Fig. 4b).
Potential mixotrophs are taxa associated with the ASVs from all BrdU experiments, regardless of whe ther BrdU experiments were conducted on the sampling date; this was applied to all samples (Fig. 5). There was a seasonal pattern to potential mixotroph abundance, with mixo trophs generally more prevalent in the summer than in the fall, winter, or spring (Fig. 5). Mixotrophs were more abundant at CR, with an average of 4072 ± 1693 potential mixotrophs ml −1 (Fig. 5a), which accounted for 51% of the average proportion of total cells counted (Fig. 5). At MH, there was an average of 301 ± 110 potential mixotrophs ml −1 (Fig. 5a), which accounted for 30% of the average proportion of total cells counted (Fig. 5b). Based on the GLM analysis, the abundance of potential mixotrophs at CR was positively related to turbidity, stratification index, and phosphate concentration and negative ly related to NO x and ammonium concentrations (Table 4). At MH, potential mixotroph abundance was positively related to temperature and salinity and negatively related to stratification index, NO x , and phosphate concentrations ( Table 5). The proportion of potential mixotrophs (to total phototroph cells counted) at CR was positively related to temperature, stratification index, and silicate concentrations and negatively related to turbidity and NO x , ammonium, and phosphate concentrations (Table 4). At MH, the proportion of potential mixotrophs (to total phototroph cells counted) was positively re lated to temperature, stratification index, and NO x and ammonium concentrations and negatively re lated to phosphate concentrations ( Table 5).
The 2 major taxonomic groups that represented the potential mixotrophs were cryptophytes and dinoflagellates, and the proportion of each group varied throughout the year and between stations (Fig. 6). On average, cryptophytes accounted for a higher proportion of the potential mixotroph abundance at both stations compared to dinoflagellates: 64.9 ± 5.6 % at CR and 79.8 ± 0.8% and MH. However, dino flagellates account ed for a higher proportion of the potential mixotroph abundance at CR (35.1 ± 5.6%) compared to MH (20.2 ± 3.6%). In general, cryptophytes were more dominant in the summer, while dinoflagellates were more dominant during winter and spring at both stations (Fig. 6).

DISCUSSION
The microscopy analysis of water samples revealed how the phototrophic assemblage differed between CR and MH. The proportion of diatoms was significantly higher at MH than at CR, while the pro portions of dinoflagellates and euglenids were sig nificantly higher at CR than at MH ( Table 2). The average composition of the phototrophic as semblage at each station was potentially explained by the stability of the wa ter column and nitrogen concentration. Diatoms are known to be better competitors in environments with a more active water column (MH), while dinoflagellates are better competitors in a more stable or stratified water column (CR; Marga lef 1978, Smayda & Reynolds 2001). Additionally, re gions with better water quality, specifically low nutrient concentrations and higher irradiance levels as seen at MH, generally support phototrophic assemblages dominated by diatoms. Conversely, regions with higher nutrients and lower irradiance levels, like those seen at CR, tend to have photosynthetic assemblages dominated by dinoflagellates (Bucha nan 2020).
The seasonal and temporal variability in the photo trophic assemblage at CR generally followed the pattern expected in a temperate estuary, with diatoms more prominent in spring and flagellates more prominent in summer (Marshall et al. 2005). At MH, dia toms tended to dominate throughout the year, but there were periods when cryptophytes would in crease in proportion. Overall, abundances of dia toms, dinoflagellates, and cryptophytes at both stations were po sitively related to temperature, reflecting that total photo troph abundance was typically higher in warmer months. The factors related to variability in dinoflagellate of microscopy data to examine whether the factors of site, season, and month were significant for differences between the phytoplankton communities. MH: Menauhant; CR: Childs River; n/a: not applicable abundance at both stations were what would be expected for a group that was more prominent in the summer: high temperature, low inorganic nitrogen (ammonium), and high stratification (negative relationship with the stratification index) (Smayda & Reynolds 2001, Li et al. 2010, Jauzein et al. 2017). However, dinoflagellates at CR had a positive relationship with the stratification in dex, which implies that their abundan ces were high er when stratification was lower; this was unexpec ted. The relationship with low ammonium concentration may be the result of dinoflagellates having taken up this form of inorganic nitrogen (Halbach et al. 2019). From our analysis at MH, the factors influencing diatom abundance differed from cryptophytes and dinoflagellates in 2 key ways: low er stratification (positively related to stratification index) and NO x concentration. Diatoms prefer low stratification and NO x , so these relationships were expected, again assuming they were taking up NO x (Glibert et al. 2016). At CR, diatom abundance was again related to low NO x concentration but not to stratification. Overall, the stratification results for CR did  Table 4. Results for the Poisson generalized linear model (GLM) with an offset analysis for select groups at Childs River (CR). The factors temperature (Temp), turbidity (NTU), stratification (S), nitrate + nitrite (NO x ), ammonium (NH4 + ), phosphate (PO 4 3− ), and silicate (SiO 2 ) were used in at least 1 of the models. The residual degrees of freedom for all GLMs was 63. (−/+): factor was positively (+) or negatively (−) related to the dependent variable. Blank space: factor was not significantly (p > 0.05) related to variability in the dependent variable or was removed due to collinearity (VIF > 3)   not make sense, likely because of the high variability in this factor. Increased sampling frequency could help improve our statistical strength and elucidate the real significance of stratification on temporal variability at CR. Regarding mixotrophy, our approach of combining BrdU experimental data with microscopy results made it possible to estimate that potential mixotrophs accounted for, on average, one-third to onehalf of the >10 μm phototrophic assemblages and that the spatial difference in those estimations was largely related to the proportion of diatoms and dino-flagellates at each site. Diatoms are the only protistan phytoplankton group with no known mixotrophic species (Flynn et al. 2013). Therefore, the proportion of potential mixotrophs >10 μm is expected to be lower at MH compared to CR since the average proportion of diatoms was significantly higher at MH ( Table 2). Many of the potential mixotrophs identified belonged to taxa in the cryptophyte and dinoflagellate groups. The most prominent group was cryptophytes, but the average proportion of cryptophytes was not significantly different between the 2 stations. This means that while cryptophytes account for the largest proportion of mixotrophs at both stations (Fig. 6), the difference in the proportion of mixotrophs between stations was not driven by cryptophytes. Instead, the microscopy-estimated proportion of mixotrophs was higher at CR compared to MH due to a combination of a lower proportion of diatoms and a higher proportion of dinoflagellates at CR (Table 2). However, looking at just the diatom:dinoflagellate ratio in a sample from any system will not provide complete insight into when and where mixotrophs dominate; all taxa should be considered. The temporal variability of potential mixotroph abundance and proportion was high at both stations (Fig. 5). It has long been hypothesized that mixotrophs dominate the plankton community when either light or nutrients are limiting (Stoecker 1998, Edwards 2019). This appears to describe both CR, a station with sufficient inorganic nitrogen concentration and potentially reduced irradiance levels, and MH, a station with lower nitrogen concentration and potentially sufficient irradiance level (Table 1). Subregions of an estuary with significant freshwater inflow, such as at CR, are associated with high turbidity, sediment and detritus input, and resuspension, all of which reduce irradiance (Kemp et al. 2005). The GLM analysis suggested that potential mixotroph abundance was related to higher turbidity, which would support a light limitation hypothesis. Assuming light was the growth-limiting factor at CR, collection of data on irradiance levels, not just turbidity, might provide more insight into when mixotrophs dominate the photosynthetic community. At MH, po -tential mixotroph abundance was related to higher temperatures and lower NO x concentrations, i.e. summer conditions. Considering that turbidity was lower at MH compared to CR, albeit not significantly (Table 1), nutrient limitation potentially played a larger role than light limitation.
Recent analysis in the southern North Sea suggested that there was little evidence of mixotrophs in turbid eutrophic temperate estuaries (Schneider et al. 2020), similar to conditions in Waquoit Bay. However, our analysis demonstrates that it is possible for mixotrophs to account for a high proportion of the photosynthetic assemblage. Our results do not contradict Schneider et al. (2020), or vice versa, because estuaries are highly heterogeneous, and it can be difficult to directly compare 2 systems. This indicates that more research on mixotrophs across estuaries could help to identify the type that favors mixotrophs and when they are an important part of the food web. For example, Chesapeake Bay is a well-studied estuarine system with robust data on the phototrophic assemblage but lacks data on the abundance and proportion of mixotrophs. Given the well-documented presence of dinoflagellates throughout the year (Mul holland et al. 2018), it would be expected that Chesapeake Bay would favor mixotrophs under appropriate conditions. However, research is needed to identify these conditions.
Our estimation of potential mixotrophs in the microscopy samples has drawbacks, but we believe it is useful in highlighting how much of the photosynthetic assemblage in temperate estuaries can be 28 Fig. 6. Proportion of estimated potential mixotrophs composed of cryptophytes and dinoflagellates at (a) Childs River and (b) Menauhant. Average data point is the average proportion of each group over the course of the study mixotrophs even if they are not actively engaged in their alternative nutrient acquisition mode. It is possible that we were simultaneously under-and overestimating the abundance of different mixotrophs. We were unable to microscopically identify most of the mixotrophic ASVs, potentially underestimating our total number of potential mixotrophs. As previously mentioned, the ASVs that were not identified in microscopic samples consisted of taxa that were too small (10−15 μm) to be identified (such as Ostreococcus spp. and Microglena sp.) and certain dinoflagellate species. This incongruence between microscopy and amplicon data has been noted in many recent studies that have compared the 2 approaches (e.g. Eiler et al. 2013, Gong et al. 2020, Rynearson et al. 2020. The abundances of the small genera were not included in our total abundance of phototrophs; therefore, the identification of potential mixotrophs was biased towards larger species, which tended to be dinoflagellates. This bias may have affected our interpretation regarding the difference in mixotroph proportion between stations, which we related to the balance between diatoms versus dinoflagellates. If we were able to identify and count smaller taxa such as Ostreococcus and Microglena, then our estimation of both potential mixotroph abundance and total photosynthetic cell abundance would likely increase. Even with this increase in the proportion of potential mixotroph abundance, large mixotrophs would still likely dominate the total biomass of mixotrophs because dinoflagellates are known for their high carbon content compared to other taxa groups (Menden-Deuer & Lessard 2000). If the proportion of potential mixotrophs was altered by including smaller species, then our understanding of factors influencing spatial and temporal variability might be different. This bias is an issue not only in our analysis but in any study that uses taxonomybased identification of mixotrophs in microscopy samples. However, our analysis (and previous analysis) is not rendered invalid by this caveat but is merely restricted to larger mixotrophs (>10 μm). It should also be noted that many of these smaller species were not previously reported as mixotrophs, and further laboratory work would help clarify this status.
In the case of larger dinoflagellate-associated ASVs, their lack of identification in microscopy samples was likely the result of misidentification or abundances being too low to be detected. While misidentification is possible, it is most likely to occur with athecate genera, such as Gymnodinium sp. and Gyrodinium sp., because of their less-defined morphologies compared to thecate dinoflagellates. In this case, Gymno-dinium sp. and Gyrodinium sp. were both mixotrophs, so the abundance of potential mixotrophs would not change, but it is possible that this would result in an overestimate if heterotrophic and phototrophic species are included in the genus (e.g. Gyrodinium). If the abundance of a species was too low to be counted in the microscopy samples, then their abundance would be too low to have an impact on the proportion of mixotrophs in a sample.
The largest source of overestimation likely came from a lack of data on which genera were actively grazing in all samples, hence the use of the term potential mixotroph. It was clear from the summer 2019 comparison of microscopy data and BrdU experiments that not every ASV was always grazing on bacteria (Table 6). There was 1 experiment (CR; July 5, 2019) in which none of the mixotrophic ASVs associated with genera in the microscopy samples were grazing (Fig. 4). Additionally, the ASV data for Micro glena uva-maris revealed an interesting pattern of mixotrophic activity; at times of high ASV abundance in the −BrdU samples at CR on June 21 and July 5, it was not identified as an active mixotroph, but it was definitively identified as mixotrophic when at lower ASV abundances. This suggests that under bloom conditions, it may grow primarily as a phototroph, but when conditions are less favorable, it ingests prey. Re gardless of whether a mixotroph was actively grazing when a sample was collec ted, knowledge of when a phototrophic community is dominated by those with the ability to switch nutrient modes is worthwhile. There is an energetic cost associated with being a mixotroph, as they must simultaneously maintain chloroplasts and feeding vacuoles (Stoecker 1998). Even if a mixotroph is not engaged in its alternative nutrient acquisition mode, it should dominate under different conditions than a pure phytoplankter. Therefore, studying potential mixotrophs still provides important information about the type of organisms in the phototrophic community.
From our analysis, it is clear that even with the limited match of 7 ASVs with 6 genera by microscopy, mixotrophic species can be an important part of the phototrophic community in temperate estuaries, at times accounting for up to 98% of the group. Very little is known about the presences of mixotrophs in estuarine systems across different temporal scales (diurnal, tidal, seasonal, annual, etc.). Using taxonomy to estimate the potential abundance of mixotrophs in plankton samples, as demonstrated here and in Schneider et al. (2020), can help to begin expanding our knowledge of this trophic mode in estu-aries. Estimates of potential mixotroph abundance by using plankton taxonomy could be readily accomplished in systems with robust plankton assemblage data, such as Chesapeake Bay. However, our work also demonstrates that estimating the abundance of mixotrophs through taxonomy is imperfect due to 2 primary challenges: not all phagotrophic phototroph taxa are known, and the conditions under which they actively graze are often unclear. We suggest that using BrdU incubations can provide a constraint on estimations of mixotroph abundance by only detecting taxa actively grazing in a given sample. We acknowledge that it is a costly and time-consuming method, but widespread estimations of mixotroph abundance based on taxonomy of plankton samples, supplemented with BrdU incubations, when possible, would rapidly increase our understanding of mixotrophic components.