Temporal and spatial variability of constitutive mixotroph abundance and proportion

Abstract Mixotrophic plankton can comprise a substantial portion of the plankton community compared to phytoplankton and zooplankton. However, there is a gap in the understanding of conditions that influence mixotroph prevalence and activity in situ because current methods often over- or underestimate mixotroph abundance. A labeled prey-tracer method was utilized to identify active mixotrophs present at two locations in a temperate estuary over a year. The tracer method was combined with light microscopy data to estimate active mixotroph abundance and proportion. This study estimated that actively grazing mixotrophic taxa were more abundant in the spring and autumn compared to summer. Dinoflagellates typically dominated the mixotrophic taxa except during autumn at the low salinity location when cryptophytes dominated. Further analysis suggested that active mixotroph abundances might not be only regulated by environmental conditions favorable to mixotrophy but, instead, environmental conditions favorable to different mixotrophs utilization of phagotrophy. By focusing on mixotrophic taxa that were identified to be actively grazing at time of sampling, this study provided a more nuanced estimation of mixotroph abundance, increasing the understanding of how mixotrophic abundance and proportion in situ are influenced by the planktonic community composition and environmental factors.


Introduction
The existence of mixotrophic plankton (organisms that combine photoautotrophy and phagotrophy) has long been recognized (Flynn et al. 2013 ).Ho w ever, it has only been within the past decade that scientists hav e a ppr eciated that mixotr ophs can comprise a substantial portion of the plankton community (Flynn et al. 2013, Stoecker et al. 2017 ).While mixotrophs are known to be pr e v alent, r esearc h on them la gs behind phytoplankton and zooplankton.Recentl y, se v er al r esearc h priorities for mixotr ophs were outlined, one of which was improved categorization of the biogeogr a phical distribution of mixotrophs (Millette et al. 2023 ).Specifically, an understanding of the spatial and temporal patterns of mixotroph abundance and proportion relative to phytoplankton and zooplankton is needed.It is important to know when mixotrophs are expected to be a prominent component of the plankton community because their role in trophic transfer and biogeochemical cycling is distinct from phytoplankton and zooplankton (Mitra et al. 2014, Ward and Follows 2016, Larsson et al. 2022 ).For example, a modeling study demonstrated that a global plankton food-web composed 100% of mixotrophs incr eased carbon sequestr ation by 35% compar ed to a food-web with no mixotrophs (Ward and Follows 2016 ).Unfortunately, without basic knowledge of what proportion of the plankton comm unity is mixotr ophic, it is impossible to assess the accuracy of this result.An important first step in predicting when and where mixotr ophs ar e an important part of the plankton food web is to understand the spatial and tempor al v ariability of mixotr oph abundance and the environmental factors associated with their variability.
Ho w e v er, ther e is pr esentl y a methodological limitation in the ability to differentiate mixotrophs from phytoplankton and zooplankton in field studies (Millette et al. 2018 ), which limits our ability to study the biogeogr a phical distribution of mixotrophs.For example, mixotr ophs ar e typicall y gr ouped with phytoplankton due to the use of photopigments (c hlor ophyll a ) as an indicator of photosynthetic ca pability.Furthermor e, gr azing activity is often attributed to only zooplankton because it is indirectly measured by looking at the change in prey concentration, rather than identifying organisms doing the grazing (Millette et al. 2018 ).Given how pr e v alent mixotr oph taxa can be in the plankton community (Domaizon et al. 2003, Unrein et al. 2007, Millette et al. 2021, Li et al. 2023 ), r esearc h needs to transition toward identifying and studying mixotrophs as their own plankton gr oup(s), a part fr om phytoplankton and zooplankton.
The ability to estimate mixotroph abundance and activity in situ is hindered by the limitations of current popular methods (Anderson et al. 2017, Li et al. 2021 ).Activ e mixotr oph abundance is typically estimated with fluorescently labeled bacteria or fluor escent micr ospher es (Ar enovski et al. 1995, Sanders et al. 2000, Czypionka et al. 2011, Anderson et al. 2017, Millette et al. 2017, Sato et al. 2017, Gast et al. 2018, Li et al. 2023 ).These methods w ork b y estimating the abundance of cells that consume fluor escentl y labeled material and contain c hlor oplasts (González 1999 ).Ho w e v er, it has been demonstrated that this a ppr oac h c hr onicall y underestimates mixotroph abundance (Anderson et al. 2017, Li et al. 2021 ) because these experiments have a short incubation period ( < 2 h) that likely does not ca ptur e all gr azing activity.Furthermore, certain species might be biased against fluor escentl y labeled particles, while others may be ov erl y biased tow ar ds them (Sanders andGast 2012 , Wilken et al. 2019 ).Recent studies have also attempted to estimate mixotroph abundance in publicl y av ailable micr oscopy datasets via taxonomy (Har a guc hi et al. 2018, Leles et al. 2019, Cesar-Ribeiro et al. 2020, Schneider et al. 2020, Li et al. 2023, Mitra et al. 2023 ).This a ppr oac h uses av ailable micr oscopy-based taxonomic data to estimate the abundance and proportion of potential mixotrophs, plankton that have been found to be capable of photoautotrophy and phagotrophy in pr e vious peer-r e vie wed studies , o v er lar ge tempor al and spatial scales.Ho w e v er, this method is likely an overestimation of the presence of mixotrophs, since it is not known if a species is utilizing mixotrophy in any given sample or location.
Recentl y, br omodeoxyuridine (BrdU) has been utilized to identify activ el y gr azing mixotr ophic taxa in field samples (Fay et al. 2013, Millette et al. 2021 ).This method targets constitutiv e mixotr ophs (CMs), or ganisms that activ el y synthesize and maintain their own c hlor oplasts, that ingest bacteria (Ghyoot et al. 2017 ).While this a ppr oac h does not ca ptur e all possible mixotr ophs pr esent in a system, it provides taxonomic information on the CMs activ el y ingesting bacteria at any given time, whic h ar e an important and often dominant type of mixotroph (Jeong et al. 2010, Edw ar ds 2019, Faure et al. 2019 ).In studies using BrdU-labeled prey, the assumptions are that mixotrophs will ingest and digest the prey, and then incorporate the labeled nucleotides in their o wn DN A during replication.Amplicon sequencing of selectiv el y r ecov er ed DN A allo ws identification of active mixotrophs in the phototrophic taxa.Those mixotrophs can be matched with taxa identified in microscopy samples to estimate the abundance and proportion of mixotrophs (Millette et al. 2021 ).So, while this a ppr oac h does not ca ptur e the full mixotr ophic community, it does capture the majority of CMs in microscopy samples.A recent study by Millette et al. ( 2021 ) utilized the BrdU method to identify and quantify CMs in microscopy samples from Waquoit Bay, MA and suggested that CMs could account for over 90% of c hlor oplast containing plankton in a sample.Ho w e v er, Br dU experiments w er e onl y conducted for six of the 23 sample dates, so most of the conclusive results were based on potential CM abundance, rather than which taxa were actively grazing in each sample.Millette et al. ( 2021 ) demonstrated that estimating CM abundance this way was likely an ov er estimation because the presence of a mixotroph does not guarantee that a mixotroph is using its alternate nutrient mode (grazing).Ho w ever, if the BrdU method was utilized e v ery time a microscopy sample is collected, then the abundance of mixotrophs actively grazing in those samples could be estimated.
In this study, BrdU experiments were conducted for all sampling dates and used to constrain estimates of CM abundance in microscopy samples to taxa actively grazing.This made it possible to taxonomically identify and estimate CM abundances to assess the variability of CMs that ingest bacteria within a temperate estuary across a whole year.The goal was to investigate the biotic and abiotic factors associated with the temporal and spatial variability of active CMs identifiable though microscopy.To accomplish this goal, there were three objectives: (1) to develop lists of CMs grazing on bacteria in two distinct parts of the York River Estuary , Chesapeake Bay , USA, (2) compare the estimations of potential versus active CM abundance and proportion, and (3) use the active CM data to investigate the biotic and abiotic factors asso-ciated with the spatial and tempor al v ariability of CM abundance and proportion.This study identified over 130 mixotroph taxa using the BrdU method that were matched up with seven taxa identified in microscopy samples.Results suggested that while CM abundance was related to environmental conditions, the abundance of active CMs was also related to the taxa present.

Stations
In total, two stations were sampled twice a month over 1 year in the York River Estuary, Virginia, USA; West Point (WP) and Gloucester Point (GP; Fig. 1 ).WP is located up-estuary, near one of the two major freshwater sources to the York River (Mattaponi River) at the WP fishing pier.GP is located closer to the mouth of the estuary, where the river flows into Chesapeake Bay, at the GP fishing pier.The WP station is c har acterized b y lo w salinity (oligo-to mesohaline) with high nutrient concentrations and high turbidity, as it is located near the estuarine turbidity maximum (Friedrichs 2009, Reay 2009 ).The GP station is c har acterized by high salinity (meso-to polyhaline) with low nutrient concentrations and high turbidity compared to WP (Friedrichs 2009, Reay 2009 ).

En vironmental da ta
All sampling occurred on an incoming tide, no later than 1 h before high tide.WP was always sampled first, followed by GP.A YSI EXO1 sonde (Xylem Inc.) was used to conduct a vertical profile of temper atur e ( • C), salinity, and turbidity (FNU: Formazin Nephelometric Unit).The data points for each profile were averaged for e v ery 0.1 m and the av er a ged data points at 0.5 m wer e used for analyses.A LI-1400 (2 π quantum sensors; Deck: LI-190SA and Water: LI-192SA) from LI-COR was used to conduct an irradiance profile of the water column and calculate the light attenuation coefficient (k d ) at a depth of 0.5 m.
A 5-l Niskin bottle (General Oceanics) was used to collect water from just below the surface.Water for nutrient analysis (30 ml of 0.45 μm-filtered water) was collected in two 20 ml acid-washed plastic vials.Once in the laboratory, no more than 2 h after samples were collected, the nutrient samples wer e fr ozen at −20 • C until analysis for ammonia (NH 3 ), nitrate + nitrite (NO x ), phosphate (PO 4 3 − ), and silica (SiO 2 ) ( μM) using a Skalar Auto Analyzer in the Virginia Institute of Marine Science Analytical Services Center (U.S. EPA 1993aEPA , 1993bEPA , 1993cEPA , 1997 ) ). Water for c hlor ophyll a analysis was collected in three 1 l clear Nalgene bottles .T he bottles were k e pt in a dark cooler until they were brought back to the labor atory.To measur e c hlor ophyll a concentr ations, 40-150 ml of the water collected was filtered onto 25 mm GF/Fs .T he filters were resuspended in 7 ml of 90% acetone for 24 h and placed in the freezer at −20 • C.After 24 h, the samples were read on a Turner Designs 10-AU fluorometer for chlorophyll a concentrations before and after acidification with 10% HCl (Arar and Collins 1997 ).

Microscopy
Additional water from the 5-l Niskin cast was collected in three 500 ml Nalgene amber bottles, immediately preserved with 5% Lugol's solution, and sealed with electrical tape .T he Lugol's samples were concentrated in the laboratory by settling for 24 h in a 500-ml beaker and then gently pipetting liquid off the top to reduce the total volume to ∼50 ml.Chloroplast containing plankton genera/species in these samples were identified and enumerated (cells ml −1 ) with a Zeiss Axio Imager.A2 microscope at 400x magnification on a Sedgewick rafter slide (Sherr and Sherr 1993 ).
A minimum of 300 cells were counted per sample.

BrdU-labeled bacterial ingestion experiments
To identify activ el y gr azing CMs pr esent at eac h station on eac h sampling date, incubation experiments were conducted using BrdU-labeled bacteria as prey.Two cultures of Photobacterium angustum (heter otr ophic bacteria) wer e gr o wn in y east extract 10 days before sampling.BrdU (800 μM final concentration) was then added to one of the bacterial cultures ( + BrdU) 3 days before sampling (Millette et al. 2021 ).The Br dU w as w ashed off the culture the day before sampling by centrifuging at 3000 r m −1 for 10 min at 4 • C. The supernatant was r emov ed without disturbing the bacterial pellet.A volume of 10 ml of sterile 1x phosphate-buffered saline (PBS) was used to resuspend pellets and they were centrifuged again.This was repeated three times and after the final wash, bacterial pellets wer e r esuspended in a total of 20 ml sterile 1x PBS and cell concentration determined using hemocytometer.The bacterial culture without Br dU ( −Br dU) w as also "w ashed" three times to ensure both cultures were treated the same .T he water collected from the Niskin cast was used to set up triplicate incubations for + BrdU and −BrdU bacterial additions.A volume of 250 ml of water was placed into a 250-ml clear Nalgene bottle and bacteria was added to a final concentration of 10 6 cells ml −1 .The bottles were then incubated for 24 h in one-layer mesh bags in the York River Estuary, near the Virginia Institute of Marine Science.At the end of the experiment, all water from each incubation bottle was collected onto 47 mm 3-μm Isopore filters and stored at −20 • C until analyzed.

Immunoprecipitation of ±BrdU DNA and sequence analysis
Nucleic acids were extracted for all samples following the hot detergent method reported by Gast et al. ( 2004) using 200 μl of lysis buffer.Extracted material from + BrdU incubations then went Br dU-labeled DN A samples from the experiments (dilution of 300 ng in 10 ul) were denatured at 95 • C for 10 min, placed on ice for 2 min, and then 10 μl of cold PBS-BSA was added.Each denatured Br dU DN A sample w as then combined with 20 μl of the blocked antibody mix and incubated overnight at 4 • C on a rotating mixer.This step allo w ed for BrdU-labeled DNA to bind to the anti-BrdU mouse antibody.The following day 40 μl of blocked beads were added to each antibody sample and incubated for at least 1 h at 4 • C on a rotating mixer to allow the BrdU-labeled DNA bound to the mouse antibody to bind to the beads .To remo ve non-BrdUlabeled DNA from the antibody/bead mixture, tubes were transferred to a magnet (DynaMag TM -2, Invitrogen 12321D) for about 5 min.Once the solution was clear, it was r emov ed and discarded, and the beads were fully resuspended in 1 ml of washing buffer (PBS with 0.05% Tween 20).The solution was allowed to clear on the magnet, follo w ed b y r emov al of liquid without bead pellet disruption.This step was repeated three times .T hen the beads were resuspended in 0.5 ml of washing buffer, the solution allo w ed to clear on the magnet, follo w ed b y r emov al of liquid without bead pellet disruption, repeated four times.To recover BrdU-labeled DNA from the beads, 100 μl of elution buffer (5 mM BrdU in PBS; Sigma B9285-250 mg) was added to each sample and incubated at 65 • C for 20 min and then at r oom temper atur e for 15 min.The solution was allowed to clear on the magnet and the liquid tr ansferr ed to a new tube for ethanol precipitation (90 μl 100% isopropanol, 35 μl 3 M NaCl, 100 μl elution) over night at −20 • C.After centrifugation, the DNA pellet was allowed to dry for se v er al minutes before resuspending in 10 μl PCR water.
Amplicons for −Br dU DN A and imm unopr ecipitated + BrdU DN A w er e gener ated thr ough PCR amplification of the V4 region of the 18S ribosomal gene using primers V418SF (5 [TCGTCGGC AGCGTC AGA TGTGT A T AAGAGAC AG] CC AGC AS-CYGCGGT AA TTCC) and V418SR (5 [GTCTCGTGGGCTCGGAGAT-GTGT A T AA GA GA CA G] A CTTTCGTTCTTGA TYRA TGA), described in Piredda et al. ( 2017 ) and modified to include 5 adapter sequences (in square brackets) for Illumina MiSeq.Each sample was amplified in triplicates using up to 3 μl template DNA, 1.25 units GoTaq Flexi DNA pol ymer ase, 2 mM MgCl 2 , 2 μl 2.5 μM dNTPs, and 2.5 μl 10x reaction buffer (25 μl total volume) with the conditions: 95 • C for 8 min; 35 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.Some of the -BrdU samples needed to be diluted in order to be amplified ( Table S2 , Supporting Information ).Amplicons were sent to the Rhode Island Institutional De v elopment Aw ar d (IDeA) Netw ork of Biomedical Research Excellence Molecular Informatics Core for libr ary pr epar ation and Illumina MiSeq (300 bp paired end; 600 cycle kit V2) sequencing.
QIIME2 was used to demultiplex, denoise, remo ve chimeras , and quality control the Illumina data.Amplicon sequence variants (ASVs) wer e gr ouped at 100% identity and taxonomy assigned using the Silva 132 database.ASVs identified as bacteria, metazoan, fungi, and macroalgae were removed from the analyzed dataset, as were those that occurred only once (singletons).For each experiment, taxa were identified as 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 av er a ge of each −Br dU ASV w as subtr acted fr om the av er a ge of the corr esponding + BrdU ASV.An ASV was considered a bacterivore if the subtracted value was positive and > 0.1% of the av er a ge total amplicon abundance, suggesting that this ASV was present at a higher proportion in the + BrdU samples.Bacterivores identified as taxa containing c hlor oplasts wer e then consider ed activ el y gr azing CMs (Fay et al. 2013, Millette et al. 2021 ).The use of this conserv ativ e a ppr oac h was to r epr esent the mor e abundant amplicons in the datasets, balance the effect of variability between incubation replicates and reduce the influence of nonspecific r ecov ery of extr emel y abundant DNA (e.g. from diatoms; Millette et al. 2021 ).

Potential and acti v e CM abundance and proportion
The abundance of potential CMs was calculated by summing the abundance of genera present in each microscopy sample that had been identified as a CM in the BrdU experiments for any sampling date at a specific location.This calculation reflects the abundance of plankton with the ability to engage in mixotrophy but not whether they were activ el y ingesting bacteria prey, or using their alternate nutrient mode, as done in Millette et al. ( 2021 ).To calculate active CM abundance, a qualita-tive identification of a microscopic taxon as an active CM was made if there were matches at the genus le v el for eac h sampling date at each station.This calculation reflects the abundance of CMs that were actively engaging in mixotrophy at the time the sample was collected.The proportion of potential and active CMs were calculated by dividing the total abundance of taxa identified as either a potential or active CM by the total phototroph abundance for that sample .P otential CM abundance and pr oportion wer e compar ed to activ e CM abundance and pr oportion to assess potential ov er estimations in the presence of CMs when not accounting for which taxa are actively ingesting.It was assumed that any CM taxa that was present but not identified to be activ el y ingesting bacteria when a sample w as collected w as onl y photosynthesizing.Further anal ysis on biotic and abiotic factors associated with the abundance and proportion of CMs was only done on active CMs.

Analysis
Quasipoisson generalized linear models with an offset (GLMs) wer e a pplied for both stations to examine envir onmental conditions associated with the abundance of diatoms , dinoflagellates , cryptoph ytes, haptoph ytes, and active CMs.Quasibinomial GLMs with an offset wer e a pplied for both stations to examine environmental conditions associated with the proportion of active CMs.Triplicate samples were k e pt se parate and sample dates with any missing data were removed from analysis (WP: n = 57, GP: n = 48).The independent v ariables (envir onmental data) included wer e temper atur e, salinity , turbidity , attenuation coefficient (K d ), and nutrient concentrations (NO x , NH 3 , PO 4 3 − , and SiO 2 ).The offset term was the number of grids on the microscope slide that was counted for each sample .T he independent variables for each model run were tested for collinearity using the "car" pac ka ge in R with the variance inflation factor (VIF) function (Fox and Weisber g 2019 ).An y envir onmental v ariable that had a collinearity VIF value of 5 or greater (Zuur et al. 2010 ) was removed so that anal yses wer e run onl y with independent variables that were not str ongl y r elated to eac h other.An y v ariables in the selected GLM with P -values less than .5 were considered to have some signal and it was reported whether this signal was positiv el y or negativ el y r elated to the dependent variable.All GLM runs were conducted in R (version 3.6.3)using the built-in linear regression (glm) function.

En vironmental da ta
T he a v er a ge physical and environmental data collected from March 2021 to February 2022 are presented in Table 1 .The aver a ge turbidity, temper atur e, and K d wer e significantl y higher at WP while the av er a ge salinity was significantl y higher at GP. Howe v er, the majority of biological (c hlor ophyll a concentr ations) and c hemical (nutrient concentr ations) data wer e not differ ent between WP and GP.The only exception to this was NO x concentr ation, whic h was significantl y higher at WP compar ed to GP. Gr a phs of the unav er a ged time series of measurements illustrate the variability of environmental factors within and between stations ( Figures S1 -S9 , Supporting Information ).

Micr oscopy-based phototr oph abundance and major taxonomic groups
T he a v er a ge ( ±SE) phototr oph abundance based upon micr oscopy was not significantly different between WP (3777 ± 305 cells ml −1 ) and GP (4883 ± 367 cells ml −1 , P -value = .43,t -test, Fig. 2 ).Ho w e v er, the tempor al v ariability of phototr oph abundance was high at both stations .T he highest phototroph abundances at WP occurred during the late summer, during the months of August and September, and remained elevated through November (Fig. 2 ).At GP, phototroph abundances were highest during the winter, in the months of January and February due to a small diatom bloom dominated by Skeletonema spp.(Fig. 2 ).The proportion of major taxa gr oups pr esent between stations and sampling dates was also highly variable.Ho w ever, the general trend for both stations was dinoflagellates dominating during the spring, diatoms during the summer, cryptophytes during the autumn, and then diatoms again during the winter (Fig. 3 A and B).At GP, the shift to w ar d diatoms dominating again during the winter was specifically due to the aforementioned Skeletonema spp.bloom (Fig. 3 B).
Based on results from the GLM analysis, at WP, diatom abundance was positiv el y corr elated to temper atur e, NH 3 and SiO 2 , and negativ el y corr elated to turbidity, NO x , and K d (Table 2 ).Dinoflagellate abundance was positiv el y corr elated to K d and negativ el y correlated to turbidity and NO x (Table 2 ).Cryptophyte abundance was positiv el y corr elated to temper atur e, NH 3 , and SiO 2 , and negativ el y corr elated to PO 4 3 − and K d (Table 2 ).Ha ptophyte abundance was positiv el y corr elated to temper atur e, PO 4 3 − , SiO 2 , and K d , and negativ el y corr elated to NH 3 (Table 2 ).At GP, diatom abundance was positiv el y corr elated to salinity, SiO 2 , and K d , and was negativ el y corr elated to temper atur e, NO x , and PO 4 3 − (Table 2 ).Dinoflagellate abundance was positively correlated to temperature and K d , and was negativ el y corr elated to salinity, NH 3 , NO x , PO 4 3 − , and SiO 2 (Table 2 ).Cryptophyte abundance was positiv el y corr elated to NH 3 , SiO 2 , and K d , and was negativ el y corr elated to temper atur e, NO x , and PO 4 3 − (Table 2 ).Haptophyte abundance was positiv el y corr elated to K d , and was negativ el y corr elated to temper atur e, salinity, NH 3 , NO x , PO 4 3 − , and SiO 2 (Table 2 ).

Mixotrophic ASVs
A total of 138 unique ASVs were identified as containing plastids and associated with the ingestion of BrdU-labeled bacteria (CMs) over the 24 sampling dates between WP and GP.A total of 109 of these ASVs were identified to at least the genus le v el.Out of the mixotrophic ASVs identified, 47 were unique to WP, 32 were unique to GP, and 59 occurred at both stations.Dinoflagellates ASVs made up the lar gest pr oportion ( ∼44%; at GP, dinofla gellates comprised ∼53% of the ASVs; Fig. 4 ) of all major taxa groups and were the most evenly distributed between the two stations.From the 61 dinoflagellate ASVs identified, 33 occurred at both stations, while 12 were unique to WP and 16 were unique to GP. Chrysoph ytes and cryptoph ytes were the least e v enl y distributed between stations.A total of 12 chrysophyte and 10 cryptophyte ASVs were unique to WP, while only one chrysophyte and two cryptophyte ASVs were unique to GP.At WP, while dinoflagellates (43%), chrysophytes (15%), and cryptophytes (14%) also accounted for a notable proportion of ASVs (Fig. 4 ).At both stations, the number of CM ASVs identified for a sampling date was highest in early spring.There was a decrease during the summer and a small increase during the autumn (Fig. 5 ).Dinoflagellate ASVs were the dominant mixotrophic taxa group throughout the whole year at GP. Ho w e v er, at WP, dinofla gellate ASVs dominated in the first half of the year with cryptophytes becoming dominant later in the year (Fig. 5 ).Chrysophyte ASVs were al ways pr ominent at WP, but ne v er dominant.uncultured marine eukaryote) (Fig. 6 ).The remaining 116 ASVs that were identified from the BrdU experiments were not associated with taxa from the microscopy samples.Most of those taxa were either too small to be accur atel y identified thr ough microscop y, or w ere ASVs with general identification (e.g.chrysoph yceae uncultured eukaryote, dinoph yceae uncultured eukaryote, and cryptophyceae uncultur ed fr eshw ater eukary ote).It is also possible that a number of these taxa were too r ar e within the system to be ca ptur ed via our light microscopy analysis.

Potential and acti v e CM abundance and proportion
At WP, the av er a ge abundance ( ±SE) of potential CMs was 1221 ± 132 cells ml −1 (Fig. 7 A), which accounted for 40% ± 3.6 of the aver a ge total number of phototrophic cells counted (Fig. 7 C).At GP, the av er a ge abundance of potential CMs was 1387 ± 151 cells ml −1 (Fig. 7 A), which accounted for 44% ± 2.6 of the av er a ge pr oportion of total phototrophic cells counted (Fig. 7 C).The abundance and proportion of potential CMs was highly variable throughout the year at both the GP (376-2824 cells ml −1 ; 4%-74%) and WP (419-3422 cells ml −1 ; 12%-83%) stations .T he abundance and proportion of phototrophic cells that were potential CMs at GP was highest during spring and autumn, while at WP, abundance and proportion of phototrophic cells that were potential CMs was highest during the spring and winter.
At WP, the av er a ge abundance ( ±SE) of activ e CMs was 388 ± 56 cells ml −1 (Fig. 7 B), which accounted for 13% ± 1.6 of the av er a ge total number of phototrophic cells counted (Fig. 7 D).At GP, the aver a ge abundance of active CMs was 274 ± 41 cells ml −1 (Fig. 7 B), which accounted for 7% ± 0.7 of the average proportion of total phototrophic cells counted (Fig. 7 D).The abundance and proportion of active CMs was highly variable throughout the year at both the GP (0-1297 cells ml −1 ; 0%-38%) and WP (0-1270 cells ml −1 ; 0%-40%) stations .T he abundance and pr oportion of phototr ophic cells that were active CMs at GP was highest during summer, while at WP, abundance and proportion of phototrophic cells that were active CMs was highest during the autumn and winter.
Cryptophytes and dinoflagellates were the only two taxonomic gr oups r epr esented in the potential and active CM abundance.At WP, there was seasonal variability in the dominant taxonomic gr oup r epr esenting activ e CMs; during the spring, dinoflagellates such as Gymnodinium spp., H. rotundata , and Scrippsiella spp.dominated the community of active CMs, while the cryptophyte Teleaulax sp. was more frequently identified as an active CM throughout the autumn (Fig. 6 A).At GP, dinoflagellate ASVs such as Gymnodinium spp.and H. rotundata were frequently identified as active CMs throughout the whole year in microscopy samples (Fig. 6 B).
Based on the GLM analysis, the abundance of active CMs at WP was positiv el y corr elated to turbidity, NH 3 , and SiO 2 concentrations, and negativ el y corr elated to temper atur e, NO x concentr ations, and K d (Table 2 ).The proportion of active CMs was positiv el y corr elated to turbidity, NH 3 , and SiO 2 concentr ations, and negativ el y corr elated to temper atur e, NO x concentr ations, and K d (Table 2 ).The abundance of active CMs at GP was positiv el y correlated to K d and PO 4 3 − concentrations, and negatively correlated to temper atur e and NH 3 and NO x concentr ations (Table 2 ).The pr oportion of activ e CMs was positiv el y corr elated to temper atur e and PO 4 3 − concentrations, and negatively correlated to salinity, NH 3 , NO x , SiO 2 concentrations, and K d (Table 2 ).

Discussion
The BrdU method was used to identify over 130 CMs that were activ el y ingesting bacteria within an estuarine system across 1 year.The BrdU data was combined with microscopy taxonomy and enumeration data to estimate the abundance and proportion of potential and active CMs that could be identified micr oscopicall y.These r esults demonstr ated that estimations of CM abundance were higher for potential CMs compared to active CMs, suggesting that the BrdU experiments can constrain estimations of abundance when used to account for taxa known to be activ el y gr azing in a given microscopy sample.Estimations of active CMs were further used in an assessment of what biotic and abiotic factors are associated with temporal and spatial variability of active CMs.Specifically, it was identified how the activ el y gr azing mixotr ophic taxa c hanged thr oughout the year and that there were distinct differences in the mixotrophic activity of some taxa at the different stations .T he results suggested that active mixotroph abundance is not only regulated by environmental conditions generally favorable to mixotrophy but environmental conditions favorable to specific taxon's utilization of pha gotr ophy.

Potential versus active CMs analysis
Measurements of potential CM abundance compared to the abundance of active CMs were on av er a ge 37% higher at GP and 27% higher at WP (Fig. 7 C).The primary reason for this difference was the presence of the cryptophyte Teleaulax sp.throughout the year at both stations.While Teleaulax sp. was continually present, it was primarily identified to be actively grazing at WP during the autumn (Fig. 6 A).Estimating potential CM abundance assumes that if a taxon has the capacity to ingest prey, that it will activ el y be ingesting all the time .T hese results indicate that at least some mixotrophic taxa are likely only using their alternate trophic mode under certain conditions .T her efor e, estimating the abundance of CMs based on who has the capacity to ingest prey can ov er estimate the pr oportion of the planktonic community that is engaged in mixotrophy at a given time.For example, a recent study that also estimated potential CM abundance and proportion in a temperate estuary (Waquoit Bay, MA, USA) using the BrdU method reported that cryptophytes wer e consistentl y a dominant part of the mixotrophic population (Millette et al. 2021 ).
The findings of this current study suggest that Millette et al. ( 2021 ) ma y ha ve o verestimated the importance of cryptophytes in the CM assemblage and underestimated the importance of dinoflagellates by focusing on the capacity for grazing.

Temporal and spatial variability in active CMs
On av er a ge, the abundance and proportion of CMs between the tw o stations w as not significantl y differ ent (Fig. 7 B and D).This result was not unexpected because mixotrophs are hypothesized to be dominant when one growth factor is limiting (Stoecker 1998 ), and each station likely had one limiting condition.The av er a ge water quality data indicated that light and nutrients wer e significantl y differ ent between WP and GP; light le v els (K d ) w ere lo w er at WP compared to GP, while NO x was significantly higher at WP compared to GP.Average NH 3 and PO 4 3 − concentr ations wer e also higher at WP, although these differ ences wer e not significant due to high annual variability ( Figures S3 and S5 ,Supporting Information ).Ov er all, this aligns with what is known about the av er a ge envir onmental conditions at these stations, one likely being light limited for phototrophic growth (WP) and the other likely being nutrient limited for phototrophic growth (GP, Friedrichs 2009, Reay 2009 ).Since both stations had one likely lim-iting factor, it was expected that the importance of CM abundance and proportion would be equally important at both stations, albeit for different reasons.
It was, ther efor e, expected that the temporal variability in active CM abundance within each station would vary based on the primary growth limiting factor, light at WP and nutrients at GP. Ho w e v er, r esults fr om the GLM anal ysis wer e inconclusiv e r elated to these factors.At the WP station, high abundance and proportion of active CMs was associated with low nutrients (low NO x ) and high nutrients (high NH 3 ), and high light (low K d ).This possibly suggested that actively grazing CMs are more prevalent at WP when light and nutrients are less limiting, and not negativ el y r elated to light, we hypothesized.At GP, a high abundance of active CMs was associated with both low nutrient (low NH 3 and NO x ) and low light (high K d ), suggesting that activ el y gr azing CMs ar e mor e pr e v alent when light and nutrients are more limiting.Ho w ever, a high proportion of active CMs was associated with low nutrients (low NH 3 and NO x ) and high light (low K d ), suggesting that activ el y gr azing CMs ar e mor e pr e v alent when nutrients ar e mor e limiting, which is the result hypothesized.Even though the proportion of active CMs at GP was related to the one environmental factor expected, most GLM analysis results do not align with the hypothesis that CMs have an adv anta ge when onl y one gr owth factor is limiting (Stoecker 1998, Edw ar ds 2019, Li et al. 2023 ).Still, there is evidence that the differences in environmental conditions between the stations might influence the av er a ge abundance of CMs.Ho w e v er, further anal ysis suggested that factors were associated with the variability in active CM abundance at each station.
The abundance of CMs at both stations was correlated to the abundance of major taxonomic groups.At GP, the abundance of active CMs was correlated to dinoflagellate abundance ( Table S3 , Supporting Information ).At WP, active CM abundance for the whole dataset was not correlated with any major taxa group, but active CM abundance for the first 6 months of sampling was correlated with dinoflagellate abundance and cryptophyte abundance for the second 6 months ( Table S3 , Supporting Information ).This suggested that when more dinoflagellates were present at GP and either mor e dinofla gellates or cryptophytes wer e pr esent at WP, depending upon the time of year, there would be mor e activ e CMs.Furthermore , the en vironmental patterns identified to be associated with a high abundance and proportion of CMs using GLM anal yses wer e similar with the patterns associated with a high abundance of dinoflagellates at GP and patterns associated with high abundances of cryptophytes at WP (Table 2 ).This suggested that the abundance of active CMs throughout the year at a specific location was related to both the environmental conditions that favored the growth of particular taxa and each taxon's utilization of pha gotr ophy.This c hallenges the assumption that mixotr oph abundance equals mixotrophic activity and means that understanding mixotrophy on a large scale r equir es assessing the conditions that favor mixotrophic activity in the dominant taxon present.As with the Teleaulax sp. in this study, the cryptophyte was often present in the microscopy samples but only identified to be activ el y gr azing bacteria at specific times and locations.Identifying dominant mixotrophic taxa present within a certain region and the conditions that cause them to switch nutrient modes are important to understanding how mixotrophy contributes to nutrient and carbon trophic transfer.The BrdU prey-labeling experiments used in this study provided a starting point for conducting this analysis.

BrdU method
The BrdU method used in this study is still r elativ el y ne w, and while it has potential, it also has dr awbac ks.First, it is costl y and time-consuming compared to most basic mixotrophy detection methods .T he imm unopr ecipitation pr ocess takes 3 days of activ e and pr ecise work, and onl y a limited number of samples can be processed at a given time.Depending on the number of samples, the processing of samples could potentially take months.For all that effort, this method exclusiv el y tar gets mixotr ophs with their own c hlor oplasts that are ingesting bacteria (CM), excluding other types of mixotrophs.Ho w ever, CMs are an important and ubiquitous mixotrophic group that commonly ingest bacteria, in addition to small protists (Jeong et al. 2010, Faure et al. 2019, Li et al. 2023, Mitra et al. 2023 ), which means this approach likel y ca ptur es the majority of CMs pr esent in a system.Furthermore, the BrdU method is targeting the same mixotrophs as most fluor escentl y labeled prey experiments, which is how most CMs ar e curr entl y identified in situ .This method might also be biased to w ar d the detection of mixotrophic dinoflagellates relative to other mixotrophic taxa due to their high ribosomal gene copy number (Millette et al. 2021 ), although a substantial number of other mixotrophic ASVs such as chloroph ytes, cryptoph ytes, and c hrysophytes wer e also identified.Furthermor e, taxa pr esent in high abundances such as diatoms (nonmixotrophic), might lead to nonspecific r ecov ery of DNA, undermining confidence in species classified as mixotrophs (Flynn et al. 2013 ).Also, it is often difficult to micr oscopicall y identify man y of the mixotrophic ASVs, potentiall y under estimating abundance of activ e CMs.Ho w e v er, most of the ASVs not identified in the microscopic samples consisted of taxa that were too small ( < 10 μm) to be identified by microscopy and were not included in reported active CM abundances or total phototroph abundances .T herefore , the active CM abundance data should be considered a representation of CM > 10-15 μm.
Although imperfect, this method identifies CM taxa in a sample activ el y ingesting bacteria.T hus , allowing for the creation of a list of active CMs that ingest bacteria within the environment being studied.This method not only allo w ed identification of over 130 CMs (80% to the genus le v el), but it also resulted in the creation of two separate lists of CMs unique to distinct parts of a temperate estuary along with correlation of the environmental conditions that supported their active ingestion of prey.By adjusting the sampling frequency or stations, future studies using this method could be expanded to examine the time of year or conditions that different CM ASVs are identified to be actively gr azing compar ed to when they are present and not grazing.This would substantiall y impr ov e estimations of the abundance and proportion of potential CMs from long-term microscopy-based taxonomic datasets.

Conclusion
This study advanced our understanding of how in situ CM abundance and proportion are influenced by not onl y envir onmental factors, but by the taxa present.It is critical to identify active CMs, the factors favoring their abundance and proportion as active CMs, and their different responses to changes in biotic and abiotic factors (Stoecker 1998, Edw ar ds 2019 ).It w as demonstr ated that measur ements of potential CM abundance sho w ed a substantial ov er estimation of the presence of CMs when compared to active CM abundance; the presence of mixotrophic taxa in the system does not necessarily mean mixotrophic activity.This highlights the limitations within studies that recategorize phytoplankton from microscopy-based taxonomic datasets to estimate mixotroph abundance, as this study shows that species were not al ways gr azing when they wer e pr esent.Anal yses with historical datasets are still very useful because they can help rapidly increase understanding of large-scale mixotroph presence and distribution, but it is necessary to better constrain the organisms classified as mixotrophs based on the conditions they are grazing under.The BrdU method can be used to help identify the conditions under which CMs are actively grazing, and more accurately use those datasets .T his wa y, our knowledge on the biogeogr a phical distribution of CM abundance can be expanded ov er lar ge time periods to better understand their contribution to aquatic food webs.

Figure 1 .
Figure 1.Map of location of fieldwork, with the two different sampling stations across the York River, VA marked with blue dot.WP Station-WP Fishing Pier and GP Station-GP Fishing Pier.
through an immunoprecipitation process to recover DNA that had incor por ated BrdU.Nucleic acids r ecov er ed fr om the −BrdU incubation were not immunoprecipitated but used directly for amplicon PCR.The first step in the imm unopr ecipitation pr ocess involv ed pr e paring an antibod y mix and a magnetic bead mix that wer e bloc ked with denatur ed bacterial DNA to help pr e v ent nonspecific binding.A total of 300 ng of P. angustum DNA per sample (10 ul per reaction) was denatured for 10 min at 95 • C, the tubes wer e tr ansferr ed to an ice bath for 2 min, and then cold PBS-BSA [1 mg acetylated bovine serum albumin (BSA) per ml 1x PBS] was added to bring to a final volume of 30 μl per reaction.The denatured bacterial DN A w as used to block the antibodies (anti-BrdU B44; BD Biosciences 347 580; 1/10 dilution: 2.5 ng μl −1 ; 10 μl each for each reaction) and the beads (Dynabeads M-280 Sheep antimouse magnetic beads; Invitrogen 11201D; 10 μl beads + 10 μl of PBS-BSA).The denatured bacterial DN A w as mixed with same volume of antibody or beads in a 1.5-ml microfuge tube and incubated at 4 • C on a rotating mixer for 1 h (antibody mix) or overnight (bead mix).

Figure 2 .
Figure 2. Total abundance of phototrophs ( ±SE) at WP and GP betw een Mar ch 2021 and February 2022.The x -axis r epr esents a year-long timeline.Each pair of bars (WP and GP) r epr esents a single sampling date along that timeline.

Figure 3 .
Figure 3. Proportion of total abundance that six plankton groups, identified through microscopy, accounted for at each sample date between March 2021 and February 2022 at (A) WP and (B) GP.The x -axis r epr esents a year-long timeline.Each bar represents a sampling date.

Figure 4 .
Figure 4. Proportion of mixotrophic ASVs of each major taxa group at (A) WP and (B) GP stations.

Figure 5 .
Figure 5. Number of active CMs for each major taxa group identified through Illumina sequencing at WP (A) and GP (B) stations for each sampling date.Each bar corresponds to the number of individual ASVs classified as being a CM based on c hlor oplast containing sequences that were ingesting BrdU-labeled bacteria.The x -axis r epr esents a year-long timeline.Each bar represents a sampling date.

Figure 6 .
Figure 6.Number of genera/species of each major taxonomic group classified as an active CM in the microscopy samples at (A) WP and (B) GP stations.The x -axis r epr esents a year-long timeline.Each bar represents a sampling date.

Figure 7 .
Figure 7. (A) Estimated abundance of potential CMs in all microscopy samples based on the genera of chloroplast-containing plankton that have demonstrated the of ingesting bacteria in previous experiments.(B) Total abundance of active CMs in microscopy samples identified through the BrdU experiments for each sampling date at WP and GP stations .T he active CMs were identified based on the genera/species of chloroplast-containing plankton that were ingesting BrdU-labeled bacteria for each sampling date.(C) Estimated proportion of potential CMs in all microscopy samples based on the genera of chloroplast-containing plankton that have demonstrated the of ingesting bacteria in previous experiments.(D) The proportion of active CMs present at WP and GP stations .T he proportional abundance of active CMs was calculated by dividing the abundance of genera/species identified as active CMs by the total abundance of phototrophs of each sampling date .T he x -axis r epr esents a year-long timeline.Each pair of bars (WP and GP) r epr esents a single sampling date along that timeline.

Table 1 .
Av er a ge ( ±SE) water quality data collected 24 times between March 2021 and February 2022 at two stations in the York River.

Table 2 .
Results for the quasipoisson and quasibinomial GLM with an offset analysis for select groups at WP and GP.