Phytoplankton diversity affects ecosystem functioning in a coastal upwelling system

ORIGINAL RESEARCH published: 26 November 2020 doi: 10.3389/fmars.2020.592255 Frontiers in Marine Science | www.frontiersin.org 1November 2020 | Volume 7 | Article 592255 Edited by: Xosé Anxelu G. Morán, King Abdullah University of Science and Technology, Saudi Arabia Reviewed by: Francisca C. García, University of Exeter, United Kingdom Savvas Genitsaris, International Hellenic University, Greece Wang Tian, North China Electric Power University, China *Correspondence: Jaime Otero [email protected] †Present address: Jaime Otero, Centro Oceanográfico de Vigo, Instituto Español de Oceanografía, Vigo, Spain Specialty section: This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Marine Science Received: 06 August 2020 Accepted: 20 October 2020 Published: 26 November 2020 Citation: Otero J, Álvarez-Salgado XA and Bode A (2020) Phytoplankton Diversity Effect on Ecosystem Functioning in a Coastal Upwelling System. Front. Mar. Sci. 7:592255. doi: 10.3389/fmars.2020.592255 Phytoplankton Diversity Effect on Ecosystem Functioning in a Coastal Upwelling System Jaime Otero1*†, Xosé Antón Álvarez-Salgado1and Antonio Bode2 1Instituto de Investigaciones Marinas (IIM-CSIC), Vigo, Spain, 2Centro Oceanográfico de A Coruña, Instituto Español de Oceanografía, A Coruña, Spain Species composition plays a key role in ecosystem functioning. Theoretical, experimental and field studies show positive effects of biodiversity on ecosystem processes. However, this link can differ between taxonomic and functional diversity components and also across trophic levels. These relationships have been hardly studied in planktonic communities of coastal upwelling systems. Using a 28-year time series of phytoplankton and zooplankton assemblages, we examined the effects of phytoplankton diversity on resource use efficiency (RUE, ratio of biomass to limiting resource) at the two trophic levels in the Galician upwelling system (NW Iberian peninsula). By fitting generalized least square models, we show that phytoplankton diversity was the best predictor for RUE across planktonic trophic levels. This link varied depending on the biodiversity component considered: while the effect of phytoplankton richness on RUE was positive for phytoplankton RUE and negative for zooplankton RUE, phytoplankton evenness effect was negative for phytoplankton RUE and positive for zooplankton RUE. Overall, taxonomic diversity had higher explanatory power than functional diversity, and variability in phytoplankton and zooplankton RUE decreased with increasing phytoplankton taxonomic diversity. Phytoplankton used resources more efficiently in warmer waters and at greater upwelling intensity, although these effects were not as strong as those for biodiversity. These results suggest that phytoplankton species numbers in highly dynamic upwelling systems are important for maintaining the planktonic biomass production leading us to hypothesize the relevance of complementarity effects. However, we further postulate that a selection effect may operate also because assemblages with low evenness were dominated by diatoms with specific functional traits increasing their ability to exploit resources more efficiently. Keywords: coastal upwelling, stability, functional diversity, taxonomic diversity, zooplankton, phytoplankton, nutrients, resource use efficiency INTRODUCTION Marine phytoplankton is responsible for roughly half of the global primary production (50 Pg C year−1;Chavez et al., 2011), contributes to nutrient cycling and regulation of climate dynamics, affects the fate of adjacent trophic levels (Richardson and Schoeman, 2004), and, ultimately, constrains fishery catches (Chassot et al., 2010). Marine phytoplankton is an extremely diverse group of organisms (De Vargas et al., 2015), and this diversity, encompassing a large variety of Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling life histories, is an essential factor that affects the whole structure of marine ecosystems (Naeem, 2012). It is indeed the diversity of functional traits that differs among and within species and taxonomic groups, the key component determining the fitness of planktonic communities along environmental gradients and influencing the functioning of pelagic ecosystems (Irwin and Finkel, 2018). Therefore, there is a need to better understand the relationship between the variability in phytoplankton diversity and its effects on ecosystem processes. In recent years, various reviews have synthesized the effects of alterations of biodiversity (B) on ecosystem functioning (EF) concluding that biodiversity has a major role in sustaining the productivity of ecosystems and their stability (Cardinale et al., 2012; Tilman et al., 2014). Most of this evidence comes from controlled experiments; however, field studies are also consistent with theory and experiments demonstrating the strong effect of biodiversity on ecosystem production even after accounting for abiotic forcing (Duffy et al., 2017). In the marine realm, several studies have been devoted also to understand the effects that marine biodiversity has on production, biomass, or on the resilience to disturbances or invasions (Stachowicz et al., 2007; van der Plas, 2019). Yet the majority of this experimental and field research has examined relationships in the benthos (O’Connor and Byrnes, 2014; Duffy et al., 2017). Regarding the pelagos, recent studies have addressed BEF relationships within natural assemblages of plankton either in fresh or marine waters. A BEF relationship in the plankton was first described by Ptacnik et al. (2008), who showed that resource-use efficiency (RUE), an ecological index that measures the proportion of supplied resources turned into new biomass (Hodapp et al., 2019), scaled positively with phytoplankton taxonomic richness in multiple Fennoscandian lakes and in the Baltic Sea. Korhonen et al. (2011) also used data from boreal lakes to connect, in this case, productivity to diversity of various plankton groups showing linear, unimodal, or nonsignificant relationships between richness and biomass production depending on the spatial scale. Furthermore, Olli et al. (2014), using long-term phytoplankton sampling, concluded that increased diversity enhanced RUE for primary producers across the brackish Baltic Sea. However, it is unclear if these BEF relationships apply to highly dynamic Eastern Boundary Upwelling Ecosystems (EBUEs). For instance, in the California Current System, an ecosystem model coupled to a circulation model showed a humpshaped relationship between diversity and productivity with portions of the diversity-productivity scatter being dependent on geographic regions (Goebel et al., 2013). Other authors did not find a relationship between phytoplankton species richness and ecosystem productivity (e.g., Cermeño et al., 2013). Typically, the majority of BEF studies have used richness as the metric of biodiversity because it is easy to manipulate in experiments and to measure in the field. However, there is also evidence of the relevance of other components of biodiversity to understand pelagic processes. For instance, some authors have studied the effect(s) of evenness on planktonic ecosystem properties showing a strong negative effect on RUE in the phytoplankton of the Wadden Sea (Hodapp et al., 2015) or on biomass and resource use in the phytoplankton of the Baltic Sea (Lehtinen et al., 2017) highlighting the importance of the identity of dominant species. Besides taxonomic diversity, biodiversity can be assessed in terms of functional diversity, that is, accounting for the expression of multiple functional traits in the community, often concluding that functional diversity can be a better predictor of ecosystem properties than taxonomic diversity as shown for phytoplankton communities in Fennoscandian lakes (Abonyi et al., 2018). Whether taxonomic or functional diversity of competing species affects ecosystem properties, it is also fundamental to incorporate trophic complexity in order to understand the effects of biodiversity across trophic levels (Duffy et al., 2007). In experimental planktonic systems, Striebel et al. (2012) showed that phytoplankton diversity increased zooplankton productivity, while Filstrup et al. (2014) showed that the effect of phytoplankton evenness on RUE switched from negative at the producer level (phytoplankton) to positive at the consumer level (zooplankton) in US lakes. Apart from average effects on ecosystem functioning, theory, experiments, and field studies predict that increasing diversity can reduce the variability of community biomass or other ecosystem properties in time and space through several mechanisms (Loreau and de Mazancourt, 2013). This would occur because more diverse assemblages containing interacting species, which respond differently to the environment are more likely to buffer the effects of perturbations conferring stability to the community and maintaining the ecosystem properties in a dynamic environment (Ives and Carpenter, 2007). Stability, however, is a complex and multifaceted concept with multiple components, such as variability, resistance, or resilience, which might be unrelated implying that the overall stability of an ecosystem might not be simply explained by one particular component (Hillebrand et al., 2018), thus leading to different biodiversity-stability relationships (Craven et al., 2018). Most analyses on biodiversity-stability relationships have been performed through experiments in terrestrial systems (Tilman et al., 2014), whereas fewer studies dealt with natural ecosystems, and the majority have focused on plants (e.g., García-Palacios et al., 2018; but see Cusson et al., 2015). In the case of natural plankton assemblages, Ptacnik et al. (2008) showed that higher levels of phytoplankton taxonomic richness implied less variability of both resource use and community composition, and Shurin et al. (2007) documented a positive relationship between zooplankton diversity and community stability in temperate lakes. Despite all these research efforts, the importance of phytoplankton diversity on stability in EBUEs, so as the shape of BEF relationships between trophic levels, and the performance of functional diversity vs. taxonomic diversity have been yet poorly addressed (but see Cermeño et al., 2013; Goebel et al., 2013; Vallina et al., 2017). Planktonic communities are highly dynamic with assemblages changing rapidly in response to circulation and fertilization patterns and other physical and environmental forcing. This is even more evident in coastal upwelling systems where planktonic assemblages might fluctuate at short-time scales with different assemblages characterizing the various phases of an upwelling cycle (Marañón, 2015). At larger spatial and temporal scales, the structure of phytoplankton and zooplankton assemblages is Frontiers in Marine Science | www.frontiersin.org 2November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling affected by seasonal changes, with species’ abundance responding to light conditions, temperature, nutrient inputs, or the presence of particular producers and consumers (Wiltshire et al., 2015). Therefore, planktonic biodiversity is affected by the physical environment, water hydrography, and biotic variables (Sarker et al., 2018). However, BEF studies have focused less on the abiotic and biotic context that might exert comparable effects to changes in species richness in mediating ecosystem properties (Godbold, 2012; van der Plas, 2019). Thus, accounting for the environmental context in BEF relationships under natural (or experimental) conditions is crucial to understand and interpret the effects that changing biodiversity has on ecosystem properties (García et al., 2018), stability (García-Palacios et al., 2018), or community performance (Schabhüttl et al., 2013). In this study, we analyzed whether the phytoplankton diversity effects on ecosystem function follows theoretical and experimental expectations in real communities occurring in EBUEs for which knowledge is limited. In doing so, we used longterm time series of phytoplankton and zooplankton community data in conjunction with meteorological and hydrographic data to examine BEF relationships in a highly dynamic coastal upwelling ecosystem. In particular, we (i) evaluated the effect of two components of phytoplankton taxonomic diversity (richness and evenness) on phytoplankton and zooplankton rates of productivity to the amount of available resources (i.e., RUE), (ii) tested whether biodiversity influences planktonic RUE variability, (iii) quantified the importance of biodiversity relative to environmental conditions in driving planktonic RUE dynamics, and (iv) evaluated the explanatory power of phytoplankton taxonomic diversity vs. functional diversity in explaining planktonic RUE. MATERIALS AND METHODS Study Area and Plankton Sampling Galicia is at the northern boundary of the Iberia/Canary current EBUE (Figure 1). Coastal winds at these latitudes (42◦to 44◦N) are seasonal; northerly winds prevail from March– April to September–October, promoting coastal upwelling, and downwelling-favorable southerly winds predominate the rest of the year. However, more than 70% of the variability in coastal winds occurs in periods of <1 month, so that the upwelling season appears as a succession of wind-stress events separated by wind-calm episodes, with a wide variety of frequencies ranging from 3 to 15 days (Álvarez-Salgado et al., 2002). Phytoplankton identification and count data were obtained from the time series project RADIALES conducted by the Instituto Español de Oceanografía off A Coruña (NW Spain, Figure 1) (Bode et al., 2009). Specifically, water samples were collected monthly with 5 L of Niskin bottles or a rosette sampler from 0-, 5-, 10-, 20-, 30-, 40-, and 70-m depths at station E2CO (water depth, 80 m; 43◦25′30′′N, 08◦26′20′′W; Figure 1) from January 1989 to December 2016 (n=313 days sampled), though sampling at 7 m ended in 1993, and from 2010 onward, samples were taken just from 3 depths between the surface and 40 m including the depth of the surface chlorophyll maximum. For each depth, samples were collected for the determination of phytoplankton abundance and inorganic nutrients and chlorophyll aconcentration following the methods described in Casas et al. (1997). Phytoplankton samples of volume 50–100 ml were preserved in Lugol’s solution and kept in the dark until analysis. Depending on phytoplankton concentration, 10–25 ml of samples was allowed to settle in the Utermöhl chamber for up to 24 h. Samples were counted following the technique described by Utermöhl (Lund et al., 1958) using a Nikon Diaphot TMD microscope until May 1997 and a Nikon Eclipse TE300 microscope until the end of the time series. A magnification of 100×was used for large forms, 250×for intermediate forms, and 400×for small forms. The entire slide was examined at 100×to account for large species, while only transects or smaller areas were examined at higher magnification. At least 250 cells were counted for each sample. Whenever possible, organisms were classified at the species or genus level. Species nomenclature followed the World Register of Marine Species (http://www.marinespecies.org). The samples were identified by two experts (M. Varela until 2010, and J. Lorenzo from 2010 to 2016). The individual biomass (in pg C) for each identified taxa was estimated from cell biovolume (in µm3) after measuring the dimensions of 30–100 cells in samples distributed over all seasons and applying conversion equations from the literature (see details in Huete-Ortega et al., 2010). Biomass (in pg C L−1) of a given species in a given day and sampled depth was then calculated as the abundance (in cells L−1) times the cell biomass. For the purposes of the present work, we used the taxa that were systematically identified at the species level in at least 10 samples along the time series. This resulted in a set of 73 taxa, 18 of which started to be identified in 2008 (Supplementary Table 1). The total biomass (obtained as the sum of the biomass of all counted species in each day and sampled depth) was significantly correlated with chlorophyll aconcentration (Supplementary Figure 1), and the biomass was dominated by diatoms (Supplementary Figure 2). The original dataset is available at PANGAEA (https://doi.org/10. 1594/PANGAEA.908815) (Bode, 2019). Additionally, we collated a series of morphological (cell size and ability to form chains or colonies), physiological (silica requirement, trophic strategy, and pigment composition), and behavioral (ability to swim) traits for each phytoplankton species (Supplementary Table 1) using our own compilation and data from the literature (Klais et al., 2017). This group of traits is of relevance for reproduction, resource acquisition, and survival (Litchman and Klausmeier, 2008). Together, they affect phytoplankton fitness and are involved in several functions such as light use, nutrient uptake, or predator avoidance. Cell size can be further considered as a key trait that affects metabolism, growth rate, and community structure among others (Marañón, 2015). Zooplankton was sampled at the same station and dates as phytoplankton by means of double oblique tows from the surface to 5 m of the bottom using a 50-cm diameter Juda—Bogorov plankton net with 250-µm (until 1997) or 200-µm (from 1997 onward) mesh size. The net was equipped with a General Oceanic Flowmeter for the calculation of water filtered and a depth recorder. Samples were preserved in 2–4% sodium boratebuffered formaldehyde. Subsamples were taken to estimate total Frontiers in Marine Science | www.frontiersin.org 3November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling FIGURE 1 | Map of study area showing the location of the plankton and hydrographic sampling station (black dot), and the location of the grid cell where the upwelling index was calculated (black triangle). zooplankton abundances (in ind ×m−3) by direct examination using a stereo microscope, and biomass (in µg DW ×L−1) by weighting dried aliquots (50◦C, 48 h). Further details can be found in Bode et al. (2012). The original dataset is available at PANGAEA (https://doi.org/10.1594/PANGAEA.908815) (Bode, 2019). Hydrographic Sampling and Analysis Concurrently with the plankton samples, vertical profiles of temperature and salinity were measured with a CTD probe (Seabird SBE-25). The salinity of the CTD was checked against the salinity of bottom water samples measured with an induction salinometer Autosal 8400A calibrated with Standard Seawater (Casas et al., 1997). Salinity was expressed in the practical salinity scale (UNESCO, 1986). Nitrate concentration (NO3 in µmol L−1) was determined by segmented flow analysis according to the standard procedures of Grasshoff et al. (1983) using a Technicon AA-II (1989–2006), a Bran–Luebbe AA3 (2007–2012), and a Seal Analytics QuAAttro 39 (2013–2016). Chlorophyll aconcentration (Chl ain mg m−3) was measured by fluorimetric analysis of acetone extracts of phytoplankton collected on 0.8-µm pore-size membrane filters (until 1992) or GF/F filters (from 1993 onward). Specific calibrations were performed to ensure the continuity of the chlorophyll series when changing from the filter fluorometer method (Parsons et al., 1984) to the spectrofluorimetric technique (Neveux and Panouse, 1987) after 2001. Nutrient and chlorophyll data series are available at: https://doi.org/10.1594/PANGAEA.885413 (Bode et al., 2018). Physical Forcing Daily upwelling index (UI in m3s−1km−1) data, a rough estimate of the volume of water upwelled per km of coastline for the period 1989 to 2016 were downloaded from: http:// www.indicedeafloramiento.ieo.es/index_UI_en.html. UI was estimated from geostrophic winds calculated from the surface atmospheric pressure fields supplied every 6 h by the US Navy Operational Global Atmospheric Prediction System (NOGAPS) model maintained by the Fleet Numerical Meteorological and Oceanography Center (http://www.usno. navy.mil/FNMOC/) in a 1◦×1◦grid centered at 44◦N 9◦W (Figure 1). This cell is representative for the physical forcing Frontiers in Marine Science | www.frontiersin.org 4November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling that determines the impact of coastal upwelling in the region (Bode et al., 2015). We used the meridional wind component, thus the UI represents the volume of water upwelled along the West–East direction with positive values of UI indicating upwelling-favorable conditions. Conversely, negative values indicate downwelling-favorable conditions. Further numerical details can be found in González-Nuevo et al. (2014). In this study, values of UI were averaged over 15 days prior to each sampling date. Data Analyses Resource-Use Efficiency and Biodiversity Estimations Phytoplankton resource-use efficiency (RUEpp) was calculated sensu Ptacnik et al. (2008) in terms of phytoplankton carbon biomass (pg C L−1) per unit of nitrate concentration (µmol L−1), i.e., RUEpp =phytoplankton biomass/NO− 3concentration. Other components such as ammonia are important for the dissolved inorganic nitrogen (DIN) pool in this region; however, this variable was not measured with the same periodicity. Nonetheless, nitrate is the main limiting nutrient for the primary production in this region (Álvarez-Salgado et al., 1997). Zooplankton resource-use efficiency (RUEzp) was calculated sensu Filstrup et al. (2014) in terms of zooplankton biomass (µg L−1) per unit phytoplankton carbon biomass (pg C L−1), i.e., RUEzp =zooplankton biomass/phytoplankton biomass. Both ratios were natural log-transformed for later modeling. Taxonomic diversity (TD) was expressed as species richness (S) and as evenness (J) (Pielou, 1966) using phytoplankton biomass as follows: J=H Hmax where Hmax =lnS and H = − PS i=1Bi Btot × ln Bi Btot where Biis the biomass of a species i, and Btot is the total biomass. Additionally, we calculated various uncorrelated multitrait-based functional diversity (FD) metrics, which capture the various aspects of functional diversity (Mouchet et al., 2010). In particular, FD was expressed as functional group richness (FGR) and functional dispersion (FDis) following Laliberté and Legendre (2010), and functional evenness (FEve) following Villéger et al. (2008). To calculate the FD metrics, the abundance matrix was based on species biomass, and the functional trait matrix contained a quantitative variable (cell biomass expressed in log units) and other six binary variables describing other morphological, physiological, and behavioral traits (see Supplementary Table 1). FGR is a dendrogram-based indicator of functional groups computed from an a posteriori classification of species based on their functional traits for which the Ward method was used to create the dendrogram of the species that was cut at nine functional groups. FDis is a distancebased metric that measures the distance to the centroid of the assemblage in the trait space, and FEve is a metric that combines the evenness of species spacing in trait space and the evenness of species relative abundances. FDis and FEve metrics were weighted by the biomass of the species, and FEve was not defined for samples with fewer than three taxa (n=6 cases). Finally, single-trait-based indices, that is, communityweighted mean (CWM) traits were also calculated to examine phytoplankton single-trait contributions to RUE. CWMs were calculated for each trait weighted by species biomass using methods implemented by Laliberté and Legendre (2010). Statistical Analyses The values of RUEpp calculated for a day iand depth jwere modeled using generalized least square (GLS) models that were formulated as follows: RUEppi,j=α+ns1(DoYi)+ns2Daysi+ns3Di,j +ns4WTi,j+ns5(UIi)+ns6Si,j+ns7Ji,j+ǫi,j (1) where αis an intercept, and nsnis a natural cubic spline describing the effect of the day of the year (DoY), i.e., the seasonality, the time trend, i.e., the interannual long-term pattern (days, i.e., consecutive days from 1989 to 2016), the depth (D), the water temperature (WT), the coastal upwelling index (UI), the phytoplankton richness (S), and the phytoplankton evenness (J). All splines had 2 degrees of freedom (df) with the exception of DoY that used 4 df. Finally, ǫi,j is a vector of errors assumed to have mean 0 and variance σ2. To evaluate the variability in the response variable, i.e., a nonconstant variance, the variance in RUEpp (σ2) was modeled as an exponential function of the environmental conditions (i.e., WT or UI) or the taxonomic diversity metrics (i.e., S or J). For instance, for WT: Var(ǫi,j)= σ2×exp (2 ×δ×WTi,j), where δis an unknown parameter to be estimated that describes the estimated change in variance with water temperature. Model fitting improvement and comparison between the possible variance covariates were evaluated using the Akaike information criterion (AIC) and likelihood ratio tests (Pinheiro and Bates, 2000). Covariability among predictors was evaluated using variance inflation factors (VIFs). The values of RUEzp calculated for a day iwere also modeled using GLS models that were formulated as follows: RUEzpi=α+ns1(DoYi)+ns2Daysi+ns3(WTi) +ns4(UIi)+ns5(Si)+ns6(Ji)+ǫi(2) In this case, to match the zooplankton oblique tows, WT was averaged over the water column, and phytoplankton carbon biomass (necessary to calculate RUEzp, see above), S and J were estimated based on the species biomass averaged over the water column. The rest of the parameters and variance model are as described for Equation (1). To enable comparison of the phytoplankton and zooplankton BEF relationships using taxonomic and functional diversity, RUEpp and RUEzp calculated for a day iwere related to phytoplankton S, J, and FD metrics estimated also for a day i, therefore, based on the species biomasses averaged over the water column as explained above for the RUEzp model. To quantify the bivariate relationships, we used reduced major axis (RMA) regression, and we did not include other covariates in these models because biodiversity was the best predictor as observed in the more detailed models (see Results below). Finally, complementary analyses were further performed. In particular, we confronted the patterns found using a measure of standing stock to quantify biomass production with other Frontiers in Marine Science | www.frontiersin.org 5November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling FIGURE 2 | Partial-effects plots from the RUEpp GLS model depicted in Equation (1). Shown are the seasonal (A) and long-term (B) changes in RUEpp, the trend with depth (C), and the relationships with water temperature (D), upwelling index (E), phytoplankton richness (F), and phytoplankton evenness (G). Bands indicate 95% confidence intervals and the rugs along the x-axes display the distribution of the data. See ANOVA table in Supplementary Table 2 and model selection of variance covariates in Supplementary Table 3. Silhouettes were obtained from http://www.phylopic.org. alternative forms of calculating RUEpp, that is, calculated sensu Ptacnik et al. (2008) in terms of chlorophyll a(mg m−3) per unit of nitrate (µmol L−1), and sensu Lehtinen et al. (2017) in terms of primary production (mg C m−3h−1) per unit of nitrate (µmol L−1). Primary production was obtained from Bode et al. (2019). These models were fitted to data covering the whole period, though restricting the phytoplankton species to the 55 taxa that were consistently identified along all three decades (Supplementary Table 1). All treatment of data and analyses were performed with the software R (version 4.0.2, R Core Team, 2020) and using the packages “nlme 3.1-149” (Pinheiro and Bates, 2000), “relaimpo 2.2-3” (Grömping, 2006), “FD 1.0-12” (Laliberté and Legendre, 2010), and “lmodel2 1.7-3” (Legendre and Legendre, 1998). RESULTS Phytoplankton Resource-Use Efficiency The model fitted to phytoplankton resource-use efficiency data (Equation 1, Supplementary Table 2) revealed that RUEpp showed a seasonal cycle peaking in late March to early April (Figure 2A) and a nonlinear long-term trend over the study period (Figure 2B). RUEpp decreased linearly with sampling depth (Figure 2C) and was positively related to water temperature (Figure 2D) and upwelling index (Figure 2E). Furthermore, RUEpp was related to taxonomic diversity scaling positively with phytoplankton richness (Figure 2F) and negatively with phytoplankton evenness (Figure 2G). Studying the importance of predictors for RUEpp based on the proportional marginal variance decomposition, taxonomic diversity, and richness (61.5%) in particular, was the best predictor in explaining RUEpp dynamics, whereas the environmental factors (WT and UI) played a secondary role (Supplementary Table 2). Including a variance model in the GLS as an exponential function of covariates resulted in better fittings (Supplementary Table 3). More specifically, the spread of RUEpp decreased with evenness that was the most optimal variance covariate. In particular, an increase in 0.2 units of phytoplankton evenness reduced RUEpp variability by 21.4% (Supplementary Table 3). Finally, the RUEpp model did not show any relevant remaining patterns in the residuals (Supplementary Figure 3). Zooplankton Resource-Use Efficiency The model fitted to zooplankton resource-use efficiency data (Equation 2, Supplementary Table 4) revealed that RUEzp showed a seasonal cycle peaking in late March to early April (Figure 3A) and a decreasing trend from 2003 onward (Figure 3B). RUEzp was positively related to water temperature (Figure 3C) and had a nonlinear relationship with the upwelling index (Figure 3D). Furthermore, RUEzp was related to phytoplankton taxonomic diversity scaling negatively with richness (Figure 3E) and positively with evenness (Figure 3F). As for the case of RUEpp, when studying the importance of predictors for RUEzp, taxonomic diversity, and richness (66.6%) in particular, was the best predictor, whereas the environmental factors played a secondary role (Supplementary Table 4). Including a variance model in the GLS as an exponential function of covariates resulted in better fittings (Supplementary Table 5). More specifically, the Frontiers in Marine Science | www.frontiersin.org 6November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling FIGURE 3 | Partial-effects plots from the RUEzp GLS model depicted in Equation (2). Shown are the seasonal (A) and long-term (B) changes in RUEzp, and the relationships with water temperature (C), upwelling index (D), phytoplankton richness (E), and phytoplankton evenness (F). Bands indicate 95% confidence intervals, and the rugs along the x-axes display the distribution of the data. See ANOVA table in Supplementary Table 4, and model selection of variance covariates in Supplementary Table 5. Silhouettes were obtained from http://www.phylopic.org. spread of RUEzp decreased with richness that was the most optimal variance covariate. In particular, an increase in 10 units of phytoplankton richness reduced RUEzp variability by 18.8% (Supplementary Table 5). Finally, the RUEzp model did not show any relevant remaining patterns in the residuals (Supplementary Figure 4). Taxonomic and Functional Diversity Effects on Resource-Use Efficiency Multitrait-based functional diversity metrics were unrelated, though they showed a certain degree of association with taxonomic diversity, especially between FDis and J, and FGR and S (Supplementary Figure 5). Figures 4,5compare the effects of taxonomic and functional diversity of phytoplankton on RUEpp and RUEzp, respectively. First, both taxonomic (Figure 4A) and functional (Figure 4C) richness had positive effects on RUEpp, while the effects were negative on RUEzp (Figures 5A,C). Second, both taxonomic (Figure 4B) and functional (Figure 4D) evenness had negative effects on RUEpp, while the effects were positive on RUEzp (Figures 5B,D). Finally, the relationship with FDis showed the same trend as for FEve for phytoplankton RUE, where more functionally similar phytoplankton assemblages had higher RUE (Figure 4E), and for zooplankton RUE, where zooplankton preying upon more functionally dissimilar phytoplankton assemblages had higher RUE (Figure 5E). Regarding the explanatory power, TD metrics had higher explanatory power (i.e., greater R2) compared to FD metrics. Within TD, richness was a better predictor (48 and 47% for RUEpp and RUEzp, respectively), whereas within FD, FEve was a better predictor (29 and 20% for RUEpp and RUEzp, respectively). All slopes were statistically significant (p <0.0001), and in all cases, elevated RUEpp values occurred when diatom biomass was higher, while,when RUEzp values were higher, phytoplankton biomass was less dominated by diatoms. Finally, regarding the phytoplankton single-trait-based indices, three CWMs had important effects on plankton RUE: the ability to use biogenic silica, the ability to swim, and the ability to form chains or colonies (Figure 6,Supplementary Table 6). For RUEpp, the largest contributor was CWMchain (Figure 6C), while CWMmotility (Figure 6B) and CWMsilica (Figure 6A) contributed less to explain phytoplankton resource use. For RUEzp, all three CWMs contributed almost equally to zooplankton resource use (Figures 6D–F), though the contribution of CWMmotility was slightly higher (Supplementary Table 6). RUEpp increased with CWMsilica and CWMchain, and decreased with CWMmotility. However, the direction of the effects was the opposite for the case of RUEzp. Complementary Analyses Models using alternative calculations of RUEpp, namely, with chlorophyll a(Supplementary Figure 6) and primary production (Supplementary Figure 7) and fitted to data covering the whole period though restricting the phytoplankton species to the 55 taxa that were consistently identified along all three decades, resulted roughly in the same patterns as described above using a measure of standing stock to quantify the phytoplankton biomass production. Frontiers in Marine Science | www.frontiersin.org 7November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling FIGURE 4 | Bivariate plots showing the effects of phytoplankton taxonomic (A,B) and functional (C–E) diversity on RUEpp. Dot size is scaled to diatom biomass, and lines show reduced major axis (RMA) fits. For each relationship, the slope (and statistical significance with parametric p-value <0.0001***) and the explanatory power are shown. DISCUSSION Phytoplankton and Zooplankton Resource-Use Efficiency The model for RUEpp showed that more phytoplankton species lead to higher biomass per unit resource. This positive relationship between phytoplankton RUE and richness concurs with earlier observations made in aquatic systems both in the field (Ptacnik et al., 2008) and in experiments (Striebel et al., 2009a). At the same time, the model revealed also a negative relationship between evenness and RUEpp, which again agrees well with previous field observations in differing aquatic systems such as is the Wadden Sea, where phytoplankton evenness was the most important driver of productivity and RUE (Hodapp et al., 2015), or in Midwestern US lakes, where phytoplankton RUE was inversely related to phytoplankton evenness (Filstrup et al., 2014). On the other hand, the model for RUEzp showed that zooplankton communities produced the least biomass per unit of phytoplankton biomass when feeding upon species-rich phytoplankton communities but dominated by few species or group of species. Again, this result concurs with previous findings documenting a negative relationship between phytoplankton evenness and the production of zooplankton biomass per unit of phytoplankton biomass in lakes (Filstrup et al., 2014). However, it differs from Filstrup et al. (2019) who found a nonsignificant role of phytoplankton richness in driving zooplankton resource use efficiency in lakes. In general, two nonexclusive mechanisms have been identified to explain why biodiversity enhances ecosystem function: complementarity, by which more diverse communities use limiting resources more efficiently through niche partitioning or facilitation, and the selection effect, by which more diverse communities are more likely to include a few dominant species with specific traits that drive the ecosystem functioning (O’Connor and Byrnes, 2014). Separating and quantifying these processes can be addressed experimentally (Loreau and Hector, 2001); however, it is rather difficult to differentiate between these two mechanisms in the field. Nevertheless, the positive effect of richness found in our study can be likely explained by the niche complementarity mechanism. This interpretation would be based on the premise that a more diverse community may include more diverse traits such as light and nutrient utilization traits reflecting a better niche differentiation in wavelength utilization and resource use (Striebel et al., 2009a,b;Behl et al., 2011). This would be important during periods of nutrient enrichment (e.g., upwelling pulses) when more and more variable resources would facilitate the coexistence of a larger number of species. On the other hand, the negative effect of evenness would Frontiers in Marine Science | www.frontiersin.org 8November 2020 | Volume 7 | Article 592255 Otero et al. Phytoplankton Effect on RUE in a Coastal Upwelling FIGURE 5 | Bivariate plots showing the effects of phytoplankton taxonomic (A,B) and functional (C–E) diversity on RUEzp. Dot size is scaled to diatom biomass, and lines show RMA fits. For each relationship, the slope (and statistical significance with parametric p-value <0.0001***) and the explanatory power are shown. indicate that certain dominant species, or group of species, would exhibit specific traits that would give them certain physiological advantages allowing for a more efficient exploitation of resources, that is, a selection effect (Filstrup et al., 2019). The simultaneous action of complementarity and selection effects would probably be associated to the heterogeneous environment that typically characterizes EBUEs. Coastal upwelling systems are highly productive regions where most of the biomass is usually composed of chain-forming diatoms especially during blooms (Sarthou et al., 2005). This fact has been shown also in our study zone where diatoms are the dominant group and assemblages change rapidly between upwelling/relaxation/downwelling phases (Casas et al., 1997). Diatoms are fast growing species capable of maintaining high nutrient uptake rates for longer periods, and exploit typical upwelling-intermittent nutrient pulses more effectively than other taxa of the same size (Marañón, 2015). Dominance primarily reflects the distribution of traits within a community and the identity of the dominant traits, thus this fact has been recognized as an important effect to explain the fate of ecosystem processes because evenness often responds more rapidly to altered environmental constraints than species richness leading to rapid responses in ecosystem functioning (Hillebrand et al., 2008). Given these premises, we suggest that the mechanistic basis for the shift in the effect of phytoplankton richness and evenness would be dependent on the diatom biomass and, more specifically, on the functional traits that diatoms have. This is evidenced by the relationships found with the communityweighted mean traits highlighting that the predominance of traits that characterizes diatom species, specially the ability to form chains and colonies, increases phytoplankton resource use. Furthermore, diatoms show high growth rates relative to their cell volume. However, nutrient traits tend to be similar among taxa for a typical cell size. This indicates that diatoms appear to be adapted to high nutrient conditions as those found within upwelling regions (Edwards et al., 2012). Indeed, phytoplankton species with greater growth rates and, in particular diatoms sampled in our region tend to respond more strongly to increased upwelling (Otero et al., 2018). Concurring with our results, a population growth model fitted to phytoplankton species over an annual cycle pointed out the fundamental importance of the selection effect in driving marine primary productivity (Cermeño et al., 2016). This effect of diatom dominance would translate up in the food web. For the consumer level, we found lower RUEzp when phytoplankton biomass was dominated by diatoms, which could be likely explained by several factors including the lower impact that mesozooplankton compared to microzooplankton has on phytoplankton (e.g., Fileman and Burkill, 2001), the inhibition that diatom exudates can exert on zooplankton grazing (e.g., Frontiers in Marine Science | www.frontiersin.org 9November 2020 | Volume 7 | Article 592255