Relative importance of bacterivorous mixotrophs in an estuary‐coast environment

Mixotrophic eukaryotes are important bacterivores in oligotrophic open oceans, but their significance as grazers in more nutrient‐rich waters is less clear. Here, we investigated the bacterivory partition between mixotrophs and heterotrophs in a productive, estuary‐influenced coastal region in the East China Sea. We found ubiquitous, actively feeding phytoplankton populations and taxa with mixotrophic potential by identifying ingestion of fluorescent prey surrogate and analyzing community 18S rRNA gene amplicons. Potential and active mixotrophs accounted for 10–63% of the total eukaryotic community and 17–69% of bacterivores observed, respectively, contributing 6–48% of estimated in situ bacterivory. The much higher mixotroph fitness outside of the turbid plume was potentially driven by increased light and decreased nutrient availability. Our results suggest that, although heterotrophs dominated overall in situ bacterivory, mixotrophs were abundant and important bacterivores in this low‐latitude mesotrophic coastal region.

Heterotrophic flagellates and microzooplankton were traditionally considered to be the primary, if not sole grazers of bacteria in aquatic food webs, but now there has been increasing recognition of the importance of unicellular phagotrophic mixotrophs (i.e., planktonic eukaryotes that are capable of photosynthesis and phagocytosis; Sherr and Sherr 2002;Stoecker et al. 2017).Model simulations suggest that mixotrophy can substantially increase primary production and energy transfer efficiency (Stoecker et al. 2009;Mitra et al. 2014;Ward and Follows 2016); meanwhile, bacterivory by mixotrophs may introduce different impacts on elemental flux and biogeochemical cycling.For instance, heterotrophs seem to respire more of the organic carbon and other elements from ingested prey to acquire energy (Hansen et al. 1997), whereas mixotrophs may utilize prey biomass primarily for nutrients, subsidizing/facilitating energy and organic carbon production through photosynthesis (Caron et al. 1993;Ptacnik et al. 2016;Mitra and Flynn 2023).
The relative prevalence and success of mixotrophs have been commonly associated with oligotrophic environments (Unrein et al. 2007;Hartmann et al. 2012;Duhamel et al. 2019;Li et al. 2022), suggesting that low nutrients may select for mixotrophy (Lindehoff et al. 2010;Edwards 2019).However, mixotrophs and mixotrophic bacterivory have also been reported to be ubiquitous in systems with relatively high nutrient supply, such as equatorial East Pacific (high light; Stukel et al. 2011) and temperate estuaries adjacent to the North Atlantic (low light; Millette et al. 2017Millette et al. , 2021)).Other studies have found higher irradiance to benefit mixotrophs, presumably because their phagotrophic abilities come at the expense of phototrophic performance (Raven 1997;Tittel et al. 2003;Edwards et al. 2022).These varied results, coupled with the observed functional diversity (Li et al. 2022) and metabolic plasticity among mixotrophs (Li et al. 2021), illustrate that there are multiple interacting variables that determine whether mixotrophs will have an ecological advantage in a particular habitat.Understanding the range of habitats and mechanisms that allow mixotrophs to thrive is, therefore, critical for interpreting their ecosystem-level biogeochemical consequences, yet remains challenging (Leles et al. 2021).
In this study, we ask what are the major mixotrophic populations, how fast they can graze, and what are the factors possibly affect bacterivory in east China sea coast receiving freshwater discharge from the Yangtze River.Hypotheses were proposed to observe a shift between mixotrophic and heterotrophic grazer communities and bacterivory contributions along that estuary-coast environmental gradient, corresponding to the changing inorganic nutrients and other biophysical conditions.To answer these questions and test our hypotheses, we developed a robust method to estimate active mixotrophs in situ via tracer experiments with fluorescent prey surrogates coupled with epifluorescence microscopic (EM) observation and high throughput flow cytometry (FC) screening.We also used 18S rRNA gene metabarcoding to evaluate overall eukaryotic taxonomy and identify potential mixotrophic bacterivores.Collectively, we aimed to investigate whether a representative coastal region influenced by river discharge favors mixotrophs.

Study sites and biochemical parameters
Water samples were collected from 8 to 15 November 2021 on board the "Zhe Yu Ke" research vessel from one estuary station and three shelf stations in the coastal East China Sea.Surface water from a 2-m depth was collected by rosette mounted on an SBE19plus CTD equipped with auxiliary sensors.Inorganic nutrients, including DIN (NO 3 ) and DIP (PO 4 ), were measured in the lab using the Auto Discrete Chemical Analyzer (Smartchem 600; AMS Alliance).Total chlorophyll a (Chl a) was filtered onto 0.7 μm Whatman GF/F glass fiber filters, extracted with acetone overnight in the dark, and measured on a Turner Designs AU-10 fluorometer.Cell abundance samples for eukaryotes and bacteria were fixed with glutaraldehyde at 0.5% final concentrations (Sigma-Aldrich) and stored at À80 C until analyzed via FC in the lab (CytoFLEX S; Beckman Coulter).

Short-term grazing experiments
Three prey surrogates were first assessed before the cruise to identify the most suitable tracers, among which the 1.0 μm green-yellow fluorescent beads produced the best detection results by FC (Supporting Information Fig. S1) and were thus chosen for the grazing experiments.During the cruise, approximately 200 mL pre-filtered (200 μm nylon mesh) seawater in triplicate was amended with 1.0 μm beads at a final concentration of $ 2.5 Â 10 5 mL À1 , accounting for 20-50% of the in situ bacteria concentrations.Samples were incubated for 1 h under 14-18 C, $ 70 μE s À1 m À2 , and 12 h : 12 h light/dark cycle.Twenty-milliliter fixed (0.5% glutaraldehyde) and stained aliquots (0.1 μg mL À1 DAPI; Biotium) were filtered onto 2 μm polycarbonate membrane filters (Millipore) for microscopic examination (Olympus BX61 inverted microscopy).

Active bacterivores identification
Under the microscope, cells with ingested beads were counted as active bacterivores, and mixotrophs were distinguished based on the presence of conspicuous chloroplasts (organelles with red autofluorescence).Potential non-constitutive mixotrophs with keptoplastic plastids were treated as heterotrophs as we could not verify their mixotrophic activities.Each sample was scanned for 130-150 fields to decrease counting errors (Lund et al. 1958).Active bacterivores were also identified on FC (CytoFLEX S; Beckman Coulter).Total pigmented eukaryotes and heterotrophic flagellates were first identified via red fluorescence (Chl a), blue fluorescence (stained nucleic acid), and size scatters, after which active bacterivores were distinguished as the same eukaryotic population, but which had also acquired the characteristic yellow-green fluorescence of ingested beads (γ em $ 505 nm).

Grazing rates estimate
The average per capita ingestion rate of beads (IR i , beads cell À1 h À1 ) in a sample was calculated using the following equation: where B i,t and B i,0 , are the total number of ingested beads observed at time t (i.e., 1 h post-incubation) and the initial time point of group i (bacterivores were binned into six size groups).N i,t are the total active bacterivores observed at time t of group i.
The estimated in situ bacterivory rate (BR), that is, the grazing impact of bacterivores in aggregate (bac.mL À1 d À1 ), was calculated using the following equation: where C i , C bac , and C b are the concentrations (cells mL À1 ) of active bacterivores (group i), total bacteria, and fluorescent beads (initial concentrations), respectively.C i was calculated by multiplying mean cell counts per field under microscopy by 3.9 Â 10 4 fields (obtained by dividing the membrane filters' area by each observational field's area) and dividing by aliquot volume (mL).

18S rRNA gene metabarcoding
We vacuum filtered approximately 1 L seawater from each sample through a 47 mm, 0.45 μm polycarbonate membrane filter (Millipore) and extracted total environmental DNA, which was used as a template to amplify the hypervariable V4 region of 18S rRNA with primers 547F and V4R (Stoeck et al. 2010).PCR products from each sample were purified and pooled in equal amounts, followed by pair-end 2 Â 250 bp sequencing using the Illumina NovaSeq platform at Shanghai Personal Biotechnology Co., Ltd.Quality check and chimeras removal of amplicon sequence variants (ASVs) were conducted using DADA2 (Callahan et al. 2016), which were further classified with two databases, the Protist Ribosomal Reference Database (Guillou et al. 2013) and SILVA (Pruesse et al. 2007) in QIIME2 (Bolyen et al. 2019).Classification disagreement between the two databases were manually curated by searching additional cultivated references using BLASTn and the NCBI GenBank (Supporting Information Table S1).ASVs that were metazoa, freshwater/land-originated (mostly fungi and some amoeba at B3), or with low numbers were excluded from further analysis.

Correction factors for 18S rRNA gene abundances
Ribosomal gene copy numbers vary substantially among eukaryotic groups, from as low as one copy in small flagellates to thousands of copies per cell in dinoflagellates (Zhu et al. 2005).This, along with biases in the amplification reaction (Gonzalez et al. 2012;Bradley et al. 2016), can result in relative gene abundances in amplicon pools that are significantly different from cell abundances (Karlusich et al. 2022).To correct for this bias, we followed the principle of a study by Martin et al. (2022) and developed our own empirical gene-to-cell correction factors (C.F.) for seven major groups that were distinguishable by morphology under microscopy (dinoflagellates, ciliates, radiolaria, cryptophytes, diatoms, other pigmented, and heterotrophic eukaryotes), which were formulated as: where %MB_gene g and %EM_cell g are the relative abundances of group g derived from metabarcoding 18S rRNA reads and from epifluorescence microscopy (EM)-based cell abundances, respectively.These C.F. were then averaged with those of Martin et al. (2022) to retrieve the final factors used in this study (results shown in Supporting Information Fig. S2).The corrected relative abundance of group g (%MB_cell g ) can be subsequently calculated using Eq. 4 and were used in all of our analysis except for Fig. 1B.

Annotation of potential mixotrophs
We assigned each genus (or species when the annotation was available) to one of four nutritional categories, that is, mixotrophy (including constitutive and non-constitutive mixotrophs), heterotrophy, autotrophy, or others, using evidence from synthesized databases (e.g., Faure et al. 2019;Mitra et al. 2023) and research articles (Supporting Information Table S2).Note that one of the abundant genus Micromonas of Chlorophyta, was assigned as autotrophs as evidence for phagotrophy was ambiguous (McKie-Krisberg and Sanders 2014; Jimenez et al. 2021).Some Radiolaria and Noctiluca were considered as others/symbiotic, all Syndiniales were assigned as others/parasitic and unknown taxa or taxa with unavailable trophic evidence were all assigned as "others/unknown."

Statistical analysis
All of our statistical analysis and visualization were carried out in R, v3.6.3 (R Core Team 2020).Pairwise t-tests were performed using basic R functions.Redundancy analysis (RDA) and Spearman correlation matrix were applied to investigate significant correlations between mixotrophic fitness and biophysiochemical variables, using the vegan and psych packages for calculation and the pheatmap package for visualization.ASV reads were first Hellinger-transformed, and all independent variables were log-transformed for normalization.Low variance inflation factor (< 10) was checked to eliminate co-linearity in the RDA.

Hydrographic and biochemical characteristics of the study region
Geographically, the estuary Sta.B3 within the plume region showed the lowest salinity (14.2) and temperature (14.2 C) and the highest nutrients (62.1 μmol L À1 NO 3 , Fig. 1.Oceanographic map with salinity in the background showing sampling sites of this study (A), 18S rRNA genes derived taxonomic composition of top 12 divisions before (B), and after gene-to-cells conversion using correction factors (C), and a heatmap constructed with Pearson distances showing trophic mode compositions in relative abundances (RA) using results in panel C but on genus level (D).Taxonomy in panels b and c were division levels based on the PR 2 database ranks, with the exception that Ochrophyta was presented in three classes (Bacillariophyta, Chrysophyceae, and Raphidophyceae).Note that unknown reads were not corrected for gene copies cell À1 , which may result in an overestimated or underestimated percentage.

Metabarcoding-derived community composition and potential mixotrophs
Before the gene-to-cell correction (in %MB_gene), amplicon data suggested that all stations were dominated by 18S rRNA genes from Dinoflagellata, except for B3, which was co-dominated by Bacillariophyta and Dinoflagellata.Other abundant groups included Radiolaria (at B13 and D7) and Cryptophyta (at D3 and D7; Fig. 1B).After correction for gene copy numbers (in %MB_cell) it is estimated that the three most abundant populations at B3 were Bacillariophyta, Chlorophyta, and Cercozoa (Fig. 1C), among which the majority genus were annotated as autotrophs and heterotrophs (Fig. 1D; Supporting Information Table S2).Sta.B13 was dominated by pigmented eukaryotes from Dinoflagellata and Haptophyta and heterotrophs from the MAST group.Cryptophyta and Chlorophyta were both highly abundant at D3 and D7, followed by Dinoflagellata and Haptophyta (Fig. 1C).Among all pigmented eukaryotes from Stas.B13, D3, and D7, the majority of most abundant genera were potential mixotrophs, such as Cryptomonadales and Teleaulax of Cryptophyta, Pyramimonas of Chlorophyta, Chrysochromulina and Phaeocystis of Haptophyta, as well as Gymnodinium and Heterocapsa of Dinoflagellata (Fig. 1D; Supporting Information Table S2).When considering both taxonomy and nutritional mode composition (in %MB_cell), the Pearson distance clustering suggested a closer association between B3 and B13 and D3 and D7, respectively (Fig. 1D).Particularly, the B transect had substantially more heterotrophs than the D transect (24-25% vs. 5-7%), while B3 presented the lowest (potential) mixotrophic proportions among all four station.

Microscopy and FC-identified active mixotrophs
Ubiquitous, actively feeding mixotrophs and heterotrophs were identified from small (< 5 μm), medium (5-10 μm), and large-sized (> 10 μm) eukaryotes under EM (Fig. 2A).Among all six size categories, 5-10 μm heterotrophs and mixotrophs were the most abundant bacterivores, contributing a mean of 71.6% and 28.4% of total bacterivory across all stations (Supporting Information *PE and HF abundance were the mean values averaged from flow cytometry and epifluorescent microscopy observations.Values in brackets are standard deviations between the two methods. Fig. S3).Averaged ingestion rates of all mixotrophs and heterotrophs varied from 1.5 to 3.8 beads cell À1 h À1 , accompanied with estimated in situ abundances and BRs of 0.3-1.8Â 10 3 cells mL À1 and 0.2-3.2Â 10 5 bact.mL À1 d À1 , respectively (Table 2).
For the estimated BRs, mixotrophs from all size classes collectively contributed significantly higher proportions outside of the plume (20% at B13 and 41-48% at D3/D7) than B3 (6%; pairwise Student's t-test; p < 0.05; Supporting Information Fig. S5a).The size category of dominant grazers did not vary across stations, that is, all from the 5-10 μm medium size range, suggesting the prey (size) is likely controlling the grazer size (Supporting Information Fig. S5b,c).Comparison of other physical and biochemical conditions across stations revealed the highest similarity between D3 and D7, followed by B13 and then B3 (Supporting Information Fig. S6a).Pearson correlation matrix suggested that bacterivory by mixotrophs was positively correlated with heterotrophic bacteria, Synechoccoccus, and abundances of potential mixotrophs (p < 0.05; Supporting Information Fig. S6b).

Discussion
Methodology to identify active mixotrophs EM and FC approaches were demonstrated to be robust for retrieving active bacterivores from natural environments, but both have limits.Microscopic observation could provide direct phagocytosis evidence but was time-consuming.FC screening was high throughput but failed to estimate grazing rates.One challenge of using the EM method is to recognize the position of prey (inside or outside of the cell).Because many mixotrophic flagellates lack cell walls, that is, shown as cells without outlines/borders under fluorescence excitation, sometimes prey appear to be, but may not be, attaching to (outside of) the cells (Fig. 2A).We have noticed the same issue during one of our previous studies using mixotrophic cultures (Florenciella sp. of Dictyochophyceae), and we partially solved the problem by staining cell membranes (CellBright 488; Biotum) to confirm that oftentimes prey was indeed ingested or in the process of entering the membrane (Li et al. 2021).In this study, we examined additional liquid samples under a bright field (have clearer cell outline) to confirm phagocytosis in common mixotrophs seen on membrane filters (Fig. 2A).Nevertheless, ingestion rates derived from EM should be annotated with caution.Other studies have also indicated that 1.0 μm beads tend to get stuck in the cell membrane during the process of phagocytosis by their mixotrophic grazers (Millette et al. 2021).
FC method generated higher estimates of mixotrophs than microscopy for our samples, agreeing to a study using a similar method of FC and acidotropic probes by Costa et al. (2022).It was either due to underestimated EM results and/or overestimated FC results, and if the latter, these overestimates were likely from omnivorous heterotrophs that ingested both pigmented eukaryotes and bacterial prey surrogates (some of the lowest panel images in Fig. 2A).Combined methods using various tools with a focus on single cells (such as single-cell sorting and sequencing) were suggested to retrieve the most reliable data for future mixotrophy research (especially for taxonomic identification; Beisner et al. 2019).

Trophic mode assignment
Annotation of putative trophic modes on the basis of taxonomy is another effective approach to study mixotrophy but can sometimes lead to biased interpretations.For example, some studies have found interspecific and intraspecific variations of trophic mode from ciliates Mesodinium sp.(Johnson et al. 2016) and dinoflagellate Karlodinium sp.(Calbet et al. 2011).Although others have suggested mixotrophy tends to broadly exist among different species within the same genus, such as Ochromonas of Chrysophyceae (Wilken et al. 2019), Florenciella of Dictyochophyceae (Li et al. 2022), and Phaeocystis of Haptophyta (Koppelle et al. 2022).Future experiments need to be carried out to verify whether these trophic variations are intrinsic (evolutionary) or due to plasticity under different physiological and environmental conditions.

Enhanced mixotrophic fitness outside of the turbid estuary
Mixotrophs compete for inorganic nutrients with autotrophs and organic nutrients (prey) with their heterotrophic competitors.The relative fitness and niches partition among them under different habitats remains an important question to answer.We chose an environment that covers light-limiting eutrophic estuary and coastal sites with higher irradiance and less nutrients.Our results suggested that mixotrophs were suppressed at the highly turbid estuary, where strict autotrophs and heterotrophs were superior competitors, but became dominant outside of the estuary plume.These niche shifts are consistent with many studies that have suggested mixotrophs benefit from low-nutrients and high-irradiance environments (Tittel et al. 2003;Stukel et al. 2011;Anderson and Hansen 2019;Edwards 2019;Edwards et al. 2022).One explanation for this is that mixotrophs are penalized for possessing dual apparatus of phototrophy and phagotrophy within the same cells, that is, they have poorer photosynthetic performance compared to autotrophs and lower ingestion rates compared to heterotrophs (Raven 1997).With enhanced irradiance, decreased nutrients, and abundant bacterial prey (e.g., sites outside of the estuary), mixotrophs benefited from both inorganic and organic nutrients, which can offset nutrient limitations and stimulate more growth of mixotrophs.Although mixotrophs may have lower ingestion rates than their heterotrophic competitors, phototrophy-derived carbon and energy (when light is unlimited) can increase the competitive ability of mixotrophs against heterotrophs (tend to be energy-limited), meanwhile suppressing prey to an equal or lower concentration and making them comparable or superior competitors (Rothhaupt 1996;Fischer et al. 2016;Edwards et al. 2023).
Overall, this study demonstrated an ecological niche for phagotrophic mixotrophs in a low-latitude mesotrophic coastal region impacted by land riverine inputs, expanding our knowledge of mixotrophic distribution and niche fitness.Finally, the rapid shifts in the partitioning of bacterivory between mixotrophs and heterotrophs happened at small geographic scales (e.g., less than 100 km), reiterating the importance of high-frequency and high-resolution sampling efforts for future research.

Fig. 2 .
Fig. 2. Epifluorescence microscopy (EM) images showing evidence of ingestion (beads in green-yellow or cyan-blue fluorescence) by representative mixotrophs (chloroplasts in red or orange fluorescence) and heterotrophs among different size groups (A).Flow cytometric plot (B) from Sta. B13 shows pigmented eukaryotes (PE) and active mixotrophs (AMF) in green dots and heterotrophic flagellates (HF) and active heterotrophs (AHF) in red dots, initially distinguished by red fluorescence (chlorophyll a) and blue fluorescence (DAPI-stained nucleic acid) as shown in Supporting Information Fig. S4.Panel (C) compares results between EM and flow cytometry (FC) at each station, including relative abundances of AMF against PE, AHF against HF, as well as %AMF and %AHF against total bacterivores.Panel (D) demonstrates the same results in panel C, but with averaged values across all stations.All scale bars in panel A are 5 μm in length.The first three images in the lower mixotroph panel demonstrate the beads' position in liquid samples (similar morphology), and the fourth and fifth images show the same cell (as in the upper panel) excited with blue (DAPI) and orange (phycobilins) fluorescence.

Fig. 3 .
Fig. 3. Composition of four nutritional modes annotated by 18S rRNA amplicon (a) and bacterivory partition between active mixotrophic flagellates (AMF) and heterotrophic flagellates (AHF) across stations (b).Log-transformed turbidity (NTU) in grey dots, and nitrate concentrations (μM) in yellow dots (both correspond to the second axes on the right) were overlay on top of bars in panel (a) and (b), respectively.Panel (c) shows the redundancy analysis (RDA) between the top 10 genera (summed across all stations) from each nutritional mode, and selected independent variables (i.e.PMF for potential mixotrophs, PE for pigmented eukaryotes and Turb.for turbidity).Panel (d) is a turbidity map with representative bacterivore illustrations to demonstrate community shift among mixotrophs (green and brown colored flagellates), autotrophs (chained diatoms) and heterotrophs (red colored flagellates).Taxa names with the "X" suffix in panel (c) mean that they could not be annotated on lower taxonomy ranks.

Table 1 .
Physical and biochemical characteristics of the study sites.

Table 2 .
Summary of grazing rates and abundances of mixotrophic and heterotrophic bacterivores in the investigated region.