What controls microzooplankton biomass and herbivory rate across marginal seas of China?

Microzooplankton are the primary herbivores and nutrient regenerators in the marine food web, but their importance is often underestimated, and the quantitative relationships between environmental factors and the biomass and herbivory rate of microzooplankton remain obscure. To fill this gap, we conducted 224 dilution experiments to measure microzooplankton biomass and herbivory rate across a vast area of the marginal seas of China. To gain the potential mechanisms controlling microzooplankton herbivory, we also use a model that combines the Metabolic Theory of Ecology and the functional responses of grazing to quantify the effects of temperature, phytoplankton biomass, and microzooplankton biomass on microzooplankton grazing rate. We estimate an activation energy of 0.51 eV of microzooplankton and found that the Holling III function best described the functional response of microzooplankton grazing with a maximal ingestion rate of 4.76 d−1 at 15°C and a half‐saturation constant of 0.27 μM N. We also find that microzooplankton biomass scales with phytoplankton biomass with an exponent of 0.77, consistent with the general 3/4 scaling law found in other ecosystems. This scaling relationship is accompanied by a shift from ciliates to heterotrophic dinoflagellates with increasing phytoplankton biomass. Our results provide empirical patterns that will be vital to parameterize and validate marine ecosystem models, particularly in China seas.

Despite their ecological importance, microzooplankton remain understudied. Routine oceanographic observations seldom monitor microzooplankton biomass or herbivory rate, although the dilution technique, an elegant method of measuring microzooplankton herbivory rate, has been developed for almost four decades (Landry and Hassett 1982). The number of observations of microzooplankton herbivory rate is around 1600 globally (Chen et al. 2012;Schmoker et al. 2013), far less than that of primary productivity (> 50,000; Buitenhuis et al. 2013). This makes validating and optimizing the grazing function of microzooplankton difficult in ocean ecosystem models. Indeed, microzooplankton have been rarely validated by observational data in mainstream ocean models (but see Buitenhuis et al. 2010). Promoted by such gap, Buitenhuis et al. (2010) suggested that "the most effective progress which can be made in defining the role of microzooplankton in global biogeochemical cycles is to make microzooplankton biomass a standard oceanographic measurement, in particular during dilution grazing experiments." We lack a thorough understanding of the relationship between microzooplankton biomass and environmental variables such as prey (phytoplankton) biomass and temperature. The slope of microzooplankton biomass against phytoplankton biomass reflects whether microzooplankton are under bottom-up or top-down control (Yuan and Pollard 2018). If microzooplankton are mostly controlled by bottom-up factors, microzooplankton biomass should increase in proportion with phytoplankton biomass. Conversely, if microzooplankton are mostly regulated by top-down effects by their predator (i.e., mesozooplankton), increases in phytoplankton biomass should lead to negligible increases of microzooplankton biomass. However, empirical evidence is mixed as to the biomass relationship between microzooplankton and phytoplankton. Hatton et al. (2015) have found that logarithmic predator biomass increases linearly with logarithmic prey biomass with the slope close to 3/4 in different ecosystems. Nevertheless, a universal scaling exponent has been challenged (Yuan and Pollard 2018). For instance, zooplankton biomass was almost constant with increasing phytoplankton biomass in eutrophic environments (Irigoien et al. 2004;Yuan and Pollard 2018). As such, it is still uncertain whether the relationship between microzooplankton biomass and phytoplankton biomass should be linear or curvilinear on a log-log scale.
There are also uncertainties as to whether temperature should affect microzooplankton biomass relative to phytoplankton biomass. Because the metabolic rate of zooplankton is believed to increase with temperature faster than that of phytoplankton, the biomass ratio of microzooplankton to phytoplankton should decrease with increasing temperature (Chen et al. 2012). However, this metabolic asymmetry has recently been challenged (Chen and Laws 2017). Thus, it remains unclear whether the microzooplankton-to-phytoplankton biomass ratio should vary with temperature or not.
In addition, we also lack clearly defined grazing functions of microzooplankton, especially for natural microzooplankton assemblages (Sandhu et al. 2019). Microzooplankton grazing rate is a complex function of the characteristics of predator and prey, such as prey quality and concentration, predator concentration and feeding strategies, and predator-prey size (Hansen et al. 1997). Among these factors, prey concentration is considered as major determinants of microzooplankton grazing rate in many biogeochemical/ecosystem models due to its relatively easy access (Buitenhuis et al. 2010;Laufkotter et al. 2016). The effect of prey concentration on microzooplankton specific grazing rate (i.e., per capita zooplankton per unit time) is described as the functional response. In general, the ingestion rate increases with increasing prey concentration linearly or nonlinearly to a plateau and the corresponding clearance rate decreases as prey concentration increases (Gentleman et al. 2003;Sandhu et al. 2019). Temperature is another factor that affects microzooplankton grazing activities and should be incorporated into the grazing functions (Hansen et al. 1997;Buitenhuis et al. 2010). Although the grazing functions of microzooplankton have been intensively studied in the laboratory (Hansen et al. 1997), few studies clearly quantified the effects of prey (phytoplankton) concentration and temperature simultaneously on natural microzooplankton assemblages (but see Buitenhuis et al. 2010), which should be worth more effort.
In this study, we fill these gaps by assembling a data set consisting of 224 microzooplankton biomass and herbivory rates that we measured through dilution experiments during nine cruises across China seas including the South China Sea, the East China Sea, and the Yellow Sea. These data expand previous studies of Chen et al. (2009Chen et al. ( , 2013, Su et al. (2014), , Zheng et al. (2015), and Liu et al. (2019) and cover the gradient from eutrophic coastal waters to oligotrophic oceanic waters. This rich data set provides us with the opportunity to address the following questions: (1) What are the most important factors affecting microzooplankton biomass? Specifically, what are the quantitative relationships between microzooplankton biomass and environmental variables such as phytoplankton biomass and temperature?
(2) What are the most important factors affecting microzooplankton community grazing rates? We will examine the empirical relationships between microzooplankton biomass and grazing rate with phytoplankton biomass and temperature. Then we will contrast these empirical relationships against theoretical predictions. We also investigate the role of mesozooplankton biomass and microzooplankton community via the fraction of heterotrophic dinoflagellates in the total community (Sherr and Sherr 2007) in affecting microzooplankton biomass and grazing rate. Our analysis will deepen the mechanistic understanding of what controls microzooplankton biomass and grazing rate in the ocean. The algorithms we have developed will also benefit parameterizing and validating ocean ecosystem models, particularly in China seas.

Methods and materials
Field data Field sampling and plankton sample analysis Microzooplankton samples and accompanying environmental data were collected during nine oceanographic cruises across the marginal seas of China including the South China Sea, the East China Sea, and the Yellow Sea (Supporting Information Table S1, Fig. 1). Samples were also collected during whole-year monthly observations at three coastal stations adjacent to the northern South China Sea (Supporting Information Table S1, Fig. 1). The stations covered a broad environment of China seas ranging from very eutrophic coastal waters to oceanic offshore waters in different seasons. The data were geographically classified into four subsets based on the study areas: the South China Sea, the East China Sea, the Yellow Sea, and Coastal stations. We grouped the data of monthly observations as one subset (Coastal stations) because the sampling stations were very close to the land and substantially affected by the river input ( Fig. 1).
We collected in situ data of temperature and chlorophyll a (Chl a) associated with microzooplankton samples. Water temperature and salinity were measured at each station by Seabird CTD probes during the cruises and a YSI multiprobe sensor (6600 or EXO2) during monthly observations. Surface seawater was collected using normal Niskin bottles attached to a CTD rosette system or an acid-washed plastic bucket. To measure Chl a concentration, aliquots of 100-500 mL seawater were filtered onto 25 mm GF/F glass fiber filters under low vacuum (0.2 μm polycarbonate membrane filters were used at stations in Hong Kong waters). All filters were stored at −80 C until further analysis. Following the protocol of Joint Global Ocean Flux Study, the filters were extracted in 90% acetone at −20 C or 4 C in the darkness and the Chl a concentrations were measured within 24 h using a Turner Designs fluorometer with a nonacidification module (Model No. Trilogy 040) (Welschmeyer 1994).
To collect the microzooplankton samples, seawater (100-500 mL) was gently siphoned from carboys into plastic amber bottles containing acidic Lugol's solution with a final concentration of 5%. The samples were stored in the dark at room temperature until further analysis. Aliquots of 10-100 mL of each sample were settled for 24 h using Utermöhl chambers (the samples with large volumes [500 mL] were processed by a two-stage settling process) and observed with an inverted microscope (Olympus IX51 or Leica Dmirb) at ×200 magnification. We primarily focused on phagotrophic protists larger than 20 μm, including ciliate and heterotrophic dinoflagellates, which mainly consume phytoplankton and are the major contributors to the grazing rate estimated by dilution technique (Calbet 2008). The general categories, including ciliates (aloricate ciliates and Tintinnids), dinoflagellates, and copepod nauplii (< 200 μm) were enumerated (Elangovan and Padmavati 2017). As the copepod nauplii cannot be spotted in most samples and their abundance cannot be accurately estimated through this method, we excluded this category when calculating the total abundance and biomass of microzooplankton. For each sample, at least 50 cells were counted.
We estimated microzooplankton biomass by summing up the biovolumes of all cells counted under the microscope. To estimate the cell volume, each cell was assigned to a geometrical shape, and their width and length were measured using the software SPOT (version 3.5) or Simple PCI6. The carbon content of each microzooplankton was estimated based on the corresponding empirical equations for acidic Lugol's fixed cells. For ciliates, a conversion factor of 0.19 pg C μm −3 was used on the aloricate ciliate (Putt and Stoecker 1989), and the equation of pg C cell −1 = 444.5 + 0.053 × biovolume was used on Tintinnids (Verity and Lagdon 1984). The biovolume of dinoflagellate was converted to cell carbon based on the equation of pg C cell −1 = 0.76 × biovolume 0.819 (Menden-Deuer and Lessard 2000). The average size of ciliates or dinoflagellates (pg C cell −1 ) was calculated by dividing their total biomass (pg C L −1 ) by the corresponding total abundance (cells L −1 ) in a sample.
We also assembled the data of mesozooplankton dry weight measured at the same time with our experiments at some stations and unpublished data). In brief, a total of 58 mesozooplankton samples were towed vertically from 200 m (or near the bottom for shallow water stations) to the surface using a plankton net (200 μm mesh size) attached by a digital flowmeter. The samples were filtered onto preweighted PC filters (20 μm, 47 mm) to measure mesozooplankton dry weight after oven-drying (24 h, 60 C). The carbon biomass of mesozooplankton was estimated by assuming carbon weight = 0.4 × dry weight (Steinberg et al. 2000).

Microzooplankton herbivory
Microzooplankton herbivory rates were estimated via the dilution technique, which is the standard method to estimate the microzooplankton herbivory rate on phytoplankton (Landry and Hassett 1982). Although the group of microzooplankton may include multicellular grazers such as copepod nauplii based on its technical definition (20-200 μm organisms), the grazing rate measured through this method was principally due to protistan grazers including ciliates and heterotrophic dinoflagellates Calbet 2008). Dilution experiments were conducted at every sampling station (Fig. 1). The measured volume of particle-free seawater, which was obtained by filtering the seawater through a 0.2 μm filter capsule (Pall Corporation) by gravity, was added into the 1.2-liter polycarbonate bottles. The bottles were then filled with the unfiltered seawater to the full capacity to get a mixture of certain percentages of unfiltered seawater with particle-free seawater. Five dilution treatments with unfiltered seawater percentages of 15%, 27%, 50%, 73%, and 100%, respectively, were established in all experiments during most cruises. In monthly observations in Hong Kong waters (2016/2017) and Xiamen Bay, two dilution treatments (15%/25% and 100%) were used, which has been proven as accurate as five dilution treatments (Chen 2015). Inorganic nutrients were added into all bottles to ensure constant phytoplankton growth. The final nutrient concentration in the bottles was 0.5 μmol L −1 NH 4 Cl, 0.03 μmol L −1 KH 2 PO 4 , 1 nmol L −1 FeCl 3 , and 0.1 nmol L −1 MnCl 2 for the oceanic stations, and 10 μmol L −1 NaNO 3 and 1 μmol L −1 KH 2 PO 4 for coastal stations. Two bottles filled with unfiltered seawater without nutrient addition were prepared for no-nutrient-addition controls. All bottles were tightly capped, incubated for 24 h in a deck incubator covered with neutral screens and cooled by running seawater during all cruises. All experimental carboys, filter capsules, bottles, and tubing were washed with 10% HCl and rinsed with distilled water and in situ seawater prior to each experiment. After incubation, Chl a samples were taken from each experimental bottle and analyzed using the method mentioned above.
The microzooplankton grazing rate was estimated following Landry et al. (2008). Assuming exponential growth for phytoplankton in each bottle, the net growth rate of phytoplankton (g, d −1 ) was calculated as g = (1/t)ln(P t /DP 0 ), where P 0 and P t were the phytoplankton biomass represented by Chl a concentration before and after incubation. D is the dilution factor (the percentages of unfiltered seawater) of each bottle. The mortality grazing rate due to microzooplankton (m, d −1 ) was estimated as the slope of the linear regression of the net growth rate against the dilution factor. At coastal stations, saturated grazing was observed occasionally. When it occurred, the grazing rate was estimated as the slope of the regression curve within the nonsaturation range (Gallegos 1989). The grazing rate was removed when the slope of the regression was positive.

Effects of prey concentration on microzooplankton biomass: Predator-prey power law
The relationship between predator and prey biomass has been described by a power-law relationship and further transformed into the linear regression model after logarithmic transformation (Hatton et al. 2015): where c is the coefficient and γ is the dimensionless scaling exponent; B z (μg C L −1 ) is the carbon biomass of total microzooplankton; P (μg C L −1 ) is the phytoplankton carbon biomass converted by multiplying the Chl a concentration by the model-derived carbon-to-chlorophyll ratio (C : Chl, gC gChl −1 ). The C : Chl ratio was estimated using a Boosted Regression Trees model (Elith and Leathwick 2017) which was constructed based on a published data set consisted of C : Chl ratio, light, temperature, and dissolved inorganic nitrogen (Chen and Liu 2010). In this model, the C : Chl ratio is a function of temperature, nutrient, and light (Supporting Information). The exponent γ < 1 means the trophic pyramid is more bottom-heavy at higher biomass and indicates bottom-up control of prey on the predator (Hatton et al. 2015). Ordinary least squares (OLS) was used to fit the log-transformed data following Hatton et al. (2015), which is commonly used in fitting bivariate power laws (Del Giorgio and Gasol 1995). We also tried to fit a second-order term in the linear regression to test whether curvilinear can fit the data better than linear: where both γ 1 and γ 2 are regression coefficients.

Effects of temperature on microzooplankton biomass
We explored the theoretical prediction of how temperature affects the ratio of microzooplankton to phytoplankton biomass. We assumed that at steady state, microzooplankton production (MP) is a constant fraction of primary production (Pμ) . This fraction is the gross growth efficiency (GGE).
where G z and μ are the growth rates (d −1 ) of microzooplankton and phytoplankton, respectively. The growth rate of both microzooplankton and phytoplankton are affected by temperature, which can be described by Boltzmann-Arrhenius equation in the Metabolic Theory of Ecology (Brown et al. 2004): in which G z0 and μ 0 are normalization constants for the growth rates of microzooplankton and phytoplankton, respectively. E p and E zg are the activation energies (eV) of the growth rates of phytoplankton and microzooplankton, respectively, which describes the temperature sensitivity of rates, k is the Boltzmann's constant (8.62 × 10 −5 eV K −1 ), T is temperature (K), and T 0 is the reference temperature (288 K). As the activation energy of microzooplankton growth rate (E zg ) is around twice of that of phytoplankton (E p ) (Chen et al. 2012), we should have: which means that the biomass ratio of microzooplankton to phytoplankton should decrease with increasing temperature if we assume that the gross growth efficiency does not vary with temperature (Chen et al. 2012).
Modeling microzooplankton grazing rate based on metabolic theory of ecology and functional response of grazing The microzooplankton grazing rate (m, d −1 ) estimated by the dilution approach is the microzooplankton community clearance rate, which is the sum of clearance rate of each grazer in a unit volume of water parcel: where F i (L grazer −1 d −1 ) is the clearance rate of the i th grazer. N is the abundance of grazers in the community. The clearance rate can be calculated by dividing the ingestion rate I i (μg C grazer −1 d −1 ) by the mean food concentration P i (μg C L −1 ) based on their definitions (Båmstedt et al. 2000).
Although the grazing activities of each grazer could be affected by many factors such as the grazing pressure by their consumers and predator-prey size ratio, we consider an ideal situation where the ingestion rate of each grazer is only dependent on prey concentration and can be described by the maximum ingestion rate and functional response as in many biogeochemical models (Hansen et al. 1997;Laufkotter et al. 2016): where I m,i (μg C grazer −1 d −1 ) is the maximum ingestion rate of the i th grazer, and f(P i ) is the functional response of each grazer describing the relationship of ingestion rate and food concentration. According to the metabolic theory of ecology, the maximum ingestion rate of each grazer can be calculated by a function of body size and temperature, thus Eq. 7 becomes: where I 0 (d −1 ) is a normalization constant; e Ez 1 kT 0 − 1 kT describes the effect of temperature on the maximum ingestion rates in which E z is the activation energy of microzooplankton grazing rate and other symbols are the same with Eq. 4; M i (μg C grazer −1 ) is the body mass of the i th grazer and α is the allometric exponent for the body mass. When the M i is in the unit of carbon content, α could approximately equal to 1. Then the sum of the α-scaled body mass ( approximatively equal to the total biomass of the grazers (B z , unit: μg C L −1 ) (Hansen et al. 1997;Chen et al. 2012). To explore the effect of prey concentration on the community grazing rate of microzooplankton, we assumed a system-level functional response of ingestion rate vs. the prey concentra- , which is the aggregate behavior of different grazers and preys under different environmental conditions ranging from eutrophic coastal water to the oligotrophic offshore waters. Then the Eq. 8 becomes: Thus the biomass-specific grazing rate (m/B z , unit: d −1 [μg C] −1 L) should be determined by a function of temperature and prey concentration: By taking logarithms on both sides, we obtain: There are several mathematical formulations for the functional response of grazers f(P). In the present study, we examined the three most common nonlinear formulations (Gentleman et al. 2003;Sandhu et al. 2019) to find the most suitable one for our data set: Holling II : Holling III : in which K p1 , K p2 , and K p3 are constants (unit: μg C L −1 ). K p2 and K p3 are half-saturation constants. P is the food concentration. As phytoplankton is the major food source of microzooplankton , and the grazing rate estimated through dilution approach is the phytoplankton mortality grazed by microzooplankton, the food concentration here is the phytoplankton concentration in carbon unit (μg C L −1 ). The fitting of nonlinear formulation combining Eqs. 11 and 12, 13, or 14 to data was implemented using the R function "nls" and package "nlstools" in the software R 3.4.3 (Baty et al. 2015). We also used generalized additive models (GAMs) to describe the variation in ln(m/B z ) with the log-transformed phytoplankton carbon biomass as a smoother term to demonstrate the patterns without parametric constraints (Wood 2006). The GAMs analysis was implemented using the function "gam" in the R package "mgcv" (Wood 2006).

Environmental controls on microzooplankton biomass
We find that microzooplankton biomass is significantly correlated with phytoplankton biomass (Spearman r = 0.61, p < 0.001; Fig. 2; detailed results for each experiment are summarized in Supporting Information Table S3). The relationship between microzooplankton and phytoplankton biomass can be described by a linear regression model on a log-log scale (Eq. 1; Fig. 3a, adjusted R 2 = 0.36, n = 223, p < 0.001). The scaling exponent γ is 0.77 AE 0.07, similar to the scaling exponent (near 3/4) found in the universal relationship between predator and prey biomass (Hatton et al. 2015). A second-order term is insignificant in the regression (Eq. 2, ANOVA, p > 0.05). The biomass of major groups of microzooplankton also increases with increasing phytoplankton biomass (Fig. 3b). The scaling exponents for the relationship between ciliates and heterotrophic dinoflagellates vs. phytoplankton biomass are 0.73 AE 0.07 and 0.83 AE 0.11, respectively.
No correlation between temperature and microzooplankton biomass is observed (Spearman r = −0.06, n = 223, p > 0.05; Fig. 2). The biomass ratio of microzooplankton to phytoplankton does not correlate with temperature either (Supporting Information Fig. S1).
In summary, we find that the relationship between microzooplankton biomass and phytoplankton biomass is well consistent with the universal 3/4 scaling law (Hatton et al. 2015) and does not show any curvilinearity as shown by Irigoien et al. (2004) and Yuan and Pollard (2018). Inconsistent with the theoretical prediction, we cannot find any temperature effect on microzooplankton biomass.

Environmental controls on microzooplankton grazing rate
As predicted by Eq. 11, temperature, the total biomass of microzooplankton and phytoplankton all play important roles in affecting the microzooplankton community grazing rate derived from the dilution experiments. The community microzooplankton grazing rate strongly increases with temperature (Spearman r = 0.31, p < 0.001), phytoplankton biomass (Spearman r = 0.29, p < 0.001), and microzooplankton biomass (Spearman r = 0.35, p < 0.001), but not correlated with microzooplankton community composition which is roughly represented by the fraction of heterotrophic dinoflagellates (HDF.ratio, Spearman r = 0.07, p > 0.05; Fig. 2). When phytoplankton biomass is low (0~about 25 μg C L −1 ), the biomassspecific grazing rate (m/B z ) increases with phytoplankton biomass (Fig. 4a, Supporting Information Fig. S2). Above 25 μg C L −1 , m/B z decreases with increasing phytoplankton biomass. m/B z also increases with temperature (Fig. 4b,  Supporting Information Fig. S2).
To investigate which form of functional response performs the best, we further fit the experimental data to Eq. 11 with three types of functional responses (Eqs. 12-14; Table 1). All parameter estimates are significant (the corresponding p values are shown in the table). The model with Holling III function explains 27% variation of biomass-specific grazing rate measurements, which is better than those with Ivlev and Holling II functions (differences of AIC are about 4; Table 1). The nonparametric GAMs fit with the log-transformed phytoplankton biomass as a smoother term also shows similar patterns of ln(m/B z ) variation (Deviance explained = 30%, n = 208; Supporting Information Fig. S2). Thus, the model rate), ln mesozooplankton biomass (Mesozoo) (n = 58), and microzooplankton community composition represented by fraction of heterotrophic dinoflagellates (HDF.ratio). The lower-off diagonal shows the scatter plots between every two factors with correlation ellipses and Lowess regression lines (red lines), the red dots are the points of the mean of x and y; the diagonal shows the histograms of every factor; the upper-off diagonal reports the spearman correlations with significant levels (*p < 0.5; **p < 0.1; ***p < 0.001). Noting that the HDF.ratio is arcsine square-root transformed ratio (arcsin ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi HDF−ratio p ).
with Holling III function was used to predicted m/B z and community grazing rate (

Changes of microzooplankton community composition
We examined whether the biomass ratio of heterotrophic dinoflagellates to total microzooplankton changed with trophic status (Sherr and Sherr 2007). We found that although both the biomass of ciliates and heterotrophic dinoflagellates increased with increasing phytoplankton biomass (Fig. 3b), (c) Relationship between the biomass of microzooplankton/phytoplankton/total prey (microzooplankton + phytoplankton) and mesozooplankton. All solid lines are fitted by the OLS linear regression. γ is the slope of the regression model (the scaling exponent in Eq. 1), n is the number of data used in the corresponding regression models. there was no correlation between the fraction of heterotrophic dinoflagellates and phytoplankton biomass (Fig. 2). The fraction of heterotrophic dinoflagellates is significantly higher at the Coastal stations than in other regions because the heterotrophic dinoflagellates are very abundant in the eutrophic waters (Fig. 6a). While in the mesotrophic and oligotrophic waters, the fraction of heterotrophic dinoflagellates is relatively low and shows no correlation with Chl a concentration (Spearman r = −0.09, n = 123, p > 0.05; Fig. 6a). In these regions, ciliates are the major consumers of phytoplankton. The average size of ciliates is positively correlated with Chl a concentration in the regions of the South China Sea, the East China Sea, and the Yellow Sea (Spearman r = 0.3, n = 151, p < 0.001; Fig. 6b).

Discussion
One major mission in ecological study is to seek for general laws underlying patterns across scales and ecosystems (Lawton 1999). It is particularly striking that some general patterns and regularities can be found in the ecosystems despite significant environmental and measurement noises. In the course of pursuing the quantitative insight into what controls microzooplankton biomass and herbivory, we have found supporting evidence for three remarkable general patterns across ecosystems which could benefit a better understanding of the role of microzooplankton in the marine food web and for parameterizing and validating biogeochemical models.
We have found that the biomass relationship between microzooplankton and phytoplankton follows a simple power-law function with a scaling exponent of 0.77 AE 0.07 close to 3/4 (Fig. 3), which is highly consistent with the universal scaling law of predator and prey biomass in diverse terrestrial and aquatic systems (Hatton et al. 2015). Similar exponent is also obtained for the relationship between mesozooplankton biomass and their total prey biomass (0.67 AE 0.12, Fig. 3c), serving as evidence for the predator-prey Table 1. Estimated parameters of models incorporating MTE with three functional responses of grazing: Ivlev, Holling II, and Holling III. I 0 (d −1 ): The maximum ingestion rate of microzooplankton at reference temperature (15 C); E z (eV): Activation energy of biomassspecific microzooplankton grazing rate; K p (μg C L −1 ): Half-saturation constant of microzooplankton grazing; K p (μM N): Half-saturation constant of microzooplankton grazing in the unit of μM N converted by a Redfield C : N ratio of 6.625; AIC: Akaike information criterion of each model. The values in brackets are the 95% CI. power law again. It has been found surprisingly that from small zooplankton to large mammals, as long as they are predators, their biomass exhibits a robust near 3/4 scaling with their total prey biomass. Furthermore, the universal scaling exponents of these large-scale patterns are very similar to the body mass allometries supported by theoretical explanations based on general features of metabolic networks (Brown et al. 2004). However, the potential mechanisms linking largescale ecosystem patterns with finer-grain processes remain elusive in ecology (Hatton et al. 2015). Although our study cannot provide any direct clues for the underlying mechanisms, the striking pattern observed in our data set serves as strong supports for the predator-prey 3/4 power law.

Model type
In the increasing pattern between logarithmic microzooplankton and phytoplankton biomass, no saturation at higher biomass was observed, which indicates that the microzooplankton biomass is mainly determined by the bottom-up control even in the eutrophic coastal waters in China seas. The negligible microzooplankton biomass increase at higher phytoplankton biomass, which was found in the previous studies (Irigoien et al. 2004;Yuan and Pollard 2018), should only occur in the extremely eutrophic environments particularly in the lakes (Yuan and Pollard 2018) or occur during phytoplankton blooms when the enlarged size and worsened quality of phytoplankton help them escape from the grazing of microzooplankton (Irigoien et al. 2005). In general, the microzooplankton biomass is mostly determined by phytoplankton biomass in the ocean, although it is also affected by the top-down control from mesozooplankton.
Mesozooplankton not only feeds on microzooplankton but also competes with them for the phytoplankton food. In the current study, we found that the mesozooplankton biomass was significantly correlated with microzooplankton biomass and grazing rate (Fig. 2), suggesting their influence on microzooplankton. Nevertheless, we also found the mesozooplankton biomass increased with phytoplankton biomass with a larger scaling exponent (Fig. 3c), which indicates that mesozooplankton may feed on more phytoplankton than microzooplankton. Under such circumstances, microzooplankton will suffer less grazing pressure from mesozooplankton but compete with them. Thus, the bottom-up control by phytoplankton should be prevalent than top-down controls. In fact, it is difficult to quantify the top-down control on microzooplankton because it is highly dynamic and depends on the biological conditions in the study area including the community compositions of mesozooplankton and the characteristics of other preys such as phytoplankton (Broglio et al. 2004). For instance, more grazing pressure on microzooplankton is found in the oligotrophic ocean, because the small-sized phytoplankton is dominated and rarely consumed by copepods (Calbet 2008). While the mesozooplankton is also found to shift their predation from phytoplankton to microzooplankton over the course of blooms as the quality of phytoplankton food decreases, which leads to a strong top-down control on microzooplankton (Löder et al. 2011). As such, more investigations are highly necessary to well define the top-down control by mesozooplankton (Calbet 2008). For the effect of temperature on microzooplankton biomass, we cannot find the biomass ratio of microzooplankton to phytoplankton decrease with increasing temperature as in previous studies (Chen et al. 2012). Conversely, temperature shows little effect on the microzooplankton biomass (Fig. 2,  Supporting Information Fig. S1). According to Eqs. 3, 5, the expected effect of temperature on microzooplankton biomass is based on one assumption that the gross growth efficiency of microzooplankton is temperature insensitive. However, the gross growth efficiency of some protists can increase with increasing temperature (Mayes et al. 1997), which would counteract the downward trend between biomass ratio and temperature. Nevertheless, how temperature affects the gross growth efficiency of protist is still controversial as decreased and constant gross growth efficiency with increasing temperature had also been reported for protists (Rose et al. 2009). Thus, studying the temperature sensitivity of gross growth efficiency of microzooplankton could be a prerequisite for understanding how temperature affects the microzooplankton biomass.
We have also found that microzooplankton herbivory is highly regulated by temperature and prey concentration (Figs. 4,5), which can be well quantified following the remarkable regularities in ecology, that is, the metabolic theory of ecology and the functional response of grazing. The metabolic theory of ecology reveals a regularity in scaling the metabolic rates of organisms with their body size and temperature. Since all organisms share the same fundamental metabolic reactions for energy transformation and biosynthesis, this regularity has been scaled up from the individual to population, community, and the ecosystem properties such as nutrient flux (Enquist et al. 2015). Despite disputes on the adequacy of mechanistic basis (O'Connor et al. 2007), the metabolic theory of ecology provides a quantitative framework allowing us to study how temperature affects the microzooplankton grazing rate. The activation energy of microzooplankton biomassspecific grazing rate estimated by the model with Holling III function is 0.51 eV (Table 1, 95% CI = 0.23-0.79 eV), which is not significantly different from the predicted value by the metabolic theory of ecology (0.65 eV). Our estimate was consistent with the results of previous regional Liu et al. 2019) and global studies (Chen et al. 2012), explicitly elucidating the general pattern that microzooplankton grazing rate increase with increasing temperature with average activation energy predicted by the metabolic theory of ecology (Figs. 4,5).
Functional response of grazing is another recognized regularity in ecology, albeit several different mathematical formulations had been proposed (Gentleman et al. 2003). In our study, although the natural microzooplankton communities in various locations are highly diverse and involved in complicated trophic interactions, the bulk biomass-specific microzooplankton grazing rate varies with phytoplankton concentration following the generally observed functional response (Gentleman et al. 2003;Sandhu et al. 2019). The biomass-specific grazing rate increases with increasing phytoplankton concentration (0~about 25 μg C L −1 ) and started to decrease above 25 μg C L −1 (Fig. 4, Supporting Information  Fig. S2). This pattern can be well depicted by the Holling III function (Figs. 4,5). Compared with the Holling II response, the Holling III response has been less used in the models in previous studies, which would overlook the response of grazing at very low food concentrations (Sarnelle and Wilson 2008). Our results suggested that Holling III response is more accurate to depict the functional response of microzooplankton grazing in nature as the prey concentration varies in a broad range among different ecosystems.
The accurate description of grazing functional response is crucial for modeling the dynamics of planktonic ecosystems. As conceptual modeling approaches usually put microzooplankton and their prey in "boxes" in the model structure, an overall functional response describing how the bulk microzooplankton grazing rate varies with prey concentration is more applicable for food web models, which entails estimation using the natural microzooplankton assemblages (Sandhu et al. 2019). Due to the scarcity of field rate measurements, the parameters used in many biogeochemical models were usually from the laboratory experiments on single species grazing. Recently, some studies have parameterized the functional responses using field data sets of dilution experiments (Buitenhuis et al. 2010). By contrast, our study not only used the natural microzooplankton communities but also incorporated the effect of temperature into the grazing functions since temperature also exerts a profound influence on microzooplankton grazing rate (Hansen et al. 1997;Chen et al. 2012). The parameters in our models, therefore, could be closer to natural situations. The maximum ingestion rate (I 0 ) in our models are all within the range of maximum ingestion rate of individual ciliates species (2.4-11.5 d −1 at 20 C), but a bit higher than that of individual heterotrophic dinoflagellates (0.26-4.08 d −1 at 20 C; Hansen et al. 1997). The high ingestion rate could be due to the fact that ciliates are dominant in the microzooplankton community in the study area as we found the fraction of heterotrophic dinoflagellates in all four regions is relatively low (Fig. 6a). To compare with the grazing rate derived from dilution approach, we converted the maximum ingestion rate to the clearance rate of microzooplankton community based on their definitions (Båmstedt et al. 2000) using the mean carbon biomass of phytoplankton and microzooplankton (121.75 μg C L −1 and 13.47 μg C L −1 , respectively). For such community, the grazing rate is 0.53 d −1 at 15 C (Holling III) which is close to the average grazing rate estimated from dilution experiments in the coastal and estuarine systems . The half-saturation constants (K p ) were converted to μM N equivalents assuming a Redfield ratio of carbon to nitrogen (C : N) being 6.625 (Table 1). The estimates are all within the ranges of K p for individual species obtained from laboratory experiments (ciliates: 0.06-2.46 μM N; heterotrophic dinoflagellates: K p2 = 0.2-11 μM N at 20 C; Hansen et al. 1997). Whereas our results suggested a lower K p (0.27 μM N) for Holling III function than that used in some ecosystem models, which were usually 0.5 μM N (Chen and Smith 2018).
Built on the metabolic theory of ecology and the functional response of grazing, our model illustrates the effect of temperature and phytoplankton biomass on the biomass-specific grazing rate quantitatively, although the variation of biomassspecific grazing rate the model can explain is relatively low ( Table 1). The explanatory power of the GAMs without parametric constraints is also limited (30%; Supporting Information Fig. S2). The explanatory power of the model could increase if more relevant factors such as the prey size and the potential effects of mesozooplankton were incorporated. However, these factors are hard to be quantified in the field experiments and usually tangling together in nature. In addition, potential sources of errors would undermine the accuracy of parameter estimates and model performance. One potential problem lies in the estimates of grazer biomass. The small metazoan such as copepod nauplii which are potentially important contributors to microzooplankton biomass in the coastal waters were not included in the grazer biomass because the sampling method hinders an exact estimate for their abundance. Similarly, due to the methodological issue, we included some mixotrophic flagellates and dinoflagellates in the total grazer biomass as they cannot be unequivocally distinguished from obligate heterotrophs based on the analysis of Lugol's-preserved samples. The mixotrophs are vague components of microzooplankton of which the contributions to the community grazing rate are difficult to be gauged. In addition, the nanoflagellates less than 20 μm were also excluded from the grazer biomass in our study. Although this grazer group does not belong to microzooplankton according to the technique definitions, they could be the major consumers of picophytoplankton and significantly contribute to the grazing mortality assessed by dilution experiments in the offshore waters. However, the percentage they account for in the total microzooplankton biomass should be small due to their small size. The other potential issue existed in the estimates of phytoplankton biomass. Due to the lack of the data of phytoplankton biomass, we used model-derived C : Chl ratio which is a function of temperature, nutrient, and light (Supporting Information). Although this approach will introduce errors into the model, it can reflect the variation of the phytoplankton biomass better than a constant C : Chl ratio. Furthermore, a common issue in estimating microzooplankton community grazing rate is that not all preys are equally accessible to all microzooplankton involved in the estimates. The model, therefore, suffers from the above methodological limitations and could be improved with more accurate estimates of grazer and prey biomass from more detailed information related to plankton community compositions.
Despite the highly diverse environmental conditions and potential measurement bias, our data set reveals certain patterns which are consistent with the regularities in ecology including the predator-prey power law, the metabolic theory of ecology, and the functional response of grazing. These patterns and regularities provide useful avenues for studying and predicting the dynamics of marine microbial food web and how it will respond to climate changes. They also pose significant challenges for modeling the functioning of ecosystems. Ecosystem model outputs should be consistent with the general patterns observed from empirical studies. For instance, the global biogeochemical model outputs should be able to reproduce the power-law relationships between predator and prey biomass. In addition, the regularities in ecology offer frameworks for quantifying the important compartments and processes in ecosystems. Based on the metabolic theory of ecology, the functional response of grazers, and the predatorprey power laws, we studied the primary factors determining the microzooplankton grazing rate and biomass in China seas and constructed models to well quantify them (Figs. 4,5). Thus, our models supported by the regularities in ecology can be used for validating regional biogeochemical models.
In addition to the factors examined under the framework of certain regularities in ecology, attempts had also been made to explore the effects of microzooplankton community composition on their biomass and herbivory. The microzooplankton community composition was found to be uncorrelated with microzooplankton grazing rate (Fig. 2). It is probably because that the fraction of heterotrophic dinoflagellates is a rough index of microzooplankton community composition without detailed information at the species level. Hence, we did not incorporate it in the above model. The fraction of heterotrophic dinoflagellates increased with increasing Chl a concentration, which is consistent with previous studies (Sherr and Sherr 2007) and implies a shift in microzooplankton community composition from ciliates to heterotrophic dinoflagellates across the gradients of trophic status. In the mesotrophic and oligotrophic regions (the South China Sea, the East China Sea, and the Yellow Sea), ciliates were the dominant groups of the microzooplankton community (Fig. 6). The community structure of ciliates in these regions also responds to the change of trophic status as their average size increased with Chl a concentration (Fig. 6b). There could be a shift in community structure of ciliates from small aloricate ciliates to large loricate ciliates (i.e., Tintinnids) with increasing Chl a concentration.

Conclusion
Our results revealed some general patterns that are consistent with the metabolic theory of ecology, the functional response of grazing, and the predator-prey power law, despite the geographical distance and highly dynamic environments, as well as methodological limitations. These general patterns provide us with common conceptual frameworks to understand the interactions between phytoplankton and their grazers and predict how the plankton ecosystem will respond to climate changes. Through this study, we gained comprehensive and quantitative insights into what factors control microzooplankton biomass and herbivory across China seas. The model derived from field measurements can be used to predict the microzooplankton biomass and grazing rates in the ocean, and to parameterize and validate ocean ecosystem models, particularly in China seas.