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

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.

Introduction

Studies on the relationships between biodiversity anpan>d ecosystem functionpan>ing (BEF) have spurred much progress anpan>d debate in ecology during recent decades (Loreau et al. 2001; Hooper et al. 2012). These BEF studies provide the theoretical bapan> class="Chemical">sis for the serious concern whether losses of biodiversity may diminish ecosystem functioning (e.g. productivity) and services (Cardinale et al. 2012). In terrestrial ecology a general consensus has been reached that, even after controlling for other confounding variables such as biomass, the effect of plant diversity, mostly represented by species richness (i.e. number of species), on primary productivity is generally positive, affirming the importance of protecting biodiversity (Cardinale et al. 2006, 2007; Tilman et al. 2014; Grace et al. 2016).

Pon class="Chemical">sitive relationpan>ships between diverpan> class="Chemical">sity and productivity can arise from both ‘selection effects’ and ‘niche complementarity’ (Loreau & Hector 2001; Cardinale et al. 2006; Loreau 2010). Selection effects may enhance the productivity of more diverse communities when these have greater probability of including the most productive, i.e. best adapted, species that tend to dominate over time, compared to less diverse communities. The complementarity effect arises when niche differentiation confers complementary resource requirements to different species at different times or places, or when species coexistence is mutually beneficial via niche facilitation (e.g. symbiosis).

While the above inn class="Chemical">sights have provided the theoretical bapan> class="Chemical">sis for understanding the effects of diversity on productivity, the roles of environmental variability have not been thoroughly investigated. Theoretical studies have proposed that more diverse communities can be more productive in sufficiently variable environments, although the presence of unproductive species may reduce productivity under low variability (Norberg et al. 2001; Smith et al. 2016).

Here we investigate BEF relationship for marine phytoplankton, the most numerous photosynthetic organisms on Earth, which contribute nearly half of global primary production (Field et al. 1998). Pelagic ocean ecosystems differ fundamentally from their terrestrial counterparts. One obvious difpan> class="Chemical">ference is that the dominant oceanic primary producers are phytoplankton, mostly smaller than 200 microns in diameter. BEF studies on marine microbial organisms that directly manipulate biodiversity are particularly scarce (Hooper et al. 2005; Krause et al. 2014). Some pioneering studies have used ocean numerical models to evaluate the effects of phytoplankton functional diversity on productivity (Goebel et al. 2014; Vallina et al. 2017). Their approach, similar to most BEF experiments, is to sample randomly subsets of different numbers of species from the total species pool and seed them into the same environment. The community is then allowed to self‐organise in the model ocean, and the results tend to confirm that species richness enhances productivity, as found in terrestrial BEF studies.

However, for microbial organisms it is widely accepted that, ‘Everything is everywhere, but the envn class="Chemical">ironment selects’ (Baas‐Becking 1934; Finlay 2002; Follows & Dutkiewicz 2011). That tenet suggests that for microbes, species richnpan>ess canpan> be practically infinite, which echoes with the familiar ‘Paradox of Planpan>ktonpan>’ (Hutchinsonpan> 1961; Kashtanpan> et al. 2014). Therefore, it is more appropriate to assume a conpan>tinuous distributionpan> for phytoplanpan>ktonpan> traits, here denoted by l. The community average growth rate canpan> then be expressed in terms of the statistics of the trait distributionpan> (Wirtz & Eckhardt 1996; Norberg et al. 2001; Merico et al. 2009): Here μcom represents the per capita growth rate (d−1) of the total community, equivalent to where P is total community biomass of phytoplanpan>ktonpan>. Thus, μcom canpan> be anpan> index for productivity, which is usually correlated with the formal definitionpan> of primary productionpan> (i.e. pan> class="Chemical">organic carbon produced per unit time) (Vallina et al. 2014b). In terrestrial studies, biomass yield is often used as a proxy for productivity. However, the fast turnover of oceanic phytoplankton may decouple their productivity from biomass.

represents the mean of trait value l, which determines the per capita growth rate or n class="Disease">fitness μ. v represents the varianpan>ce of l, a proxy of diverpan> class="Chemical">sity. The second derivative of μ (), appearing in both eqn 1a and 1c can be understood as the effect of competition on both community productivity and diversity. Intense competition leads to narrow peaks of fitness around the optimal trait and reduces both productivity and diversity, known as ‘Competitive exclusion’. Eqn 1b states that the rate of change of mean trait is proportional to trait diversity, analogous to Fisher's fundamental theorem of natural selection (Fisher 1930). Eqn 1b also captures the selection effect described above: greater diversity allows the community to retain more species that differ in some functionality, thereby enhancing community productivity under environmental fluctuation, known as the ‘insurance effect’ (Yachi & Loreau 1999). Therefore, eqn 1 provides an ideal theoretical framework for investigating BEF relationships for microbial organisms having nearly continuous trait distributions. Note that eqn 1 can be easily extended to two or more traits (Wirtz & Eckhardt 1996; Savage et al. 2007). For the sake of simplicity, we herein assume that size is the only master trait for phytoplankton, because many key traits that quantify aspects of phytoplankton physiology, such as nutrient uptake and photosynthesis, vary systematically with size (Litchman & Klausmeier 2008; Finkel et al. 2010; Edwards et al. 2012, 2015; Marañón 2015). Phytoplankton size structure is also an important determinant of community respiration (Del Giorgio & Williams 2005) and the efficiency of the biological pump, i.e. carbon export from the euphotic zone (Laws et al. 2000). We assume a trade‐off between maximal growth and adaptation to low resource (nutrient or light) availability. Small phytoplankton can be considered as ‘gleaners’ adapted to oligotrophic environments, while fast‐growing intermediate‐size ‘opportunists’ thrive in resource‐rich environments (Grover 1990; Barton et al. 2010; Smith et al. 2016; Vallina et al. 2017). Therefore, we expect this tradeoff to result in a complementarity effect in temporally and spatially variable environments, leading to an overall positive effect of size diversity on productivity (Loreau 2010). This community‐based approach is also computationally advantageous compared to resolving discretely a finite number of species. It therefore allows a wider range of set‐ups and numerical experiments compared to species‐based approaches.

Here we investigate the efn class="Chemical">fects of phytoplanpan>ktonpan> pan> class="Chemical">size diversity on primary productivity in the North Pacific using the continuous trait‐based approach described above. We model the dynamics of total phytoplankton biomass, mean size (log cell volume) and size variance in a three‐dimensional ocean circulation model using realistic physical forcing for the North Pacific, covering from equatorial to subarctic regions. The North Pacific can be broadly classified into several major biogeographic provinces (Ducklow 2003; Moore et al. 2013). In the oligotrophic gyre, permanent stratification limits upward supply of nitrogen into the euphotic zone and therefore primary production remains low. In the subarctic North Pacific, characterised by low temperature and light, nitrate concentrations remain consistently high and primary production is potentially limited by iron availability. The equatorial Pacific is another high‐nitrogen‐low‐chlorophyll (HNLC) region where equatorial upwelling maintains high nitrate concentrations and production is limited by iron. Due to the growth limitations by either nitrogen, iron, or light, picophytoplankton (smaller than 3 μm in diameter) tend to dominate in the open ocean, while larger cells can become important in coastal waters (Odate 1996; Fujiki et al. 2014).

It needs to be emphan class="Chemical">sised that the complicated intertwining of diverpan> class="Chemical">sity, productivity, and the environment in nature as well as in our model ocean poses a major challenge for BEF studies (Huston 1997; Grace et al. 2016). A significant bivariate correlation between productivity and diversity does not constitute unambiguous evidence for a positive or negative BEF relationship. To generate a diversity gradient independent of environmental effects, we employ two approaches to sustain different levels of diversity in the model. The first is to vary the ‘trait diffusion’ (TD) coefficient (u), which is the probability that the offspring of individuals from one trait class evolve into other trait classes via genetic mutation or trans‐generational plasticity (Merico et al. 2014). The second is to vary the zooplankton ‘kill‐the‐winner’ (KTW) grazing coefficient (a ), which describes how the zooplankton feeding preference changes with prey abundance (Vallina et al. 2014a,b, 2017; Wirtz 2014). By choosing different values of the TD and KTW parameters respectively, to generate diversity gradients, we can indirectly separate environmental effects on phytoplankton productivity from the effects of size diversity per se (Grace et al. 2016). By combining three‐dimensional (3D) ocean modelling and idealised simulation experiments, we test the following hypotheses:

Higher diversity (induced by mutation rates or density‐dependent feeding preferences) should lead to higher primary production in general.

The effects of size diversity on productivity should depend on which environmental factor (nitrogen, iron or light) is most limiting for phytoplankton growth.

The two approaches that resolve discrete species vs. moments of the continuous distribution should generate conn class="Chemical">sistent patterns.

Materials and Methods

North Pacific model

We constructed a plankton ecosystem model within which phytoplankton size follows a conpan>tinuous distributionpan> anpan>d then coupled the ecosystem model with a three‐dimenpan> class="Chemical">sional hydrodynamic model of the North Pacific (Shchepetkin & McWilliams 2005). The ecosystem model was built on a typical nitrogen‐based, Nutrient‐Phytoplankton‐Zooplankton‐Detritus (NPZD) plankton model with the addition of an iron cycle and a lognormal distribution for phytoplankton size (cell volume) (see Fig. S1 in Supporting Information). We quantify size diversity in terms of the variance of log‐transformed cell volume, following previous studies (Wirtz 2014; Acevedo‐Trejos et al. 2015; Smith et al. 2016). The model details have been reported in Chen & Smith (2018) with the only difference that only one zooplankton compartment is included in the present study. Here we briefly describe the main features that are most relevant for the present study. The ecosystem model has eight tracers (dissolved inorganic nitrogen (N, mmol N m−3), phytoplankton (P, mmol N m−3), zooplankton (Z, mmol N m−3), detrital nitrogen (D, mmol N m−3), dissolved iron (fer, μmol m−3), detrital iron (DETFe, μmol m−3) and two raw moments of phytoplankton biomass distribution ( and where is mean log volume (μm3) and is the variance ((log μm3)2) of log volume)). Bruggeman (2009) has shown that the raw moments of the phytoplankton biomass distribution can be treated as typical tracers subject to advection and diffusion.

For a given size class, phytoplanpan>ktonpan> per capita growth rate (μ, d−1) depends onpan> temperature, light (I, W m−2), N anpan>d pan> class="Chemical">fer, following the Liebig Monod‐type function:with all the parameters including maximal growth rate (μ), half‐saturation constant for nitrogen (K ), initial slope for the photosynthesis‐irradiance curve (a ) and half‐saturation constant for iron (K ) depending on cell size: where μ0, K , K , a , a μ, βμ, a , a and a are parameters independent of size (Table S1). Here μ is assumed as a unimodal function of l, reflecting higher respiratory costs in picoplankton, which gives nanoplankton an advantage under nutrient‐replete conditions (Chen & Liu 2010; Wirtz 2011; Marañón et al. 2013). The half‐saturation constants increase with size, which favours small sizes in oligotrophic environments. Hence, this parameterisation constitutes a trade‐off between maximal growth and adaptation to low nutrient availability for marine pico‐ and nano‐phytoplankton that are the dominant primary producers in the open ocean. Compared to nutrients, the effect of light on size is weaker (i.e. a  < a ), which still confers some advantage to small phytoplankton under light limitation (Edwards et al. 2015).

Carbon based net primary production (NPP, mg C m−3 d−1) on the community level is calculated as:in which Q is the nitrogen‐to‐carbon ratio (mol N (mol C)−1) that increases with ambient N or fer (assuming iron availability limits nitrogen uptake under iron limitation; Morel 1987):where Q min and Q max are the minimal and maximal nitrogen‐to‐carbon ratios respectively.

Further model details relevant to n class="Chemical">size diverpan> class="Chemical">sity, the TD coefficient u and the KTW coefficient a , are given in the Supporting Online Materials and Chen & Smith (2018). To analyse the effect of the diversity enhancing coefficients on total NPP, we also decomposed the differences of integrated NPP among model runs using the chain rule:where is the NPP integrated over all model grids (dV stands for the volume of each grid) from 0 to 260 m over an annual cycle. means the differences of the quantity P normalized by the differences in the TD or KTW coefficient u.

To further understand the contributions of various factors affecting the lonpan>g‐term average (i.e. expectationpan>) of phytoplanpan>ktonpan> growth rate at the meanpan> pan> class="Chemical">size (eqn (2)), which is the major quantity affecting NPP (see Results), in a seasonally variable environment, we applied the technique of Taylor expansion around the seasonal mean to the second‐order (Wirtz 2000; Mandal et al. 2014):in which , N 0 and I 0 represents the annual mean phytoplankton mean size (), limiting nutrient N (nitrogen or iron), and light level (I) respectively, in each model grid. and represent the temporal variance of x and covariance between x and y (x or y represents , N or I) respectively. The first term on the right side can be treated as a constant. Because our model results show that the seasonal variations of light, nitrate and iron are relatively insensitive to the coefficients (u) of TD or KTW (data not shown), the effect of u on the long‐term expectation of phytoplankton growth rate can be simplified to:

We varied the values of TD or KTW coefficients from 0 to 0.1 to generate six diversity gradients (u or a  = 0, 0.01, 0.03, 0.05, 0.07, 0.1) so that a total of 11 simulations were run in three‐dimension. After each simulation was run for 10 years, the total annual NPP integrated throughout the euphotic zone (from surface to 260 m) over the entire model domain was calculated for the final year.

Idealised simulation experiments

We conducted idealised simulationpan> experiments to compare our conpan>tinuous pan> class="Chemical">size‐based approach with the conventional approach of controlling the level of species richness, based on a simple phytoplankton‐nutrient model that has often been applied to investigate the effects of environmental variability on species coexistence (Supporting Online Material; Huisman 2010; Loreau 2010). For both models, the external nutrient supply varied as a seasonal sinusoidal function. Two amplitudes of the sinusoidal function were set up to simulate different levels of environmental variability. For the continuous model, we fixed the size diversity at different levels and, for each diversity level, took a random value between −2.73 log μm3 and 15.2 log μm3 (equivalent to 0.5 and 200 μm in diameter) as the initial mean size (log cell volume) of the community, running a total of 50 replicate simulations at each diversity level. In the case of controlling species richness (n, 1 ≤ n ≤ 10) with the discrete model, at each richness level, we sampled n species with sizes randomly distributed between −2.73 log μm3 and 15.2 log μm3 and ran 50 replicates. Each model configuration was run for 10 years and the annual NPP was calculated for the final year.

Results

The model was able to reproduce the large‐scale patterns of nitrate, chlorophyll a and NPP (Fig. S2). Total NPP over our model domain was estimated between 13.0 and 13.3 Pg year−1 (Fig. 1), roughly one‐fourth of annual NPP for the global ocean (Field et al. 1998). Total NPP was inversely related to u and a , although the variations were negligible. Total NPP was 2% higher in the lowest diversity treatment (u = 0 and a  = 0) compared with the highest diversity treatments (u = 0.1 or a  = 0.1).

Figure 1

Modelled annual primary production integrated for the euphotic zone (0–260 m) summed over the whole model domain against the coefficients of ‘trait diffusion’ or ‘kill‐the‐winner’.

Modelled annual primary production integrated for the euphotic zone (0–260 m) summed over the whole model domain against the coefficients of ‘trait diffun class="Chemical">sion’ or ‘kill‐the‐winner’.

The spatial patterns of annual mean n class="Chemical">NPP conpan>firmed the negligible difpan> class="Chemical">ferences between the highest (u = 0.1 or a  = 0.1) and the lowest (u = a  = 0) diversity treatments, although, as expected, size diversity was substantially lower in the lowest diversity treatment than in the highest diversity treatment (Fig. 2). The mean phytoplankton sizes and nitrogen biomass were similar among all treatments.

Figure 2

Modelled annual mean patterns of size diversity, mean size and primary production (NPP) for three diversity treatments. The first row is the lowest diversity treatment without any TD or KTW. The second and third rows are the treatments with the largest TD and KTW coefficients.

Compared to the lowest diversity treatment, the highest diverpan> class="Chemical">sity treatment agreed better with observed size‐fractionated chlorophyll patterns in terms of picophytoplankton (< 2 μm) and size variance, while both treatments similarly matched mean size data (Fig. 3). With no TD or KTW to sustain diversity, the model severely underestimated size variance (Fig. 3c).

Figure 3

Simulated patterns of (a,b) fractions of picophytoplankton (< 2 μm), (c,d) mean size and (e,f) size variance of high (u = 0.1) and low diversity treatments. Colour contours indicate model data density. Red dots indicate observational data from Marañón et al. (2012).

To understand why enhancing size diverpan> class="Chemical">sity by TD or KTW has negligible effects on NPP, we calculated and integrated the components in eqn (6) for the lowest and highest diversity treatments respectively (Table 1). The three most important quantities were the differences of growth rate , size diversity and nitrogen‐to‐carbon ratios (Q N) at mean (i.e., dominant) size respectively. However, the positive effect of size diversity on growth induced by high mutation rate (u) was overall counterbalanced by the negative effect of including more unproductive species (i.e. higher v), because was usually negative at mean size. The third factor, Q N, eventually determined the negative, albeit small, effect of size diversity on NPP. This is because high Q N was usually associated with high growth rate , reducing organic carbon production per unit nitrogen biomass. That is, although higher diversity tended to enhance , it also enhanced Q , which counteracted the effect of diversity on NPP.

Table 1

Decompositions of contributions of different components (unit: Pmol C year−1 d (ln μm3)2) in eqn (6)

∫PQ_NΔμΔudVdt	4.6
−∫PμQ_N^2ΔQ_NΔudVdt	−2.8
∫μQ_NΔPΔudVdt	−0.12
∫vPΔ∂^2μQ_N∂l^2ΔudVdt2	−0.07
∫v∂^2μQ_N∂l^2ΔPΔudVdt2	−0.03
∫P∂^2μQ_N∂l^2ΔvΔudVdt2	−4.6

Decompon class="Chemical">sitions of contributions of difpan> class="Chemical">ferent components (unit: Pmol C year−1 d (ln μm3)2) in eqn (6)

We examined the spatial patterns of the effects of pan> class="Disease">TD coefficient u on phytoplankton growth rates and NPP by calculating the logarithmic ratios of the annual , carbon‐based growth rate , community average carbon‐based growth rate and NPP, each from the highest (u = 0.1) to the lowest diversity treatment (u = 0) and integrated throughout the euphotic zone (Fig. 4). The positive effect of size diversity on , which prevails over most regions, is strongest along the fronts between the central oligotrophic regions and the two adjoining subarctic and equatorial regions (Fig. 4a). In regions where light and iron limitations are more important (Fig. S4), size diversity only weakly affects . The effect of Q N weakens the net effect of size diversity on carbon‐based growth rate (Fig. 4b). Further considering the effect of unproductive species within the community (i.e. the effect of v), we found that the positive effect of size diversity on community‐integrated carbon‐based growth rate was reduced more. However, in the subarctic and equatorial Pacific, where light limitation was strongest (Fig. S4), this decrease was absent or even reversed (Fig. 4c). Because the summed negative effects overweighed the summed positive effects, the net effect of size diversity on NPP was slightly negative for the whole North Pacific (Fig. 4d).

Figure 4

Spatial patterns of log ratios of annual mean (a) phytoplankton growth rate at mean size, (b) carbon‐based growth rate at mean size, (c) carbon‐based phytoplankton growth rate of the whole community and (d) carbon‐based primary production (NPP) integrated from surface to 260 m in the diversity treatments with the highest vs. the lowest coefficients of trait diffusion.

Spatial patterns of log ratios of annual mean (a) phytoplankton growth rate at mean size, (b) pan> class="Chemical">carbon‐based growth rate at mean size, (c) carbon‐based phytoplankton growth rate of the whole community and (d) carbon‐based primary production (NPP) integrated from surface to 260 m in the diversity treatments with the highest vs. the lowest coefficients of trait diffusion.

To better understand the mechanisms of diversity efpan> class="Chemical">fects on phytoplankton growth rate , we plotted the spatial patterns of the three components of the right side of eqn (8) within the surface mixed layer, as their difference between the high and low diversity treatments (Fig. 5). The contributions of to differences were negative along the fronts separating the central gyre and the adjoining north and south areas, while the contributions of (covariance between mean size and nutrient) were positive, corresponding to the coefficients of variations of surface nitrate. The effects of (covariance between mean size and light) were relatively minor.

Figure 5

Spatial patterns of (a,b,c) differences of different components of R.H.S. of eqn (8) including seasonal variances of mean size and covariances between mean size and light/nutrient in surface waters between the high‐ and low‐diversity treatments. (d) Coefficients of variation (= SD/Mean) of surface nitrate concentrations.

Spatial patterns of (a,b,c) difn class="Chemical">ferences of difpan> class="Chemical">ferent components of R.H.S. of eqn (8) including seasonal variances of mean size and covariances between mean size and light/nutrient in surface waters between the high‐ and low‐diversity treatments. (d) Coefficients of variation (= SD/Mean) of surface nitrate concentrations.

We ran a series of idealised n class="Chemical">simulationpan> experiments to better understanpan>d the difpan> class="Chemical">ferences between our continuous trait‐based approach, in which size diversity is the key diversity metric and the more typical approach that takes species richness as the diversity metric. For both levels of environmental variability, we found that the median NPP increased from the lowest to an intermediate level of size diversity (0.05 (log μm3)2) and then slightly decreased with increasing size diversity afterwards (Fig. 6a). When we used the more typical approach of randomly sampling a given number of species, we obtained the familiar increasing trend of NPP with species richness, for which NPP reached a plateau beyond the richness of 4 (Fig. 6b; Goebel et al. 2014; Vallina et al. 2017). For both models, at low diversity levels, NPP values are higher under low environmental variability.

Figure 6

Simulated relationships between primary production and diversity using (a) continuous trait‐based approach that assumes a lognormal distribution for phytoplankton size and (b) discrete species approach. Imposed box‐and‐whisker plots show the median (thick lines), 25th and 75th percentiles (box edges), data ranges and outliers (open circles).

n class="Chemical">Simulated relationpan>ships between primary productionpan> anpan>d diverpan> class="Chemical">sity using (a) continuous trait‐based approach that assumes a lognormal distribution for phytoplankton size and (b) discrete species approach. Imposed box‐and‐whisker plots show the median (thick lines), 25th and 75th percentiles (box edges), data ranges and outliers (open circles).

Discussion

Our result that increan class="Chemical">sing pan> class="Chemical">size diversity has negligible or even negative effect on NPP across the whole North Pacific seems to contradict many previous reports of positive effects of species richness on ecosystem productivity (Cardinale et al. 2006, 2007; Tilman et al. 2014), including studies using ocean circulation models (Goebel et al. 2014; Vallina et al. 2017). This is somehow unexpected given the already well described theoretical mechanisms including both ‘selection effects’ and ‘niche complementarity’ and experimental evidences for the positive effect of diversity on productivity (Hooper et al. 2005; Cardinale et al. 2012; Tilman et al. 2014). However, as argued in Hillebrand & Matthiessen (2009), the understanding of diversity effects on productivity would be superficial if a causal linkage between individual functional traits and ecosystem functioning is not adequately set up. We propose that eqns (6), (7), (8) can be used as a quantitative framework to understand diversity effects on productivity in variable environments and then the specifics conditions for negligible or negative diversity effects on productivity to occur can be inferred.

Covariance between trait and environment

Equation 1a reveals that the positive efpan> class="Chemical">fect of diversity on productivity at the community level must counteract the generally negative effect of diversity on the growth rate at the mean (best‐adapted) trait, because the second derivative of growth rate should be negative under most conditions. In the absence of positive species interactions or niche differentiation under static conditions as in our case, diversity can enhance productivity exclusively in a variable environment. More specifically, based on eqn (8), only the covariance between the trait and the environment can contribute to positive effects of diversity on productivity. The positive effect of size diversity on productivity under environmental variability involves a selection effect. Compared to a low diversity community, in a high diversity community the dominant trait class tracks the environment faster, i.e. approaches closer to the optimal trait (Fig. S6). Favourable environments (high nutrient or light) select fast growing opportunists of large size (within the phytoplankton size spectrum relevant for the North Pacific), while unfavourable environments select small gleaners that are relatively insensitive to resource shortage. This temporal or spatial niche partitioning may also be viewed as a result of spatiotemporal complementarity increasing overall production, analogous to the ‘storage effect’ induced by the covariance between the environment and competition, which helps to maintain diversity in variable environments (Chesson 2000).

Importance of trait characteristics

Equation (8) makes clear that different functionpan>al growth dependencies, onpan> both traits anpan>d envpan> class="Chemical">ironmental conditions, can have different effects on BEF relationships. For example, light can affect both phytoplankton size and productivity. The covariance between light and mean size () turns out to be less important than because the combined growth dependency is small. Therefore, despite sometimes large differences of between different diversity treatments, the overall light effect on productivity is relatively small compared to the nutrient effects (Figs S7 and S8). This is because most production takes place within the surface mixed layer with relatively high light but low nutrient availability. Furthermore, the optimal size class at low nutrient levels is similar to that under light limitation, so that moderate changes in light levels hardly affect the size distribution.

Loreau (2010) emphan class="Chemical">sised that the asynchronpan>ous responpan>ses of species in variable envpan> class="Chemical">ironments are the key to temporal complementarity, which generates a positive diversity effect on productivity. By contrast, synchronised responses can be understood as niche overlapping. This implies that niche differentiation can promote diversity maintenance and positive diversity effects on productivity, particularly under environmental variability. Vallina et al. (2017) pointed out that the complementarity effect may be negligible for niches with open‐ended forms such as nutrient uptake; i.e. nutrient uptake is not inhibited at high concentrations. This is because even in the absence of the most productive species, a less efficient species can still occupy its niche (utilize the nutrients) at a rate only slightly slower. In the present study, both nutrient uptake and light acquisition are open‐ended functions with positive first derivatives and negative second derivatives. Species responses to both nutrient or light are more synchronised at high resource levels, typical of light, than at low resource levels, typical of nutrient, making nutrient more important than light for determining diversity‐productivity relationships. Much stronger complementarity effects are expected for closed‐ended niches (e.g. unimodal functions) such as optimal temperature or light, which impose greater disadvantages on sub‐optimal species (Vallina et al. 2017). That is, on a continuous trait space, fitness decreases more steeply away from the optimal trait for a closed‐ended niche than for an open‐ended one (fig. 2 in Vallina et al. 2017). Hence for a closed‐ended niche in a dynamic environment, greater trait variance/diversity v is more likely to enhance fitness (up to a point), by increasing the speed at which the mean trait (eqn 1b) tracks its optimal value, resulting in stronger covariance between trait and environment.

The master trait cell size mostly correlates with nutrient uptake strategies (Litchmanpan> et al. 2007; Edwards et al. 2012; Marañón et al. 2013). Thus, pan> class="Chemical">similar to nutrient uptake traits, the lack of a strong complementarity effect of size diversity on productivity is therefore not surprising. For future studies it would certainly be worthwhile to include more traits, such as optimal growth temperature and light that are not strongly dependent on size, particularly in areas with substantial seasonal or spatial fluctuations of temperature or light (Goebel et al. 2014; Vallina et al. 2017; Bestion et al. 2018).

Environmental variability and negative effect of trait diversity on productivity

n class="Chemical">Since early BEF studies, it has already been realised that at high diverpan> class="Chemical">sity levels, increasing diversity can reduce productivity particularly when the environment is not sufficiently dynamic, which can lead to a unimodal relationship between productivity and diversity (Hector et al. 1999; Norberg et al. 2001; Hillebrand & Matthiessen 2009; Tilman et al. 2014; Smith et al. 2016). As such, it is not surprising to find negative effects of size diversity on community productivity, which usually occur in areas with negligible and surface nitrate varies little (Figs 4 and 5). Conversely, in dynamic environments, high diversity is needed to achieve selection effects sufficiently strong to generate positive covariance between the mean size (trait) and nitrate (environment) and thereby to enhance productivity (Figs 5 and 6).

Saturating diversity levels for productivity

Schwartz et al. (2000) pointed out that productivity may saturate at quite low diversity levels, less thanpan> half of total species richnpan>ess. Our idealised pan> class="Chemical">simulation experiments showed that the positive effect of size diversity on productivity is evident only at very low levels of size diversity (< 0.05 (ln μm3)2), above which this effect becomes negligible (Fig. 5). This effect is also evident in the discrete case, in which productivity saturates around richness of 4. The low level of size diversity that saturates productivity also relates to the trait characteristics (e.g. the open‐ended niche form) discussed above.

In the 3D model ocean, even in the unrealistic cases without n class="Disease">TD or KTW to sustain diverpan> class="Chemical">sity, advection and diffusion still maintain levels of size diversity greater than 0.05 (ln μm3)2 in many areas (Figs 2 and 3). Vertically, deepening mixed layers during fall can entrain communities living at depth where nutrients are plentiful but light is limiting, into surface waters where nutrient can be limiting (Chen & Smith 2018). Horizontally, high diversity can also emerge along ocean fronts where mixing is active (Barton et al. 2010). Thus, three‐dimensional water mixing can be an important mechanism to sustain diversity and consequently relatively high productivity. However, some additional mechanism such as TD or KTW is required to reproduce the observed levels of phytoplankton size diversity based on size‐fractionated chlorophyll measurements (Fig. 3). Therefore, we expect that within realistic ranges of TD or KTW parameters, the weak negative effect of size diversity on productivity should be robust.

Stoichiometry effects

Another important factor contributing to the reduction of integrated n class="Chemical">carbon‐based pan> class="Chemical">NPP with increasing diversity is the variable nitrogen‐to‐carbon ratio of phytoplankton, which offsets the enhancement of growth rate at high diversity (Table 1). This underlines the importance of a clear definition of productivity (per capita growth rate) as distinct from absolute rates of production. Our results are consistent with the view that the emergent pattern of phytoplankton stoichiometry can result from diversity and, vice versa, stoichiometry can also affect competition and diversity patterns among phytoplankton (Bonachela et al. 2015). Classical competition theory has rarely accounted for flexible stoichiometry. Future studies are needed to extend the size‐based approach to multiple traits including flexible stoichiometry. This will be complicated, because minimal and possibly also maximal nutrient ratios correlate with size, while the actual C : N : P : Si : Fe : Chl ratios mostly reflect the environmental history of cells. Furthermore, diversity needs to be incorporated into a coherent and robust modelling framework that can be practically useful for assessing how global change and anthropogenic activities affect biodiversity and ecosystem functioning in the ocean.

Concluding remarks

While our first assessment reveals a negligible effect of phytoplanpan>ktonpan> pan> class="Chemical">size diversity on productivity, the mechanisms underpinning BEF relationships should be common to all ecosystems. Particularly, our analysis suggests the importance of interactions between environmental variability and trait distributions. The environment can encompass several dimensions, such as light, temperature and different nutrients in the case of autotrophs and food concentrations and qualities in the case of heterotrophs. The effective trait space is therefore multi‐dimensional, with each dimension relating to one or more environmental axes (Savage et al. 2007). The details of each function relating traits and the environment may have far‐reaching implications for the intricate interactions among the environment, biodiversity and ecosystem functioning. Thus, it is expected that the overall diversity effect on productivity depends on the spatial scale considered, with higher probability of positive effect of diversity on productivity to be observed in environmentally more dynamic regions such as ocean fronts. While our exercise provides an initial step, better mathematical and modelling tools are needed to disentangle such complexity. On the experimental side, marine experiments similar to the terrestrial Cedar Creek and BIODEPTH experiments that examined the effects of diversity per se on productivity (Tilman et al. 1997; Hector et al. 1999) are also much needed to validate the results of theory and numerical experiments.

Authorship

BC and SLS designed the study. BC wrote the model codes, ran the simulations, and wrote the first draft of the manuscript. All authors discussed the results and contributed to revisions.

Conflict of Interest

The authors declare that they have no conflict of interest.

Click here for additional data file.