Nonrandom species loss in phytoplankton communities and its effect on ecosystem functioning

The effects of species diversity on ecosystem functioning have been broadly studied, mostly considering random artificial assemblages. However, natural communities are shaped by ecological interactions and environmental conditions often leading to nonrandom species extinctions. Here, we manipulated a natural phytoplankton community by generating a taxonomic diversity gradient based on rare species exclusions and conducted a mesocosm experiment to investigate the diversity effects on ecosystem functioning (resource use efficiency and biomass) under two nutrient levels. We quantified two functional traits (size and photosynthetic pigments) to evaluate the relation of functional diversity and ecosystem functioning. In a second experimental phase we simulated temperature fluctuations to assess the role of diversity on temporal stability of ecosystem functioning. We did not find a significant effect of diversity on ecosystem functioning and the temporal stability of ecosystem functioning, regardless of nutrient level. These results indicated that loss of biomass caused by rare species extinctions was compensated by the species retained in the diversity gradient. Phytoplankton size diversity was positively related to diversity, but this was not transferred into a positive diversity effect on ecosystem functioning. Additionally, the loss of species did not result in a loss of pigment diversity. The lack of a biodiversity–ecosystem functioning (BEF) relationship in our study may be due to the weak coupling of functional and species diversity and a low manifestation of functional diversity under the evaluated conditions. We emphasize that more realistic biodiversity loss scenarios in experiments can yield different results from those in classical BEF research paradigms.

The relationship between biodiversity and ecosystem functioning (BEF) has been broadly studied in aquatic ecosystems in order to understand the possible effects of species loss on ecosystem processes and the generation of management tools to mitigate these effects Gamfeldt et al. 2015;O'Connor et al. 2017). The positive effect of diversity on ecosystem functioning, as well as its maintenance over time, has been related to the assumption that communities with higher number of species comprise higher variability in function (Cadotte et al. 2011;Weithoff and Beisner 2019). Higher functional diversity or variability promotes complementarity in resource use or niche occupancy enhancing or sustaining ecosystem functions (Ptacnik et al. 2008(Ptacnik et al. , 2010Striebel et al. 2009b;Vallina et al. 2017). Trait diversity has been used to evaluate phytoplankton functional diversity and explain the relationship between species diversity and ecosystem functioning (Litchman et al. 2010;Behl et al. 2011;Hodapp et al. 2016). In this sense, phytoplankton photosynthetic pigments and cell size have been described as import functional traits. While pigment composition is linked to light absorption and can be relevant for carbon fixation, size distribution might be relevant for nutrient uptake strategies influencing primary production (Striebel et al. 2009a;Litchman et al. 2010;Acevedo-Trejos et al. 2018).
However, observations in freshwater ecosystems show that species diversity and trait or functional diversity are not always coupled (Kruk et al. 2017;Cardoso et al. 2017). On one side, natural communities are shaped by environmental filtering processes that can result in communities composed by organisms similar in function independent of the number of species in the community (Schwaderer et al. 2011;Vallina et al. 2017). On the other side, environmental variability might influence species diversity providing heterogeneous conditions where species with different traits can coexist and express (Hillebrand and Matthiessen 2009;Smith et al. 2016).
Under temporal environmental variations, communities with high diversity are expected to buffer temporal changes maintaining ecosystem functioning by presenting higher temporal complementarity and functional redundancy (Cardinale et al. 2007;Ptacnik et al. 2010;Reich et al. 2012). Thus, the positive relation between taxonomic and functional diversity is affected by environmental characteristics and may drive the diversity effect on ecosystem functioning (standing stock and temporal stability) (Cadotte et al. 2011;Weithoff and Beisner 2019).
Modifications of nutrient availability generate strong changes of phytoplankton biomass and richness as well as in the relation of these variables (Dokulil and Teubner 2000;Tian et al. 2017). Phytoplankton species diversity and biomass are positively correlated under low nutrient concentrations, while a negative correlation is expected under high nutrient levels (Vallina et al. 2014;Tian et al. 2017). Thus, high nutrient concentrations may increase phytoplankton biomass but can generate a decrease in diversity favoring the dominance of the most competitive species (Elliott et al. 2006).
Despite the wide evidence available for phytoplankton BEF relationships, most studies use field observational data or experiments with randomly assembled communities composed of species that differ in origin and without a common ecological history (Hillebrand and Matthiessen 2009;Gamfeldt et al. 2015). These artificial assemblages might have an indirect effect on trait diversity, e.g., when the species included in the experiment are explicitly chosen to represent variable taxonomic groups they also increase the diversity of traits in the community. Experiments of this type have been criticized for omitting the nonrandom assembly of species and order of species extinctions common in natural ecosystems (Smith and Knapp 2003;Srivastava and Vellend 2005).
Different approaches have been developed to analyze effects of more realistic nonrandom species loss, considering aspects as sensitivity of species to environmental stress, body size, and rarity of species (Solan 2004;Engel et al. 2017). Multitrophic experiments applying nonrandom species approaches indicate that loss of rare species has a strong negative impact on ecosystem processes (Bracken and Low 2012;Pendleton et al. 2014). But, studies that consider one trophic level suggest that common productive species are more important than rare species in maintaining ecosystem functioning (Smith and Knapp 2003;Solan 2004). In this scenario, the loss of dominant species production is not replaced by less abundant species when they are removed, while rare species loss does not affect ecosystem functioning (Smith and Knapp 2003;Walker and Thompson 2010;Yoshihara et al. 2019). Experiments conducted using ordered species losses have been used in terrestrial plants (Smith and Knapp 2003;Yoshihara et al. 2019), bacteria (Roger et al. 2016), benthic (Solan 2004), and fish communities (Pendleton et al. 2014). However, to our knowledge, BEF relationships based on nonrandom species extinctions have not been analyzed in phytoplankton communities.
In the present study we manipulated a natural phytoplankton community to generate a species diversity gradient by inducing exclusion of rare species, and tested the relationships among this diversity gradient, trait diversity, and ecosystem functioning under different nutrient conditions at a mesocosm scale. Rare species are often more sensitive to environmental changes getting extinct faster because of their small population sizes (Pimm et al. 1988). Thus, the manipulation of a natural community to create a nested rare-species gradient can be used as a more realistic scenario to evaluate the effects of species diversity on trait diversity and ecosystem functioning under controlled conditions (Smith and Knapp 2003). This approach avoids possible site-dependent effects but includes treatments with similar species richness as in natural systems. As diversity effects on ecosystem functioning might be mediated by environmental variability (as mentioned before), we included variability in terms of temperature fluctuations and assessed the ecosystem functioning temporal stability under different levels of diversity and nutrient concentrations. For that, in a second experimental phase (within the same setup) we evaluated the phytoplankton responses to simulated temperature fluctuations. Considering a scenario of nonrandom loss of species in a natural community, we hypothesize that: Hypothesis 1: Communities with higher species diversity comprise greater variability in function and therefore higher functional trait diversity. Consequently, we expect phytoplankton size and pigment diversity to reflect the species diversity gradient.
Hypothesis 2: Rare species loss has a weak effect on ecosystem functioning as the deterministic sequence of species exclusion by rarity implies that the common species performing (for function) in a higher extent prevail in all treatments and might compensate for functions lost by rare species extinction. The effect of species loss on ecosystem functioning is nutrient level dependent.
Hypothesis 3: Species diversity stabilizes ecosystem functioning under changing environmental conditions. Diverse communities are more likely to contain species favored by the different environmental conditions maintaining function over time. Even considering a rare species gradient, we expect a positive effect of diversity on the temporal stability of ecosystem functioning as different species can be favored under environmental fluctuations, thus changing their relative contribution in the community composition. Higher nutrient level may destabilize ecosystem functioning counteracting the diversity effect.

Experimental setup
We conducted an indoor mesocosm experiment including 12 units with 600 liter capacity, build-in rotor to prevent wall Gerhard et al.
Phytoplankton nonrandom species loss growth, temperature sensors, LED lights that simulate daylight and manual mixing discs (Gall et al. 2017). The experiment had a total duration of 43 days subdivided into two phases ( Fig. 1). In the first phase (from day 1 to 27), we tested the effects of phytoplankton diversity and nutrient supply on ecosystem functioning (estimated as standing biomass and resource use efficiency) and the relation between phytoplankton taxonomic and functional trait diversity (estimated as pigment diversity and size diversity). The treatments included a rare-species diversity gradient of six levels (increasing from D1 to D6) and two nutrient conditions (low and high) in which the nitrogen (N) and phosphorous (P) concentrations were manipulated (Fig. 1). In the second phase (from day 28 to 43) we evaluated the temporal stability of ecosystem functioning estimators across the diversity gradient and nutrient level. For that, we included environmental variability generating four serial temperature fluctuation events in all treatments ( Fig. 1).
To develop the phytoplankton diversity treatments, a natural community from a German lake (53 33 0 05 00 N; 7 58 0 49 00 E) was collected at the end of the summer (2017). The lake water was filtered through a 53 μm mesh to remove the zooplankton grazers. The diversity treatments were generated by applying a dilution method to simulate nonrandom species loss resulting in a rare-species gradient (Engel et al. 2017;Hammerstein et al. 2017 The indoor mesocosms were filled with distilled water and the components of the phytoplankton WC growth medium (except for N and P) were added to avoid the limitation by these elements (Guillard and Lorenzen 1972). N and P were added as NaNO 3 and K 2 HPO 4 to generate two nutrient concentration levels but maintaining the N : P ratios in a narrow range. The realized total nutrient concentrations during the experiment showed a mean of: 20.7 μmol L −1 (SD = 3.6) and 45.0 μmol L −1 (SD = 4.9) for N, and 2.3 μmol L −1 (SD = 0.7) and 3.7 μmol L −1 (SD = 0.8) for P, in the low and high nutrient levels.
The phytoplankton inoculation and nutrient addition were done at the beginning of the experiment (8 October 2017). Each phytoplankton diversity level was used to inoculate two mesocosms, one corresponding to the low nutrient conditions and the other to the high nutrient conditions. The same amount of algal biomass was added to each mesocosm through the correction of the inoculum volume based on optical density data (inoculum volumes between 80 and 150 mL). The mesocosms were mixed manually every day in the morning (before the lights turned on) to avoid sinking losses and ensure homogeneous conditions. The day-night cycle was of 12 : 12 h at a light intensity of about 300 μmol photon m −2 s −1 . The temperature in the units was constant at 20 C in the first experimental phase (temperature similar to the lake temperature in summer). The temperature fluctuations were simulated between sampling dates during the second phase of the experiment. Temperature was increased gradually by 3 C within 24 h and afterwards decreased to 20 C again. This was repeated four times (see Fig. 1) during the second experimental phase. Thus, deterministic temperature fluctuations were simulated to compare ecosystem functioning responses under temperature variations within a natural range for lakes (phase II) to constant experimental conditions (phase I).

Samplings and laboratory analysis
Optical density and in vivo chlorophyll a (Chl a) were measured daily to monitor the phytoplankton biomass. Optical density (wavelength 440-450 nm) was measured using a custom-tailored device and Chl a was measured with hand fluorometer (AquaFluorTM; Turner Designs). The first sampling was done 1 week after the phytoplankton inoculation (day 7). Afterwards, samplings were performed every 4 days until the end of the experiment (day 43). On every sampling day, 10% (60 liter) of water was replaced with water containing the design nutrient concentrations. Samples were taken to analyze nutrients and phytoplankton communities.
Samples for pigment composition, particulate organic carbon (POC) and particulate organic nitrogen were filtered onto acid-washed pre-combusted glass-fiber filters (Whatman GF/C). Pigments were extracted with 10 mL ethanol (90%) and the absorbance spectrum was measured photometrically (AquaMate Plus UV-Vis, Thermo Scientific) (400-700 nm) for theidentification and quantification according to Thrane et al. (2015). Filters for particulate organic carbon and nitrogen were dried at 60 C and measured using an elemental analyzer (Flash EA 1112, Thermo Scientific). Samples for the dissolved nitrogen fractions (NO x − and NH 4 + ) were filtered (0.45 μm surfactant-free cellulose acetate syringe filters) directly into Epicaps. The samples were poisoned with HgCl 2 and analyzed photometrically, using a multiscan GO microplate spectrophotometer (Thermo Science). NO x − (NO 3 − and NO 2 − ) was determined following the method described by Schnetger and Lehners (2014) and ammonium (NH 4 + ) was analyzed after a modified version of Benesch and Mangelsdorf (1972) method. Total phosphorus was measured to control P concentration in the experimental units over time by molybdate reaction of unfiltered water samples after digestion with potassium peroxydisulfate (K 2 S 2 O 8 ) solution (Wetzel and Likens 2013).
The size distribution of the phytoplankton community was determined using a cell counter (Beckman Coulter, Z 2 ) (except for the first sampling because the equipment was not available). Cell sizes (as equivalent spherical diameter, ESD) were distributed in 255 classes and the cell abundance of each size class was obtained. The equivalent spherical diameter is obtained through electrodes that sense the volume of electrolyte displaced by each cell.
To analyze the phytoplankton communities 30 mL sample were fixed with 1% lugol solution. Samples corresponding to the inoculated communities (day 0), start sampling (day 7), the end of the first experimental phase (day 27) and the end of the second experimental phase (day 43) were counted and identified under inverted microscope based on Utermöhl (1958). Phytoplankton was identified to the species level and morphospecies were used when clear assignment of a species name was not possible. Mesozooplankton absence was verified by counting 20 liter water samples that were filtered through 105 μm mesh in every sampling day (taken from the exchanged water).

Statistical analysis
All statistical analyses and figures were performed in R version 3.5.0 (R Development Core Team 2018) and a level of α = 0.05 was considered for statistical significance in the analyses.

Treatment effects on community species diversity
Data of phytoplankton species composition corresponding to the first experimental phase were used to evaluate how taxonomic diversity reflected the dilution gradient over time in the nutrient levels. For that, phytoplankton species diversity was estimated using species richness and the inverse Simpson diversity index based on phytoplankton species abundance data (ID spp ) with vegan package. We selected the ID spp because of its robustness for diversity estimation integrating richness and evenness (Chase and Knight 2013), but also richness because the rare-species loss used to develop the diversity gradient is reflected by the number of species. Because we only had community composition information for two samplings (day 7 and day 27), repeated-measures ANOVAs were used to analyze the phytoplankton species diversity responses to the treatments (diversity and nutrients) and its changes over time (n = 24). The diversity gradient was included as continuous variable assigning numbers from one to six for each level (increasing from D1 to D6), and nutrient level was used as factor with two categories (low and high). Species richness and ID spp were ln-transformed to reduce over dispersion and improve the fit to a Gaussian distribution of the residuals.

Treatment effects on community functional diversity (H1)
Data corresponding to the first experimental phase were used to test H1. Pigment composition and size diversity of the phytoplankton communities were used to evaluate functional diversity. These trait diversity estimations were assessed as community weighted variability to analyze the functional diversity in the community, not as absolute values. We calculated the inverse Simpson index to estimate pigment diversity based on the extracted pigment concentrations (ID pig ). We used the cell size and abundance obtained with the cell counter to calculate the weighted mean cell size (S w ) and its weighted variance (V w ) as estimation of where A is the abundance of each size class i, S is the size of each size class i, and n is the total number of size classes (255 in our case). The responses of pigment diversity, mean size and size diversity to the treatments were analyzed using linear mixed models (following the protocol recommended by Zuur et al. (2009)). We included two and three-way interactions considering the diversity gradient, nutrient level and time as fixed variables. For the random component we tested models including random intercept (of each mesocosm), random slope (of each mesocosm over time), and both random intercept and slope together. Because variables were measured over time in the same units (repeated measures design) we screened for temporal autocorrelation through the autocorrelation function (acf) and the analysis of residuals. Due to all variables showed significant autocorrelation we included the residual AR-1 autocorrelation structure in the models. Bestfit model selection was performed using Akaike Information Criterion (AIC) and analysis of residuals. The explained variance was estimated as marginal and conditional R 2 (R 2 m and R 2 c, respectively) using r.squaredGLMM function in the MuMIn package. Because R 2 values were low, we ran the models excluding diversity treatment D1 (D1 showed particular patterns in the community traits measurements). Explained variance increased in the new models, but treatments effects were similar in the models including or excluding D1 (see results section). Data of the six samplings, corresponding to the first experimental phase, were used for the pigment analysis (n = 72), while data for size of the first sampling (day 7) were missing (n = 60). All models were performed using lme function (nlme package).

Treatment effects on ecosystem functioning (H2)
Data from the first experimental phase were used to analyze the ecosystem functioning responses to the treatments and test H2. Variables considered as indicators for ecosystem functioning were the standing biomass (using POC, Chl a and optical density as estimations) and resource use efficiency (RUE). The phytoplankton RUE was calculated as the ratio between standing biomass (POC) and the total concentration of the limiting nutrient (Ptacnik et al. 2008). Because N was the limiting nutrient in all experimental treatments (tested by a nutrient limitation bioassay, see results in Supplementary Information Appendix 1 Fig. A5), we used total nitrogen (particulate organic nitrogen plus dissolved nitrogen as NO x − and NH 4 + ) to calculate RUE.
Different estimations of phytoplankton standing biomass were compared with Spearman correlations using the whole dataset (n = 120). POC was positively correlated to optical density (r = 0.90; p < 0.001). Chl a showed also a significant positive correlation with POC (r = 0.78; p < 0.001) and optical density (r = 0.79; p < 0.001), but with lower coefficients. To analyze the treatment effects over time on the standing biomass and the RUE, we followed the same procedure as used for the functional diversity data (n = 72 for each variable). Final models included AR-1 autocorrelation structure and R 2 estimations were high for the selected models. Spearman correlations were used to identify possible relations between phytoplankton biomass estimations and trait diversity, independent of the taxonomic diversity during the first experimental phase (n = 72).

Treatment effects on ecosystem temporal stability (H3)
We tested H3 by analyzing the treatment effects on the temporal variability of standing biomass and RUE. For this, we used the dataset of both experimental phases (phase I n = 72 and phase II n = 48) which represented temporal constant conditions (constant temperature) and variable conditions when temperature fluctuations were applied. We extracted the residuals of the linear mixed models used for H2 including the complete data set and tested their variability between experimental phases. The variance of the residuals was calculated for each mesocosm in each experimental phase and GLMs were performed to analyze changes in variability according to the nutrient level, the diversity gradient and the experimental phase (n = 24). To identify possible changes in the phytoplankton community composition in response to the temperature fluctuations we conducted an analysis of similarity (ANOSIM) using Bray-Curtis index with 999 permutations between experimental phases (day 27 vs day 43). Relative species abundances were used for the analysis to avoid the effect of changes in biomass over time.

Diversity gradient generated by dilution (inoculum)
The phytoplankton communities used to inoculate the mesocosms (at day 0) contained a total of 89 species and were dominated by chlorophytes in all cases (between 60% and 99% of the relative abundance) (Supplementary Information Appendix 1 Fig. A4). The richness corresponding to the inoculated diversity levels ranged from 7 to 52 and the Inverse Simpson diversity (ID spp ) ranged from 1.1 to 12.3 from the lowest to the highest diversity treatment (D1 to D6) (Supplementary Information Appendix 1 Figs. A2, A3).

Treatment effects on species diversity
During the first experimental phase, the phytoplankton species richness in the mesocosms reflected the general pattern of the diversity gradient (significant diversity effect, F 1,8 = 28.7; p < 0.001, Fig. 2); but decreased over time showing lower richness at the end of phase I (day 27) than at the beginning of the experiment (day 7) (significant time effect F 1,8 = 6.5; p = 0.03, Fig. 2). The ID spp also showed a tendency to increase with the diversity gradient treatment, but statistical significance was marginal (F 1,8 = 4.5; p = 0.066, Fig. 2), and changes over time were not found (Fig. 2). Thus, the initial diversity gradient generated by rare species loss was principally reflected by richness over time. Changes in the measured diversity during the experiment were expected as a response to the treatments. The initial diversity manipulation affected not only the number of species but also the community composition via inter-specific interactions (e.g., competition) and potential responses to nutrient availability. No significant effects of the nutrient level were found on the species diversity indices.

Treatment effects on community functional diversity (H1)
While mean size and size diversity showed significant treatment effects, pigment diversity was not affected by diversity and nutrients (Table 1; Fig. 3). Because trait diversity showed particular patterns in treatment D1, models including and excluding D1 were compared. While only slight differences were detected in the model outputs for size diversity with and without D1, the model explained variance (R 2 c) increased when D1 was excluded. Size diversity showed a positive effect of the diversity gradient. The diversity effect corresponded to a main effect in the model excluding D1, but a diversity-time interactive effect in the model including D1 (Table 1; Fig. 3b). This is because the lowest diversity level (D1) presented high size diversity which increased over time, deviating from the general pattern (Fig. 3). The high nutrient level showed higher size diversity than the low nutrient level, and the difference increased over time (Table 1; Fig. 3). Like size diversity, the phytoplankton mean size showed a positive relation with the taxonomic diversity gradient (Table 1; Fig. 3c). Thus, higher abundance of larger cells was present across the species diversity gradient increasing the mean size, but also increasing size diversity as smaller cells were maintained. Phytoplankton pigment diversity, size diversity and mean size increased over time in the treatments (Table 1; Fig. 3).

Treatment effects on ecosystem functioning (H2)
None of the variables used to estimate phytoplankton biomass (POC, Chl a, optical density) nor RUE showed a significance response to diversity as main effect or diversity interactive effects. A significant positive nutrient effect was detected for Chl a over time, while RUE increased with low nutrients (Table 1; Fig. 4; Supplementary Information Appendix 1 Fig. A6). All biomass estimators and RUE increased over time (Table 1) as expected considering that a small volume of phytoplankton was inoculated in the mesocosms. Although the interactive effects of diversity and nutrients were not significant (Table 1), ecosystem functioning estimators showed different patterns between nutrient treatments. These patterns were similar for RUE, POC, and Chl a (Fig. 4, Supplementary Information Appendix 1  Fig. A6). Under low nutrient concentrations, the diversity treatment D2 presented the highest RUE and biomass. Under high nutrient concentrations the lowest RUE and biomass occurred in D1 and increased or was maintained across the diversity levels with D6 presenting the highest values at the end of the experiment ( Fig. 4; Supplementary Information Appendix 1 Fig. A6). Contrary to POC and Chl a, the biomass measured as optical density was highest for D2 in both nutrient levels ( Supplementary Information Appendix 1 Fig. A6). Correlations between trait functional diversity and biomass estimators only showed a significant positive correlation between POC and  pigment diversity (r = 0.41; p < 0.001) (size diversity and POC, e.g., were not significantly correlated with a correlation coefficient and p value of r = −0.1; p = 0.44).
Treatment effects on ecosystem temporal stability (H3) Residual variance of biomass and RUE did not show differences between experimental phases (Table 2), suggesting that the variability in standing biomass and RUE did not change between experimental constant and fluctuating conditions. RUE and Chl a residual variance responded to the nutrient level. RUE variability increased under low nutrients (significant nutrient effect) while Chl a variability was higher under temperature fluctuations and high nutrients (significant nutrient × phase effect, Table 2).
The phytoplankton community composition was dominated by chlorophytes in the different treatments at the different samplings (between 83% and 99% during the experiment). The ANOSIM did not show significant differences of species composition between experimental phases using relative abundance data (R = −0.05; p = 0.9) (Supplementary Information Appendix 1 Fig. A7). However, the species abundance between the experimental phases showed a decrease in dominant species abundance after temperature fluctuations (points below the 1 : 1 slope) (Fig. 5).

Rarity diversity gradient
In this study, we showed that the dilution-method used prior to the setting up of the experiment was successful for generating a phytoplankton taxonomic diversity gradient from a natural community simulating nonrandom species loss (Engel et al. 2017;Hammerstein et al. 2017). However, the taxonomic diversity gradient obtained did not have a positive effect on ecosystem functioning and ecosystem functioning temporal stability (supporting H2 but contrary to H3). Considering that we used a nonrandom species loss approach based on rarity, our results support the idea that the loss of rare species have a weak effect on total producers' biomass as these losses are compensated by common species (Smith and  Table 2. Results of GLM models for residual variance of the variables corresponding to the first and second experimental phases. F 1,16 and p values (between brackets) are presented. Significant effects (p<0.05) are highlighted in bold. Div, diversity treatment; Nut, nutrient level; Phase, experimental phase; POC, particulate organic carbon; Chl a, chlorophyll a; OD, optical density; RUE, resource use efficiency.

POC
Chl a OD RUE  Yoshihara et al. (2019) showed that random species losses in terrestrial plants have stronger effects on yield than nonrandom losses, suggesting that experiments based on artificial random assemblages might be conditioned to find significant diversity effects on ecosystem functioning (Srivastava and Vellend 2005;Gamfeldt et al. 2015). In our experiment, the lowest diversity treatment (D1) showed lower RUE and standing biomass compared to the other diversity levels. This might be explained by the fact that most species were removed in this treatment, not only the rare species but also common productive species. Solan (2004) proposed that ecosystem functioning is affected by the order of species extinction when this sequence involves the loss of traits related to the analyzed function, which might not be the case in our experimental setup where species loss was driven by rarity. We did not find significant interactive effects of diversity and nutrients on biomass as expected in H2, but the diversity treatments responded differently under low and high nutrient levels. The second diversity level (D2) represented almost a monoculture of Monoraphidium contortum and showed an unexpected high standing biomass, especially in the low nutrient treatment (see Fig. 4 and Supplementary Information Appendix 1 Fig. A6). The dominance by a single species can result in high biomass when individual species present particular traits that favor high efficiency of the limiting resource use (species identity effect) (Cardinale 2011;Lewandowska et al. 2016). M. contortum has been described as a highly productive species under N limitation in monocultures (Bogen et al. 2013), but also might dominate natural communities when N is limiting (Ferragut and de Campos Bicudo 2012). Since our treatments were N limited and we found that D2 showed the highest RUE and biomass under low nutrients, our results support previous findings suggesting that M. contortum was successful in the use of nitrogen, especially under low nutrient concentrations. However, high nutrient concentrations showed a tendency to increase biomass in treatments with higher diversity over time.

Treatment effects on ecosystem functioning mediated by traits
We obtained different relationships between the taxonomic diversity gradient and the measured traits: while size diversity showed a positive relation with taxonomic diversity, pigment diversity was not coupled to taxonomic diversity. This means that contrary to our expectations in H1, the generated taxonomic diversity gradient did not necessarily result in a gradient of trait diversity. The potential positive effect of diversity on ecosystem functioning implies that with higher species diversity a higher diversity in traits and consequently in function is present (Cadotte et al. 2011;Weithoff and Beisner 2019). Thus, the absence of a positive BEF relationship in our study might have been influenced by the lack of a relationship between the diversity of pigments (trait related to biomass production) and the species diversity. This was supported by the positive correlation between pigment diversity and particulate organic carbon suggesting that pigment diversity influenced carbon assimilation, and might have restricted the detection of positive BEF relationships as pigment diversity (i.e., functional diversity) was not coupled to taxonomic diversity (Cadotte et al. 2011;Fontana et al. 2018).
Contrary to our results, previous experimental studies found a positive relationship between species diversity and pigment diversity (Striebel et al. 2009a). However, phytoplankton species and functional diversity are not always coupled (Kruk et al. 2017;Cardoso et al. 2017;Fontana et al. 2018). In natural ecosystems, environmental filtering processes can restrict the functional diversity present in the communities (Schwaderer et al. 2011;Vallina et al. 2017). Thus, different results in the relation between species and functional diversity obtained in field and experimental studies might be explained by the presence of higher phylogenetically related species and similar functionality in natural compared to artificial communities (Behl et al. 2011). In our experiment, similar pigment diversity across the species diversity gradient could be influenced by the dominance of chlorophytes in all the experimental units (Supplementary Information Appendix 1 Fig. A7). Because pigment composition is phylogenetically constrained (Rowan 1989;Lenning et al. 2004), the relation between species and pigment diversity might not be evident when most of the species in the community belong to the same taxonomic class. Supporting this, a general analysis of freshwater phytoplankton indicated that traits tend to be conserved implying that phylogenetically related species are more likely to present similar functional traits (Bruggeman 2011).
In our experiment, size diversity was positively related to taxonomic diversity and showed a positive nutrient concentration effect. As phytoplankton size is related to nutrient uptake, organisms that differed in size can respond differently to resource availability favoring complementarity in resource use and biomass production (Litchman et al. 2010;Marañón 2015;Acevedo-Trejos et al. 2018). However, we did not find significant responses of RUE and biomass to diversity that could be related to the size diversity increase, and only Chl a increased with high nutrients. These results suggested that even if the size diversity and mean size were positively related to taxonomic diversity, no potential positive size-mediated effect on ecosystem functioning was detected. Accordingly, we did not find a significant correlation between size diversity and biomass. Phytoplankton models propose that diversity in traits associated to nutrient uptake strategies have a weak effect on production (Vallina et al. 2017;Chen et al. 2019); but also that positive phytoplankton size diversity effects on production can be mediated by environmental variability (Smith et al. 2016;Chen et al. 2019). Thus, the low nutrient supply variability in our experiment (where treatments had different total concentrations but nitrogen was the limiting nutrient in all treatments over time) might have influenced the absence of association between phytoplankton size diversity and biomass.

Treatment effects on ecosystem functioning temporal stability
Contrary to H3, environmental variability generated by temperature fluctuations did not increase the variability in ecosystem functioning in comparison with the constant temperature phase, and no general changes in relation to diversity were detected in our experiment. As the richness gradient in the treatments was maintained until the end of the experiment, the lack of diversity effects cannot be explained by a change in these levels (see Supplementary Information Appendix 1 Figs. A8, A9). RUE variability increased under low nutrients and the responses of Chl a variance to the nutrient level change between the experimental phases. Thus, a stabilizing effect of diversity and destabilizing effect of higher nutrients on ecosystem functioning as shown by Ptacnik et al. (2008) in lakes were not identified as a general pattern in our experiment. These results were supported by no significant differences in species composition between the experimental phases (constant vs. fluctuating conditions), suggesting that the communities maintained the variability in function under fluctuating conditions without changes in the community composition. We found a decrease in abundance of highly dominant species after the temperature fluctuations, but it was not enough to significantly change the composition of the communities and the biomass responses in the timeframe applied. As we used a natural community for the experimental setup, species could be adapted to natural temperature variations so that the fluctuation events did not present a strong change for the species present in the community.

Experimental considerations
We used an experimental design that minimized common laboratory experiment limitations based on artificial assemblages (Srivastava and Vellend 2005). However, other limitations were present in our approach for BEF relationships evaluation. The experiment was conducted in controlled closed systems where species dispersal was not possible as occurring in natural ecosystems. This might have caused the decrease in richness detected over time in the experimental units as possible migration effects on diversity persistence were not included (Hodapp et al. 2016). Furthermore, even if we tried to incorporate environmental variability simulating temperature fluctuations, the experimental systems had low environmental variability in comparison with natural systems. As environmental heterogeneity has been suggested to mediate positive diversity effects on ecosystem functioning (Ptacnik et al. 2010;Hodapp et al. 2016;Smith et al. 2016), the lack of high environmental variability (e.g., number of environmental variables considered, magnitude of variation) could have limited our findings.

Conclusions
In this study we use for the first time a nonrandom species loss approach to evaluate species diversity effects in combination with nutrient supply on ecosystem functioning of a natural phytoplankton community. We found that the diversity gradient generated by a rare species loss did not show an effect on ecosystem functioning (RUE and standing biomass) or temporal stability of ecosystem functioning under different nutrient concentrations. We proposed that patterns showed by the trait functional diversity in addition to identity effects, might have influenced the lack of a positive BEF relationship. Beyond experimental limitations, the novel experimental design presented in our work combined with the poorly understanding of mechanism explaining BEF relationships in phytoplankton communities highlight the importance of the patterns presented in our study. The generation of new approaches to assess the relative importance of different aspects of diversity on ecosystem functioning considering more realistic scenarios is crucial to evaluate ecological responses to perturbations and environmental fluctuations. However, progress is needed in applying these new approaches to phytoplankton communities.