How important is diversity for capturing environmental-change responses in ecosystem models?

. Marine ecosystem models used to investigate how global change affects ocean ecosystems and their functioning typically omit pelagic plankton diversity. Diversity, however, may affect functions such as primary production and their sensitivity to environmental changes. Here we use a global ocean ecosystem model that explicitly resolves phytoplankton diversity by deﬁning subtypes within four phytoplankton functional types (PFTs). We investigate the model’s ability to capture diversity effects on primary production under environmental change. An idealized scenario with a sudden reduction in vertical mixing causes diversity and primary-production changes that turn out to be largely independent of the number of coexisting phytoplankton subtypes. The way diversity is represented in the model provides a small number of niches with respect to nutrient use in accordance with the PFTs deﬁned in the model. Increasing the number of phytoplankton subtypes increases the resolution within the niches. Diversity effects such as niche complementarity operate between,


Introduction
Ocean ecosystems are under pressure from global environmental change and an increasing human demand for natural resources. As a consequence of often ultimately anthropogenic perturbations, a loss of diversity has been observed in a variety of different ecosystems including marine environments (Butchart et al., 2010). Increasing evidence suggests that such a diversity loss coincides with a reduction in ecosystem functioning such as primary production or nutrient use (Cardinale et al., 2011;Reich et al., 2012). Losses in functioning are small in highly diverse systems and tend to increase with decreasing diversity. Diversity might thus potentially aid in sustaining an ecosystem's established functioning through periods of environmental change.
Experimental evidence for effects of diversity on ecosystem functioning predominantly originates from terrestrial and benthic ecosystems, while studies on pelagic communities are still scarce (Cardinale et al., 2011;Ptacnik et al., 2010). Here we use an ocean ecosystem general circulation model to approach investigating diversity effects in the pelagic ocean on the global scale. There has been a long debate on the required level of complexity within ecosystem models (Ward et al., 2013). Friedrichs et al. (2007) found that more complex models ported more successfully between different regions than simpler models. However, an earlier paper (Friedrichs et al., 2006) also suggested that the difference in responses between ecosystem models was often dwarfed by the responses between different physics. However, an intercomparison of 10 ecosystem models to 21st century climate change (Bopp et al., 2013) showed very different regional responses in primary production even though they all showed regional warming trends that were almost everywhere robust among the different models. These models only agreed on the sign of primary production change in very few locations. How much the change was related to the difference in complexity of the ecosystem models, and how much to different assumptions about the biogeochemical and grazing environments, is not clear.
The majority of such studies of marine ecosystem changes with an altered climate have employed at most a handful of phytoplankton functional types (PFTs). Here we instead employ a novel global coupled ocean ecosystem modelling approach  which resolves phytoplankton diversity within four PFTs and allows investigating some aspects of diversity effects on ecosystem functioning without observational and experimental limitations. Plankton diversity in this model is represented by variability in nutrient-, light-and temperature-dependent phytoplankton growth. This approach clearly reduces natural plankton diversity to a limited number of traits in the model. Models incorporating other aspects of diversity, however, such as algal mixotrophy (Ward et al., 2011) or spectral light use (Hickman et al., 2010), have not yet been applied to the global scale. Our model is thus representative of global biogeochemical ocean circulation models currently used for investigating biogeochemical fluxes on large spatial and temporal scales (e.g. Bopp et al., 2013).
We examine whether the magnitude of changes in primary production arising from environmental change depends on the level of phytoplankton diversity as indicated by the number of phytoplankton subtypes and PFTs in the model. We relate our findings to the underlying structure of ecological niches implicit in the model formulation, and thereby assess the model's inherent ability to capture diversity effects such as niche complementarity (Tilman et al., 1997;Loreau, 1998) and selection effects (Aarssen, 1997;Huston, 1997). Niches in the classical sense are created by variability in ecological factors and can be identified by the species fitness as a function of the magnitude of the respective ecological factor (Hutchinson, 1957;MacArthur and Wilson, 1967;Schoener, 1988). In plankton communities, in addition to this environmental dimensionality the traits of the species -for example regarding resource uptake, tolerance width to environmental conditions or mobility -shape further niches, thereby creating trait dimensionality (Ptacnik et al., 2010). Within this framework, diversity effects are thus constrained by both environmental and trait dimensionality.
One of the predicted consequences of global change on the upper ocean is an increased stratification and an associated reduction in mixing processes, and thus nutrient supply to the surface mixed layer (e.g. Sarmiento et al., 2004). Here we impose an instantaneous reduction in vertical mixing for nutrients and plankton as idealized environmental change to examine its consequences for diversity and productivity in a model ocean. We employ simulations with different levels of prescribed phytoplankton diversity and different numbers of PFTs in order to examine if diversity and the way to resolve it affects the simulated primary production changes.

Methods
The model employed is the Darwin ocean ecosystem model coupled to the MITgcm general circulation model  in the configuration used by Prowe et al. (2012a). In the standard setup, the model simulates the dynamics of four nutrients (phosphorus, nitrogen, iron, silica), 78 phytoplankton subtypes, two zooplankton functional types, and dissolved and particulate organic matter. Phytoplankton subtypes are assigned to one of four PFTs (large diatom and small Prochlorococcus analogues, other large and other small phytoplankton). Cell size determines sinking speed, palatability, maximum growth rate and the basic level of the halfsaturation concentration for nutrient uptake of the different PFTs. Within each PFT, subtypes are further distinguished by randomly assigned parameter values for nutrient-, light-and temperature-dependent growth. Half-saturation constants for phosphorus uptake are randomly assigned within a given range around a basic value characteristic of the PFT. Corresponding half-saturation concentrations are adopted for the other nutrients using the same fixed stoichiometric ratios for all phytoplankton subtypes. Temperature dependencies of individual subtypes are characterized by different optimal temperatures chosen at random from a range of −2 to 30 • C. Grazing on phytoplankton by a small and a large zooplankton functional type is formulated as a Holling type 3 functional response. Details can be found in Dutkiewicz et al. (2009) and Prowe et al. (2012a).
Phytoplankton diversity or richness is measured at each grid point at each time step as the number of phytoplankton subtypes exceeding a low threshold concentration of P th = 10 −8 mmol P m −3 . As the standard simulation we select one member of the ensemble of simulations with n = 78 initial phytoplankton subtypes used by Prowe et al. (2012a). Three simulations with reduced initial diversity are obtained by randomly selecting subpopulations of n = 30 subtypes from this setup. A configuration with each PFT represented by only one phytoplankton subtype (n = 4) with the optimumtemperature function replaced with a simpler temperature dependence (represented by a factor of 1.04 T [ • C] ) allows simulating PFT-based ocean ecosystem models that do not resolve diversity within PFTs, comparable to models currently employed for global change simulations (e.g. Bopp et al., 2013). This so-called Eppley temperature function integrates the diversity in temperature responses resolved explicitly in the simulations with many subtypes. Although previous studies suggest that capturing individual temperature responses might significantly influence model behaviour (e.g. Moisan et al., 2002), preliminary studies indicate that PFT distribution is not influenced significantly in this model. The use of the Eppley curve thus allows us to focus our attention on the nutrient uptake aspect rather than on environmental temperature changes (see below). Two companion simulations -one without the Prochlorococcus and one without the other-small PFT (n = 3p and n = 3o, respectively) -reveal consequences of reducing PFT diversity, which indicates functional diversity in the model. An overview of the different configurations is shown in Table 1.
Simulations are run offline for 20 yr with physical forcing from the ECCO-GODAE state estimates (Wunsch and Heimbach, 2007). Environmental changes are prescribed as a sudden reduction in mixing (vertical eddy diffusivity k e reduced by 50 %) after 10 yr of simulation. The change in k e in the offline model does not affect temperature, stratification or the depth of the mixed layer, but only the mixing of nutrients, phytoplankton and the other biogeochemical tracers. Excluding physical ocean changes allows us to separate direct effects on the ecosystem from more complex responses driven by physical feedbacks. This idealized setup also allows us to simulate diverse phytoplankton communities explicitly without compromising spatial or temporal resolution. Annual average diversity and primary production in the upper 100 m are compared to simulations with no reduction in mixing after another 10 yr of simulation. Primary production is estimated from phosphate uptake using a fixed molar C : P ratio of 106 : 1.

Results
The reduction in k e leads to an overall decrease in nutrient supply, and thus in primary production (PP) and net community production in the lower latitudes and increases in the Southern Ocean and the North Atlantic where PP is limited by light during large parts of the seasonal cycle ( Fig. 1). Figure 2 shows the corresponding absolute and relative differences of diversity and PP in year 20 between the simulations with reduced k e and the respective simulations with unchanged k e . For the standard setup with n = 78 phytoplankton subtypes, diversity measured as the number of coexisting subtypes overall changes by less than 5-10 % upon reducing k e , with increases at higher latitudes and decreases from roughly 40 • S to 40 • N. Relative PP changes are larger, with reductions around 10 % in the tropics and subtropics and reductions exceeding 20 % around 30 • S and 30 • N where PP is lowest. At higher latitudes, PP changes are generally small and, in the zonal average, positive only in the Southern Ocean.
Higher PP under global warming in this region has been inferred both from observations (Behrenfeld et al., 2006) and from coupled ocean ecosystem models (e.g. Bopp et al., 2001). In these studies, it can be explained by higher light availability for photosynthesis in a shoaling surface mixed  layer. The effect of higher light availability at higher latitudes is reflected at least partly by lower phytoplankton losses due to mixing in the short-term simulations presented here. This increases PP (Fig. 1) and also the potential for export of organic matter to the deep ocean, particularly in regions that do not become nutrient depleted during the seasonal cycle. A PP decrease in the remaining ocean reflects aggravated nutrient limitation and closely corresponds to a decline in export production.
In the model, PP (= i µ i P i ) may change because of changes in the specific growth rate (µ i ) and/or the biomass (P i ) of the individual phytoplankton subtypes i. Growth rate changes reflect the direct effect of changed nutrient concentrations or changed light conditions for each individual subtype. Biomass changes in turn reflect both changes of total biomass across all phytoplankton subtypes as well as shifts in the community composition. Biomass changes indicate a shift in the balance between gains and losses, whereas changes in community composition point towards a shift in competition between subtypes. In the oligotrophic low latitudes (about 30 • S to 30 • N; here shown for an Atlantic section), PP changes are almost entirely driven by a reduction in total phytoplankton biomass (Fig. 3). Community composition shifts and changes in the individual growth rates have only negligible effects. In these regions, the reduced supply of nutrients from deeper layers into the productive surface layer decreases gains in relation to losses for all phytoplankton subtypes alike without large changes in individual growth  Fig. 2e) due to the reduction in mixing, and decomposition into effects of total biomass, growth rate, or composition changes along 25 • W. The mixing reduction affects the specific growth rate (red) via nutrient and light conditions as well as the individual biomass of phytoplankton subtypes via effects on both total biomass and community composition (solid blue). Effects on the individual biomass arise from changes in total biomass (dashed light blue) and shifts in community composition, i.e. the relative biomass of each phytoplankton subtype (dashed dark blue). Displayed are zonal averages from the simulation with n = 78 subtypes. See the appendix for details on the decomposition.  Fig. 3. The lower growth rate is more than compensated by an overall increase in biomass. We find no pronounced difference in the response of simulated PP and diversity to a reduction of k e between the simulations with n = 78 subtypes and the three simulations with n = 30 subtypes. Although the former simulation has overall higher PP at higher latitudes (cf. Prowe et al., 2012b), relative changes in both properties are of similar order of magnitude and display similar zonal patterns among the simulations. Global average PP changes vary by only 1 % among the model runs employing different levels of diversity. Differences within the ensemble of three simulations with similar initial (n = 30) diversity are comparable with differences between simulations with different initial diversity (Fig. 2). "Diversity" in the generic simulation with each PFT represented by just one subtype (n = 4) is affected by mixing changes only in the Southern Ocean. In all other regions all PFTs are present (i.e. they exceed the biomass threshold for diversity of P th = 10 −8 mmol P m −3 ), although typically one or two PFTs dominate a region. Changes in PP in response to reduced mixing are generally comparable to changes predicted for the simulations resolving diversity within the PFTs. Southern Ocean PP increases more strongly in the simulation with low diversity (Fig. 2e, f). Here, surface nutrients tend to be higher in the low-diversity simulation compared to the more efficient nutrient utilization in the high-diversity runs. This leaves more potential for a PP increase upon a mixing reduction in the lower-diversity simulations. In the oligotrophic lower latitudes, PP decreases more strongly in the low-diversity runs, although the essentially always complete drawdown of surface nutrients is not affected by the reduction in mixing. In contrast to the simulation with four PFTs, the n = 3p simulation with three PFTs and Prochlorococcus omitted shows higher relative PP changes between 20 and 50 • N, i.e. in regions where Prochlorococcus dominates in the other simulations.
The fact that all model communities with the same number of PFTs react in a similar way to the mixing reduction independent of diversity within the PFTs reflects the dimension of the trait space resolved by the assemblages, which is relevant in these reduced-nutrient-supply simulations. The biogeography of each subtype is limited by the optimum temperature for growth, which defines suitable habitats on the global scale independent of PFT assignment. However, the half-saturation concentrations for nutrient uptake (e.g. K PO 4 for phosphate) are assigned randomly within each PFT so that lowest values distinguish Prochlorococcus analogues, slightly higher values characterize the other-small PFT, and both large PFTs have high half-saturation concentrations. This pattern is the basis for PFTs in the generic n = 4 simulation, and is generally maintained in all other simulations (Fig. 4). The relative biomass fraction indicates that each assemblage fills only two main niches in terms of nutrient uptake traits with one or few dominant subtypes, with low-K PO 4 subtypes at high temperatures and low light (deep chlorophyll maxima in oligotrophic low latitudes) and high-K PO 4 subtypes mostly at lower temperatures and high light (spring-bloom conditions in higher latitudes). This demonstrates that the warm, oligotrophic regions select for subtypes with lowest K PO 4 ("gleaners"; Dutkiewicz et al., 2009) via resource competition, while in the colder seasonally mixed oceans high growth rates identify dominant "opportunist" subtypes. The number of subtypes per bin and fraction of total biomass emphasizes that all assemblages can by model design only resolve a few nutrient niches with similar characteristics given by the PFT definition (Fig. 5). In the low-K PO 4 niche, the assemblage is strongly skewed towards subtypes with lowest K PO 4 , while in the high-K PO 4 niche such a shift is not found.
Contrary to the n = 4 simulation, omitting the Prochlorococcus PFT, which has lowest K PO 4 levels, in the n = 3p simulation reduces the trait space resolved in the model by lifting its lower boundary for nutrient uptake. This reduction is most relevant in the oligotrophic regions (Fig. 2), which now are populated by other subtypes with higher K PO 4 . In contrast, omitting the PFT representing other small phytoplankton (n = 3o), which has intermediate K PO 4 levels, does not affect the lower boundary of the trait space and agrees well with both the n = 4 simulation and the simulations resolving diversity within PFTs. Consequently, a reduction in nutrient supply caused by reduced mixing affects PP changes more strongly in the n = 3p run compared to the other simulations.

Discussion
Analyses of observations suggest a decline in ecosystem functions resulting from potential diversity losses and a stronger decline at lower overall diversity (Cardinale et al., 2011). Our model, in apparent contrast and for a relatively small diversity range tested so far, does not simulate such an effect of the diversity level on the magnitude of productivity changes under reduced mixing. Specifically, resolving phytoplankton diversity within plankton functional types in the way done in this model does not affect the sensitivity of PP to environmental change, although explicitly resolving diversity can improve the representation of community structure and seasonal succession (Prowe et al., 2012a). As models currently employed for global change simulations often use setups similar to our simulation with four PFTs (e.g. Bopp et al., 2013), our results indicate that omitting diversity within PFTs as captured in our approach does not appear to compromise such predictions.   30 types (b-d). Colour shading indicates the fraction of total biomass for all subtypes within a K PO 4 bin. Vertical lines indicate the PFTs in the n = 4 and n = 3 simulations (see Fig. 4 for details).
The model setup employed here provides two main niches for nutrient uptake, which agree with the characteristics of the two small and large PFTs: one low-nutrient hightemperature niche and one high-nutrient niche mostly realized at lower temperatures. These niches are resolved by different phytoplankton subtypes characterized by a random combination of trait values for temperature, light and nutrient use within the boundaries of the niches according to the assigned PFT. The Holling type 3 grazing will facilitate diversity within the two niches (Prowe et al., 2012b). Additional sensitivity experiments have shown that results are qualitatively similar for simulations employing a Holling type 2 grazing functional response which promotes competitive exclusion. This emphasizes that more diverse assemblages do not occupy more niches, because there are only two main global niches defined by nutrient use resolved in the model.
On regional scales, silicate requirements for diatoms and nitrate availability for other small phytoplankton (Prochlorococcus in these simulations are assumed to not be able to assimilate nitrate) can generate additional niches. Temperature creates geographical niches as the optimum temperature of each subtype determines whether positive net growth is possible or not under a given temperature regime. In any given region, however, the temperature traits do not create complementary niches in terms of nutrient uptake, as they do not involve a trade-off in our model. Similarly, the Holling type 3 grazing functional response does not shape complementary niches because the susceptibility to grazing is, in our current model, not linked to any costs, for example related to grazing defences. Consequently, one should pose the question of what actually defines different niches in an ocean ecosystem model and in the real world, and how these relate to the PFTs and trade-offs resolved in the model.
Reducing PFT diversity may correspond to a reduction in the relevant trait space resolved and demonstrates that in current ocean ecosystem models the variety of PFTs determines the niches resolved. In the framework of Ptacnik et al. (2010), here it is the trait dimensionality resolved in the phytoplankton community that governs diversity effects on PP against the background of environmental dimensionality simulated by the model. Higher trait dimensionality, for example related to cell size (Ward et al., 2012), might potentially allow larger diversity effects, if the traits give rise to coexistence or competition along an axis of environmental dimensionality. In any given location, i.e. in a given physical and ecological environment, this condition would require model formulations based on trade-offs between different traits. Furthermore, while niches may change dominance under changing environmental conditions, the rigid structure of this modelling approach does not allow for new niches to be populated.
Small diversity effects on PP changes are evident from the differences between simulations with n = 30 subtypes, and are related to the identity of subtypes present. The n = 30c simulation reacts to the mixing reduction in a way more sim-ilar to the standard simulation with n = 78 subtypes than to the other simulations with same nominal diversity. This indicates that in this model type identity plays a more important role than overall diversity and raises the question of whether the number of coexisting phytoplankton subtypes is the most appropriate measure for relating diversity to ecosystem functions. While there is a higher probability in more diverse assemblages that a specifically suited type is present (Aarssen, 1997;Huston, 1997), our experimental setup is not aimed at identifying this effect.
The similarity of responses among the simulations with different diversity also reflects the fact that the diversity of the phytoplankton community is formulated in terms of a small number of traits that might not include the most relevant traits under changing environmental conditions. Including, for example, nitrogen fixation as an additional trait adding PFT diversity may become more relevant in future compared to current conditions and alter our results. Moreover, our model setup does not take into account coexistence through trade-offs in stoichiometric ratios for nutrient requirements (Göthlich and Oschlies, 2012), and thus omits potential effects of stoichiometrically imbalanced nutrient supply on both diversity and productivity (Gross and Cardinale, 2007;Cardinale et al., 2009). Also, other aspects of plankton diversity such as algal mixotrophy (Hartmann et al., 2012;Ward et al., 2011) or spectral light use (Stomp et al., 2007;Hickman et al., 2010) may play different roles in current and future climates. Future investigations covering an extended diversity range and employing a larger variety of traits, tradeoffs and different random phytoplankton assemblages may help to confirm or revise our findings.
Reduced mixing due to enhanced stratification of the upper ocean and thereby lower supply of nutrients into the surface layers is one of the predicted consequences of a warming climate (e.g. Sarmiento et al., 2004;Steinacher et al., 2010). Here we employ an idealized reduction in mixing to investigate the model's ability to capture phytoplankton diversity effects on PP related to environmental changes. The sudden reduction in nutrient supply implied by our idealized model setup is not realistic in terms of expected timescales of real ocean change. A more realistic setup requires a model configured to perform climate change simulations over centennial timescales (e.g. Dutkiewicz et al., 2013), possibly with lower spatial and temporal resolution, or fewer state variables. Our approach also excludes effects of deep ocean equilibration, which may increase the nutrient supply to the upper ocean particularly in upwelling regions on longer timescales. Phytoplankton community structure and diversity adapts to the changed environmental conditions well within the 10 yr timescale. The reductions in PP observed in the 10th year after reducing mixing are of the same order of magnitude as predictions from centuryscale simulations (Steinacher et al., 2010;Taucher and Oschlies, 2011). Thus the sensitivity experiments performed here provide a reasonable framework for investigating model behaviour in terms of potential diversity effects. However, they do not take into account the temperature changes that will accompany the environmental change. Warming may enhance PP (Taucher and Oschlies, 2011;Dutkiewicz et al., 2013), possibly mitigating some of the decrease from a reduction in nutrient supply.

Conclusions
One of the hypothetical diversity effects under environmental change is that less diverse communities show a larger decrease in primary production than more diverse communities, because with more species present it is more likely that the existing niches can be filled also under different environmental conditions. Our simulations comprise diverse communities with up to 78 phytoplankton subtypes grouped in four PFTs and distinguished by temperature, light and nutrient use. However, the less diverse communities occupy a similar trait space as more diverse communities, so that the model essentially only resolves two niches with respect to nutrient use. Consequently, primary production changes are independent of overall diversity within the PFTs, as more diverse communities do not occupy more niches. The magnitude and sign of these changes is affected by type identity, but it is not linked to the number of subtypes. Our results thus show that niches as currently represented in typical ocean biogeochemical models are directly related to the PFTs defined. Which effects of diversity on ocean ecosystem functioning under environmental change are captured by such models depends on the kind and number of trade-offs between different PFTs. Adding diversity in traits other than nutrients and light use as examined here, such as size, may reveal stronger diversity effects. Capturing changes in the niches occupied under changing environmental conditions might thus require formulations allowing niche plasticity or niches not occupied under current conditions. The effect of a mixing reduction on total primary production (PP) can be decomposed into effects of changes in total biomass, biomass of individual subtypes, or community composition (cf. Fig. 3) according to the following equations. Here, PP std and P std denote the total PP and total phytoplankton biomass, respectively, and pp std i is the specific growth rate (in d −1 ) of phytoplankton subtype i in the standard run. Consequently In the model simulations, the individual effects to very good approximation sum up to the total effect, so that PP ≈ PP gro + PP phy ≈ PP gro + PP bm + PP comp .