For microzooplankton, rates of prey ingestion are influenced by the time involved in both prey capture and prey digestion. The former includes encounter rates, prey handling and prey avoidance and/or rejection, the latter depends on how quickly prey can be processed through the digestive system before the predator becomes satiated (Flynn and Mitra, 2016). Addition of surrogate prey may affect both aspects and, even if selected against by the predator, the additional handling time involved will reduce the ingestion rate of natural prey and hence their rates of grazing mortality.

Demonstrations of the functional response of protozoan and metazoan zooplankton involve measurement of ingestion rates of grazers at increasing concentrations of a single prey type, including in some cases, artificial prey such as latex beads (Fenchel, 1980; Jonsson, 1986; Verity, 1991). In a procedure that mimics determination of the functional response, the proposed saturation approach involves adding increasing abundances of surrogate prey to incubations of natural communities and measuring the observed growth response of the natural prey of interest. At increasing levels of surrogate prey addition, as grazers become “saturated” by handling and ingesting surrogate prey, the apparent growth rate (μ_app) of the natural prey increases and approaches a maximum value that represents their specific growth rate in the absence of grazing (μ_max). The difference between μ_app with no addition of surrogate prey, effectively the net growth rate (μ_net), and μ_max, provides an estimate of the level of mortality due to grazing (m)

with the assumption that the specific growth rate of the prey is not affected by the addition of surrogate prey.

Analagous to plant growth in response to increasing amounts of a limiting nutrient, the increasing growth rate of prey in response to increasing surrogate prey was represented by an asymptotic regression model, of similar form to the Mitscherlich response model (Mitscherlich, 1928). The model was used to derive values of μ_net, μ_max and m from the relationship between μ_app and level of addition of surrogate prey,

where C is the surrogate prey concentration and a is proportional to the relative rate of increase in μ_app with increasing C. Note, u_max is referred to as “gross growth rate” but does not account for other causes of mortality, such as viral infection. Model fitting was carried out using the non-linear least squares regression function (nls) in R (R Core Team, 2021).

Proof-of-concept tests were performed using laboratory cultures of a single isolate representative of the picophytoplankton and a typical protistan predator. Specifically, the experiments used the mixotrophic chrysophyte Ochromonas danica (CCMP 1391), ∼ 6–9 μm in diameter, as a predator, and the generally abundant picoeukaryote Micromonas pusilla (CCMP 485), ∼ 1–2 μm in diameter, as prey. The cultures were maintained in sterile L1 media in 50 ml borosilicate glass test tubes in a 21^°C incubator with a 14 h light/10 h dark cycle and light levels of ∼ 90 μE m^–2 s^–1. Prior to the start of the experiments, O. danica was transitioned from K media and rice to only sterile L1 media and then fed every other day on M. pusilla. The prey, M. pusilla was maintained in semi-continuous growth through regular transfers (2–3 days). In 3 experiments, similar ratios of M. pusilla (target 3 × 10⁵ ml^–1) and O. danica (target 2,000 ml^–1) were generated in tubes containing a total volume of 40 ml. Fluorescent polystyrene microspheres (beads) of 2 μm in diameter (Fluoresbrite Plain YG microspheres, Polysciences, Inc., Warrington, PA), were used as surrogate prey.

To minimize clumping, the beads were blocked in a solution of 1% bovine serum albumin overnight then centrifuged for 5 min at 2,000 rpm, after which the pellet was resuspended in 0.2 μm filtered seawater. A new solution of beads was prepared at the start of each experiment. Duplicate tubes were prepared for 6 different surrogate prey concentrations ranging from 0 to 7.5 × 10⁶ ml^–1. During the incubations, the experimental tubes were rotated at 1.2 rpm on a plankton wheel to keep particles in suspension. Two subsamples of 1 ml were removed from each tube at the start (T₀) and the final time point (T_F) after ∼24 h, for flow cytometric analysis.

The experimental format established in the culture experiments was then adapted for use in natural waters. From 16th July to 15th August, 2019 (days of the year 197–227) and 8th June to 26th July 2021 (days 165–207), a series of saturation experiments were performed that focused on determining the growth rates and grazing mortality of picoeukaryote and Syncechococcus sp. populations. Gulf of Maine coastal seawater was collected from the Damariscotta River Estuary at Bigelow Laboratory’s dock off East Boothbay, Maine, United States (43.8604° N, 69.5781° W). Water for all experiments was collected within 1 h of high tide at 1 m depth using a 5 l Niskin bottle and was gravity filtered through 200 μm mesh to remove zooplankton. Several Niskin casts of seawater were combined in an acid washed carboy before being siphoned into 600 ml acid washed polycarbonate bottles. For each saturation experiment, 12–14 bottles were used to generate a range of levels of surrogate prey addition. Fluorescent polystyrene microspheres of 2 μm in diameter, treated as for the laboratory experiments, were added to the bottles either as duplicate treatments or in a continuous series of abundance. Maximum surrogate prey abundance ranged from 1.7 × 10⁶ beads ml⁻¹ on day 197 to 3.8 × 10⁶ beads ml⁻¹ on day 226, in 2019.

Bottles were incubated for 24 h in a flow-through incubator on the dock under a nylon mesh that removed ∼40% of the surface PAR. In an initial test of the flow-through incubator, five 600 ml bottles were filled with seawater from the same carboy and with no further treatment, were incubated for 24 h. The counts of the T₀ abundance of picoeukaryotes from the five bottles showed a coefficient of variation of 1.3%, that increased in the T₂₄ counts to only 4.1%, indicating consistent growth among the 5 replicate bottles. Following all experiments, the seawater to which beads had been added was filtered through a 0.45 μm capsule filter to recover the beads before discarding the water.

For comparison, a series of dilution experiments was carried out in parallel to the saturation experiments, using a modification of the approach of Landry and Hassett (1982). The same water used for the saturation experiments was used to prepare triplicate bottles of the whole water (100%) and duplicate bottles at four different levels of dilution (10, 25, 50, and 75%) in 600 ml polycarbonate bottles. Diluent was prepared from the 200 μm pre-screened seawater by gentle filtration through a 0.45 μm mini-capsule filter (Pall Life Sciences). The dilution experiments were set up following the saturation incubations to allow time for flow cytometric analysis of each set of subsamples and so, T₀ and T₂₄ were delayed by ∼ 2 h relative to those of the saturation experiments.

A ZE5 Cell Analyzer (Bio-Rad, Hercules, CA, United States) was used to measure optical properties of single cells from each sample. Particles were excited with a 488 nm blue excitation laser (100 mW). Data acquisition was triggered on forward scatter (FSC). Signals were recorded from detectors with bandpass filters for right angle light scatter and fluorescence emission in red (692 nm / 80 nm band pass) indicative of chlorophyll a, orange for phycoerythrin (593/52 nm), and green (525/35 nm). To ensure accurate calibration of the flow cytometer, ZE5 QC beads (Bio-Rad, Hercules, CA, United States) were run daily. Data files were analyzed from logarithmic dot plots based on fluorescence and characteristic light scattering properties (Durand and Olson, 1996) using FlowJo 10.6 Software (Becton Dickinson & Company, San Jose, CA, United States) (Supplementary Figure 1).

For culture and natural community experiments, surrogate prey were gated using FSC and green fluorescence. Despite the blocking procedure, some clumping of beads did occur but, in this application, doublet and triplet clumps were treated as single surrogate prey particles. For laboratory experiments, O. danica was gated on side scatter and red fluorescence and the M. pusilla populations were determined using forward scatter and red fluorescence. For field experiments, picophytoplankton were identified based upon cell size, using forward scatter as a proxy, and red fluorescence. Picophytoplankton were classified, as either Synechococcus or picoeukaryotes based on the presence or absence of the accessory pigment phycoerythrin, respectively. Using these gating criteria, the number of cells in each identified population was enumerated and converted to cell abundances based on the processed sample volume.

The single prey and predator laboratory culture experiments confirm the theoretical prediction. Apparent growth rate (μ_app) of M. pusilla increased with the addition of beads and tended towards a horizontal asymptote at increasing levels of surrogate prey addition, indicating saturation of ingestion in O. danica (Figure 1). Fits of Equation 2 to the relationship between μ_app and the surrogate prey abundance, provided highly significant estimates of μ_max (p-value < 0.001), μ_net (p < 0.05) and a (p < 0.05) in each of the experiments (Table 1). Derived values of m decreased from 0.60 to 0.33 d^–1 over the three experiments, despite the reduced prey to predator ratios. Although a distinct population of O. danica could be counted initially, as they grazed on M. pusilla and ingested beads, the population became more difficult to define in the flow cytometric signal and it was not possible to obtain final abundances and reliable growth rates for the predator.

Initial prey to predator ratios were approximately 200, 180, and 160 in experiments A, B, and C, respectively (Table 1). At the levels of abundance used, growth and grazing rates were almost balanced in the cultures with no bead addition, resulting in similar, relatively low rates of μ_net for M. pusilla in the three experiments (Table 1).

In general, the abundances of freely suspended beads decreased during the incubations by an average of 25, 3, and 25% in A, B, and C, respectively. The decrease possibly results from a combination of ingestion by predators, sinking out of suspension despite the mixing during incubations and sticking to other particles. Proportionally higher losses tended to occur at higher levels of bead addition but this was not always the case.

The dynamics of the two picophytoplankton populations showed variable patterns indicative of differences in their rates of growth and the grazing pressure they experienced. During the 2019 summer period, the abundance of the picoeukaryote population increased by approximately 50% from 43,100 cells ml^–1 on day 197 to 60,600 cells ml^–1 on day 211 and then gradually declined over the following 2 weeks (Figure 2A). In contrast, the abundance of Synechococcus increased more than 10-fold from 2,860 on day 197 to 30,100 on day 211 and remained at a relatively high level to day 226 (Figure 2A). Similar abundances of each group were observed in the summer of 2021 (Supplementary Tables 1, 2). Note, the mouth of the Damariscotta River Estuary experiences high levels of tidally driven water exchange and this may also have influenced population dynamics.

As illustrated for 2019, increased bead addition resulted in increased μ_app, tending toward a horizontal asymptote at bead abundances that exceeded ∼2.0 × 10⁶ ml^–1 (Figures 2B–H). The growth response to bead addition indicates increasing saturation of the existing grazer population and reduction of grazing mortality in the picophytoplankton populations. Fits of Equation 2 to the relationship between μ_app and the surrogate prey abundance, provided highly significant estimates of μ_max (p < 0.001), μ_net (p < 0.01) and a (p < 0.05), for each of the taxonomic groups in all seven experiments, except for Synechococcus on day 206 (Table 2). Values of m ranged from 0.26 to 0.43 d^–1 for picoeukaryotes and 0.23 to 0.42 d^–1 for Synechococcus. Estimates of gross growth in the absence of grazing (μ_max) indicate a highly productive population of picoeukaryotes that grew at rates of over 1 d^–1. These rates of μ_max gradually declined over the period to close to 1 doubling d^–1, by day 227 (Table 2). Synechococcus also showed high rates of μ_max initially, that stabilized to between 0.47 and 0.53 d^–1 from days 212 to 227. Similar rates of growth and mortality were measured in the four experiments successfully conducted at the same location during the summer of 2021 (Supplementary Table 2).

Paired t-tests showed significant differences in μ_net (p = 3.5 × 10^–5), m (p = 1.3 × 10^–2) and μ_max (p = 1.9 × 10^–4) occurred between Synechococcus and picoeukaryotes in the combined dataset of saturation experiments. Rates of grazing mortality were the most comparable of the coefficients between the two taxonomic groups, while higher rates of gross growth among picoeukaryotes were largely responsible for the significant differences in net growth rates (Supplementary Figure 2).

In general, the rates from each approach show a similar range of values for μ_max, μ_net and m (Figure 3 and Supplementary Table 1). Successful, parallel experiments allow direct comparison on days 197 and 204 in 2019 and on days 174, 182, 196, and 207 in 2021. It should be noted that the dilution incubations were started approximately 2 h after the saturation incubations, which may explain the slight differences in μ_net between experiments on the same days (Table 2 and Supplementary Tables 1, 2). In the combined dataset of comparable experiments, there was no statistical evidence that the mean difference between experimental approaches was significantly different from zero (paired t-tests, p-value < 0.05). This is apparent in the lack of clear trends between approaches when estimates of μ_max and m are compared for picoeukaryotes or Synechococcus, although large variations were apparent on certain days, particularly days 197 and 182 (Figure 4).

Note, the raw data from both the saturation and dilution experiments is included in the Supplementary Material.

The saturation approach confirmed intriguing differences in the dynamics of co-occurring picoeukaryote and Synechococcus populations. The most notable difference being the considerably higher rates of u_max in the picoeukaryote population, equivalent to between ∼1 and 2 doublings d^–1 compared to generally < 1 doubling d^–1 for Synechocccus (Table 2, Figures 3, 4, and Supplementary Figure 2). Similar, higher rates of gross growth and grazing mortality of picoeukaryotes compared to Synechococcus, have been observed in other temperate and sub-tropical coastal waters (Worden et al., 2004; Kimmance et al., 2007; Fowler et al., 2020).

At the light levels used for the incubations, approximately 60% surface PAR, μ_max of both populations exceeded m, resulting in positive net growth rates (Figure 3). If it is assumed that the same coastal waters were sampled on each date, then the picoeukaryote population showed an increase in abundance between days ∼190 and 211 of 2019, equivalent to a net rate of growth of 0.05 d^–1 (Figure 3A). The similar pattern in Synechococcus abundance lagged the picoeukaryote population by approximately 5 d but increased at a faster net rate between days 193 and 211 of 0.15 d^–1 (Figure 3A). Even during these periods of maximal in situ net growth, the μ_net values derived from the incubations were considerably higher for both populations (Figure 3 and Table 2). The closer balance between gross growth and grazing mortality in the water column as a whole, may be a consequence of reduced overall light levels on gross growth rates caused by mixing of photoautotrophic cells through the water column, while grazing mortality rates may be more constant with depth.

The approach as presented here, was focused on specific groups of picophytoplankton whose abundance could be precisely quantified by flow cytometry. In this respect, the method may not be a suitable means of obtaining total phytoplankton grazing mortality in the way that the dilution approach has been applied. How chlorophyll concentrations changed in response to increasing bead addition was not measured but may provide estimates of grazing and growth rates when a single size spectrum dominates the phytoplankton. As illustrated previously and here, the dilution approach can also be used to generate phytoplankton group-specific information that can be directly compared to results from the saturation approach (Figures 3, 4; Landry et al., 1995; Worden and Binder, 2003).

The saturation approach relies on several assumptions that need to be considered and more thoroughly tested for applications in specific environments. The primary assumption is that the addition of surrogate prey does not influence μ_max of the natural prey. In preliminary laboratory experiments, 2 μm beads were added to six tubes containing only M. pusilla over a range of 0–6 × 10⁶ beads ml^–1, to test for an effect of the microspheres on growth. Under the illumination and mixing regime generally used for the culture experiments, growth rates of M. pusilla averaged 0.48 ± 0.03 d^–1 across the six tubes, with no apparent trend, indicating minimal influence of the surrogate prey on growth (not shown). Similar tests were not carried out for natural waters, where the separation of prey and predators is more difficult to accomplish.

There is the possibility that in nutrient limited waters, the reduced grazing caused by surrogate prey additions may decrease nutrient recycling rates and as a consequence, potentially reduce estimates of m and μ_max. A possible advantage of the saturation approach is that the biomass of the natural population is not altered across the incubation series as part of the initial experimental set-up. Population size does vary, nonetheless, between bottles in the series over the course of the incubations. In the natural water experiments, the maximum variations in population size between incubations exhibiting net growth (zero bead addition) and gross growth (μ_max) within an experiment, were 1.7-fold for picoeukaryotes on day 204 and 1.9-fold for Synechococcus on day 197 (Table 2). Thus, variation in nutrient availability between treatments of the saturation approach is likely to be lower than across the 10–5-fold difference in initial biomass that is required in dilution experiments.

A key aspect of the experimental design is judgment of the quantity of surrogate prey required to illicit a saturation response, while minimizing potential artifacts from too many additional particles. In the coastal waters of the Gulf of Maine, 2 μm bead additions of > 2 × 10⁶ ml^–1 saturated the grazers of picoeukaryote and Synechococcus populations, sufficiently to generate significant values of m and μ_max (Figure 2). It seems reasonable to expect the behavioral response of phagotrophic protists to be consistent in similar water types but this remains to be established. Preliminary studies in oceanic waters suggest grazers of Prochlorococcus and Synechococcus populations may tend toward saturation at bead additions of < 1 × 10⁶ ml^–1. In theory, coefficient a of Equation 2 varies in proportion to the rate of increase in μ_app relative to C, the surrogate prey concentration. It may be possible to use values of a, or an analogous coefficient, to judge how effectively different surrogate prey types or sizes lead to saturation of predators. As more experiments are carried out, it should be possible to standardize the surrogate prey type and levels of addition for different water types.

The asymptotic regression model interpretation of the response of μ_app to increasing surrogate prey, provided highly significant values of μ_max (p < 0.001) and μ_net (p < 0.01) and hence, estimates of m (Table 2 and Supplementary Table 2). However, multiple predator types may consume picophytoplankton populations, each having a different “functional response” to increased surrogate prey abundance and the interpretation of the observed experimental results may be considered over-simplistic. Future studies could usefully examine more complex model interpretations of the saturation approach that consider, for example, predator and prey motility and encounter rates, the process and duration of prey capture and selectivity by the predator (Flynn and Mitra, 2016).

The governing role that dissolved compounds play in microbial interactions is becoming ever more apparent (Strom et al., 2007; Stocker et al., 2008). Indeed, a contributing factor to inconsistent results generated by dilution approach is the alteration of seawater chemistry through the filtration step to create the diluent required in the experimental set-up. One specific example, is the release of polyunsaturated aldehydes to the dissolved phase and the resulting increased inhibition of phytoplankton apparent growth with increased levels of dilution (Stoecker et al., 2015). An advantage of the saturation approach is that modification of seawater chemistry may be kept to a minimum. Although the bovine serum albumin treatment of beads and subsequent resuspension in filtered seawater appeared to be sufficient for these experiments, more thorough cleaning protocols for artificial surrogate prey are likely to be required for low nutrient and oceanic waters. Nonetheless, it may be possible to investigate grazer-mediated production of dissolved organic compounds, inorganic nutrients and possibly, trace elements by using the approach to generate a gradient of grazing pressure.