Effect of phytoplankton size diversity on primary productivity in the North Pacific: trait distributions under environmental variability

Abstract While most biodiversity and ecosystem functioning (BEF) studies have found positive effects of species richness on productivity, it remain unclear whether similar patterns hold for marine phytoplankton with high local richness. We use the continuous trait‐based modelling approach, which assumes infinite richness and represents diversity in terms of the variance of the size distribution, to investigate the effects of phytoplankton size diversity on productivity in a three‐dimensional ocean circulation model driven by realistic physics forcing. We find a slightly negative effect of size diversity on primary production, which we attribute to several factors including functional trait‐environment interactions, flexible stoichiometry and the saturation of productivity at low diversity levels. The benefits of trait optimisation, whereby narrow size distributions enhance productivity under relatively stable conditions, tend to dominate over those of adaptive capacity, whereby greater diversity enhances the ability of the community to respond to environmental variability.

The equation for phytoplankton biomass is: where µ com is the phytoplankton specific growth rate (d -1 ) of the whole community (i.e. integrated over the whole size spectra). The operators Adv() and Diff() describe three dimensional water advection and diffusion, respectively. The Arrhenius function , = Holling type III functional response of zooplankton grazing. I max is the maximal per capita ingestion rate (d -1 ), and K P the half-saturation constant of grazing.
Following Merico et al. (2014) and Smith et al. (2016), the equations for µ com , l, and v are: the second derivative of zooplankton clearance rate (d -1 ) against phytoplankton size and can be derived from the log-normal prey size distribution underlying our moment approximation and a log-normal size size grazing kernel consistent with empirical data and mechanistic considerations (Wirtz 2014). This grazing kernel F[ QFQ \]^J introduces the selectivity s as the inverse feeding breadth. The KTW mechanism implies that the optimal prey size (l opt ) of zooplankton grazing always coincides with the mean (log) size of the phytoplankton size distribution.
The second order derivative of the specific rate then fulfils the condition (Smith et al. 2016): in which α g > 0 means zooplankton is engaged in selective grazing and has a feeding preference on abundant prey. 6 The total ingested food by zooplankton is divided into net growth, excretion into the inorganic nitrogen pool, and defecation of unassimilated food into the detritus pool (Buitenhuis et al. 2010).
Zooplankton mortality is set to be proportional to the squares of its biomass and is also converted into detritus pool. As such, the dynamics of zooplankton follow: where NGE is the net growth efficiency of zooplankton. m z is the mesozooplankton mortality coefficient Detritus is converted to DIN at a rate (R dn , d -1 ) that has the same temperature sensitivity with zooplankton grazing. Detritus is also assumed to have a sinking rate (W d , d -1 ) that linearly relates to surface temperature (SST) (Gregg 2008).
where unass represents the fraction of unassimilated food by zooplankton. The sources and sinks of fer largely follow DIN with an additional source (atmospheric deposition; Fe depo ) and sink (scavenging; fer scav ) (Nickelsen et al. 2015): To translate between nitrogen and iron in phytoplankton and zooplankton, a constant fer:N ratio (R fer_N ) of 0.0265 is assumed. The data of monthly atmospheric deposition of total soluble iron are extracted from the Scenario III in Luo et al. (2008). Iron scavenging rate (fer scav ) is composed of both linear scavenging rate (k scm , d -1 ) and particle absorption rate (k sc , d -1 (µM N) -1 ) (Nikelsen et al. 2015): in which Fe prime is the concentration (nmol L -1 ) of free iron: where k eq (unit: (mol lig L -1 ) -1 ) is the equilibrium constant between free iron and ligands and is assumed to depend only on temperature: Note that T is in absolute temperature (K). l fe is the total iron ligand concentration that is assumed constant (0.6 nmol lig L -1 ).
The equation for D Fe is: Phytoplankton growth rate (µ) depends on temperature (T, K), light (I, W m -2 ), DIN and fer: The trait parameters . , K N , K fer, and α c are all dependent on cell size l: Phytoplankton chlorophyll-to-carbon (θ, g Chl (mol C) -1 ) and nitrogen-to-carbon (Q N , mol N (mol C) -1 ) ratios are calculated from ambient light and nutrient levels: where θ min and θ max are minimal and maximal Chl:C ratios, respectively. Q min and Q max are minimal and maximal N:C ratios, respectively. Total Chl a concentrations (Chl, µg L -1 ) and net primary production (NPP, µgC L -1 d -1 ) integrated over the whole size range can be calculated as: Light levels (I z ) at depth z were calculated based on PAR 0 and Chl a concentrations following the Beer-Lambert law: in which K w and K chl are the attenuation coefficients for seawater and Chl a, respectively. To realistically estimate the average light field that a phytoplankton cell should experience in the surface mixed layer (Franks 2015), the ambient light level for phytoplankton is calculated as the average light throughout this layer. All the parameter descriptions and values are shown in Table S1.
The initial conditions (January 1 st ) of inorganic nitrogen were interpolated from the WOA2013.
The initial conditions of phytoplankton biomass were estimated from the SeaWIFS chlorophyll monthly climatology with the chlorophyll-to-carbon and nitrogen-to-carbon ratios calculated from the ambient nitrate, light, and temperature following the model developed by Pahlow et al. (2013). The vertical profiles of chlorophyll were calculated based on surface values following Morel and Berthon (1989). The initial conditions of zooplankton and detritus were assumed 0.78 and 0.2 times the phytoplankton biomass, respectively, following previous model results (Ward et al. 2012). Initial and v were set globally set to 1 log µm 3 (approximating an equivalent spherical diameter of 1.7 µm) and 1 (log µm 3 ) 2 , respectively. The initial Fer conditions were extracted from the output of the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) ocean biogeochemical model that has been validated against iron fertilization experiments (Aumont & Bopp 2006). DETFe was initialized as 10 -3 (mol Fe (mol N) -1 ) of phytoplankton biomass.
We admit several limitations of our model. For example, N 2 fixation is not taken into account. We already know that in terrestrial BEF experiments, the complementarity between N 2 fixing legumes and other plants is an important factor leading to the positive BEF relationship. While N 2 fixation is certainly important in some parts of the North Pacific, the active process of denitrification due to the low oxygen environment essentially maintains the subtropical North Pacific limited by nitrogen. This suggests that we need to incorporate denitrification as well as oxygen cycle into the model if we plan to include N 2 fixation, which will add more complexity and uncertainty. Another limitation is that although we included two different nutrients, dissolved nitrogen and iron, phytoplankton species with different sizes do not differ in their uptake ratios of N: Fe according to the common parameterizations of size dependency of nutrient uptake on nitrogen and iron (Ward et al. 2012). Thus, we are essentially simulating a phytoplankton community competing for a single set of nutrients, which however does not imply that different size classes cannot coexist.

S.1.3 Comparisons between model outputs and observational data
We compared model outputs of dissolved inorganic nitrogen (DIN, µM), Chl a concentrations (µg L -1 ), net primary production (NPP, µg C L -1 d -1 ) and fractions of picophytoplankton (<2 µm) from the simulation using the largest "trait diffusion" (TD) coefficient with observational data (Fig. S2). The outputs from the simulation using the largest "kill-the-winner" (KTW) coefficient are largely similar with those using the largest TD coefficient. Reducing the coefficients of TD or KTW does not have a large impact on bulk properties such as nutrient, Chl a, or NPP, but has some effects on size structure (see The subarctic North Pacific is characterized by low temperature, relatively low light, high nitrate, large phytoplankton size and high biomass. Although the atmospheric dust deposition is not low in the subarctic Pacific, abundance in nitrate increases iron limitation of phytoplankton growth. But in general, light makes the main limitation factor in the subarctic and Chl:C and N:C ratios are generally high.
The subtropical Pacific the most oligotrophic area in the whole North Pacific, characterized by low nitrogen concentration due to strong water column stratification. The phytoplankton community is dominated by picophytoplankton. Light levels usually allow for sufficient productivity, also due to a shallow mixed layer. Due to the dual effects of strong light and low nitrogen, Chl:C and N:C ratios are low.
The eastern equatorial Pacific is characterized by strong upwelling, which leads to deep mixing and transport of nitrate into surface waters. This also leads to the paradoxical pattern that although surface light is strong in the equatorial regions, phytoplankton can be light-limited due to deep mixing. Also due to the scarce dust deposition, iron instead of nitrogen generally limits phytoplankton growth and particularly phytoplankton size. Therefore, the eastern equatorial Pacific is similar to the subarctic in many aspects except for temperature. Consequently, the highest primary production is found in this area, consistent with the output of the cbpm algorithm.

S.1.5 Temporal changes of phytoplankton size and productivity in two selected hotspots
To more thoroughly understand the effects of size diversity on productivity, we select two hotspots where the effects of size diversity on productivity were the most positive (hotspot A) and negative (hotspot B), respectively (Fig. S5). The hotspot A was located within the North Pacific subtropical gyre, while the hotspot B was in the eastern equatorial Pacific.
At hotspot A, the QVQ (i.e. the growth rate of the dominant species) in the high diversity treatment was usually higher than those in the low diversity treatment, suggesting that the size of the dominant species was closer to the optimal size in the highest than in the lowest diversity treatment (Fig.   S6). This behavior is inherent to our model equations and is consistent with the insurance hypothesis that the high diversity confers the phytoplankton community to more rapidly adapt to the changing environment by maintaining more species (Yachi & Loreau 1999;Norberg et al. 2001;Merico et al. 2014;Smith et al. 2016;Vallina et al. 2017). One side effect of this was that the higher µ was often associated with a more negative & J P &Q J (i.e. the difference of growth rates between productive and unproductive species is greater when growth rates are higher), because higher growth tends to induce stronger interspecific competition (Hutson 1979). In addition, the community becomes more diverse, containing more unproductive species (i.e. larger v), and thereby reducing the productivity at the community level.
In spite of the negative impacts of larger v and more negative was still larger in the high diversity treatment, although the differences were small (Fig. S6g).
At hotspot B, despite of the differences in the temporal trajectories of mean sizes due to the different size diversity, growth rates were almost identical between the two treatments ( Fig. S6j-l). Due to the higher size diversity and negative & J P &Q J values, the high diversity treatment had the slightly lower community productivity than the low diversity treatment, but the differences were negligible (Fig. S6n).
The key difference between the oligotrophic hotspot A and mesotrophic hotspot B is that in A, the sensitivity of µ to mean size was greater than in B so that faster adjusting gave greater increases in growth rate µ. These adjustments reflect a model community that has to optimize the tradeoff between achieving the highest productivity by retaining less unproductive species and the capacity to adapt to environmental fluctuation by increasing diversity. In general, the adaptive capacity is preferred when environmental fluctuation amplitudes are large (Smith et al. 2016).

S.1.6 Spatial patterns of covariances and second partial derivatives in surface waters
In Fig. S7, we compare the spatial distributions of temporal variances of ( Q O ), covariances between and light ( Qoe ), nitrate ( Qu ), and iron ( Qwxy ) in surface waters in both low-and high-diversity treatments.
Compared to the low diversity treatment, Q O increased substantially along the fronts between the central gyre and the two adjoining subarctic and equatorial Pacific. Note that because °JP°Q J QVQ E was usually negative, increases of temporal variances of mean size in the high diversity treatment actually reduced phytoplankton growth rate. The enhancement of Qwxy in the high diversity treatment was evident only in a few coastal areas. Qu also increased substantially along the fronts in the high-diversity treatment and also decreased evidently in the subarctic Pacific where phytoplankton growth is mostly limited by light. Comparatively, the enhancement of Qu was mostly pronounced in the light-limiting subarctic Pacific, while it mostly decreased in the western Pacific east of Japan.

S.1.7 Idealized numerical experiments
In the continuous case of the idealized numerical experiments, the model resolves only two compartments, nutrient (N) and phytoplankton (P).
where D is the nutrient turnover rate (d -1 ). µ com is the phytoplankton community-based growth rate (d -1 ).
As in the 3D simulations, we assume that phytoplankton size (log cell volume l, ln µm 3 ) follow a normal distribution with the mean and variance v ((ln µm 3 ) 2 ). The average growth rate of the community can then be calculated as: where is the growth rate at mean size, which depends on nutrient concentration N.
in which both µ m and K N are functions of l, similar to the 3D simulations: Nutrient inflow N 0 (t) is a seasonally varying sinusoidal function: where S 0 is the time-averaged nutrient supply. A and λ are the amplitude and period of the nutrient fluctuation, respectively.
In the simulations of investigating relationships between primary production and species richness (R), we randomly sample R species (1 ≤ R ≤ 10) with the sizes between 0.5 µm and 200 µm. The dynamics of nutrients and phytoplankton follow: where P i and µ i are the biomass and growth rate of the i th species. µ i depends on the size of the species as in Eq. (32-34).
The following parameters were used in the simulation: The initial N and total P are 1 and 0.001 µM N, respectively.  (2003)  chlorophyll-to-carbon ratios, phytoplankton mean ESD, the variance of log cell volume, growth rate at mean size, whether iron is limiting (1 indicating iron limitation and 0 indicating nitrogen limitation), and the extent of light and nutrient limitations (1 means no-limiting and 0 means 100% limiting and no growth). Note that we set that phytoplankton growth can be limited only by iron or nitrogen, whichever is more limiting. But the phytoplankton growth rate is a multiplicative function of light and nutrient terms, indicating both light and nutrient can be limiting.

S.2 Supporting
27 Fig. S5 Logarithmic ratios of NPP between the high-and low-diversity treatment (the same as Fig. 4d) overlaid with two hotspots.
28 Fig. S6 Seasonal changes of dissolved inorganic nitrogen (DIN), dissolved iron, mean size, size diversity, growth rate at mean size, second derivative of growth rate at mean size, and community productivity at two selected hotspots where negative and positive effects of diversity are most evident (see Fig. S5).

29
Fig. S7 Spatial comparison of temporal (a, b) variances of mean size and covariances (c, d) between mean size and light, (e, f) between mean size and nitrate, and (g, h) between mean size and iron within the surface mixed layer in the high-and low-diversity treatments. °JP°Q°oe QVQ E ,oeVoe E within the surface mixed layer in the high-and low-diversity treatments. X 0 represents the annual mean value of variable X.