Phytoplankton Diversity Effect on Ecosystem Functioning in a Coastal Upwelling System

Otero, Jaime; Álvarez-Salgado, Xosé Antón; Bode, Antonio

doi:10.3389/fmars.2020.592255

ORIGINAL RESEARCH article

Front. Mar. Sci., 26 November 2020

Sec. Marine Ecosystem Ecology

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.592255

Phytoplankton Diversity Effect on Ecosystem Functioning in a Coastal Upwelling System

$\nJaime Otero &#x;$ Jaime Otero¹^*^†

Xosé Antón Álvarez-Salgado¹

Antonio Bode²

¹Instituto de Investigaciones Marinas (IIM-CSIC), Vigo, Spain
²Centro 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.

Introduction

Marine phytoplankton is responsible for roughly half of the global primary production (50 Pg C year⁻¹; 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 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 hump-shaped 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 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 long-term 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).

FIGURE 1

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).

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 a concentration 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 μm³) 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⁻¹) of a given species in a given day and sampled depth was then calculated as the abundance (in cells L⁻¹) 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 a concentration (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 borate-buffered formaldehyde. Subsamples were taken to estimate total zooplankton abundances (in ind × m⁻³) by direct examination using a stereo microscope, and biomass (in μg DW × L⁻¹) 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 (NO₃ in μmol L⁻¹) 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 a concentration (Chl a in mg m⁻³) 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 m³ s⁻¹ km⁻¹) 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 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 (RUE_pp) was calculated sensu Ptacnik et al. (2008) in terms of phytoplankton carbon biomass (pg C L⁻¹) per unit of nitrate concentration (μmol L⁻¹), i.e., RUE_pp = phytoplankton biomass/ ${NO}_{3}^{-}$ concentration. 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 (RUE_zp) was calculated sensu Filstrup et al. (2014) in terms of zooplankton biomass (μg L⁻¹) per unit phytoplankton carbon biomass (pg C L⁻¹), i.e., RUE_zp = 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 = \frac{H}{H_{m a x}}$ where H_max = lnS and $H = - \sum_{i = 1}^{S} \frac{B_{i}}{B_{t o t}} \times l n \frac{B_{i}}{B_{t o t}}$ where B_i is the biomass of a species i, and B_tot 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 distance-based 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, community-weighted 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 RUE_pp calculated for a day i and depth j were modeled using generalized least square (GLS) models that were formulated as follows:

\begin{array}{l} {RUEpp}_{i, j} = α + n s_{1} ({DoY}_{i}) + n s_{2} ({Days}_{i}) + n s_{3} (D_{i, j}) \\ + n s_{4} ({WT}_{i, j}) + n s_{5} ({UI}_{i}) + n s_{6} (S_{i, j}) + n s_{7} (J_{i, j}) + ϵ_{i, j} & (1) \end{array}

where α is an intercept, and ns_n is 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 σ². To evaluate the variability in the response variable, i.e., a nonconstant variance, the variance in RUE_pp (σ²) 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) = σ² × exp (2 × δ × WT_i,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 RUE_zp calculated for a day i were also modeled using GLS models that were formulated as follows:

\begin{array}{l} {RUEzp}_{i} = α + n s_{1} ({DoY}_{i}) + n s_{2} ({Days}_{i}) + n s_{3} ({WT}_{i}) \\ + n s_{4} ({UI}_{i}) + n s_{5} (S_{i}) + n s_{6} (J_{i}) + ϵ_{i} & (2) \end{array}

In this case, to match the zooplankton oblique tows, WT was averaged over the water column, and phytoplankton carbon biomass (necessary to calculate RUE_zp, 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, RUE_pp and RUE_zp calculated for a day i were 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 RUE_zp 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 alternative forms of calculating RUE_pp, that is, calculated sensu Ptacnik et al. (2008) in terms of chlorophyll a (mg m⁻³) per unit of nitrate (μmol L⁻¹), and sensu Lehtinen et al. (2017) in terms of primary production (mg C m⁻³ h⁻¹) per unit of nitrate (μmol L⁻¹). 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 RUE_pp 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). RUE_pp decreased linearly with sampling depth (Figure 2C) and was positively related to water temperature (Figure 2D) and upwelling index (Figure 2E). Furthermore, RUE_pp 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 RUE_pp based on the proportional marginal variance decomposition, taxonomic diversity, and richness (61.5%) in particular, was the best predictor in explaining RUE_pp 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 RUE_pp decreased with evenness that was the most optimal variance covariate. In particular, an increase in 0.2 units of phytoplankton evenness reduced RUE_pp variability by 21.4% (Supplementary Table 3). Finally, the RUE_pp model did not show any relevant remaining patterns in the residuals (Supplementary Figure 3).

FIGURE 2

Figure 2. Partial-effects plots from the RUE_pp GLS model depicted in Equation (1). Shown are the seasonal (A) and long-term (B) changes in RUE_pp, 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.

Zooplankton Resource-Use Efficiency

The model fitted to zooplankton resource-use efficiency data (Equation 2, Supplementary Table 4) revealed that RUE_zp showed a seasonal cycle peaking in late March to early April (Figure 3A) and a decreasing trend from 2003 onward (Figure 3B). RUE_zp was positively related to water temperature (Figure 3C) and had a nonlinear relationship with the upwelling index (Figure 3D). Furthermore, RUE_zp was related to phytoplankton taxonomic diversity scaling negatively with richness (Figure 3E) and positively with evenness (Figure 3F). As for the case of RUE_pp, when studying the importance of predictors for RUE_zp, 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 spread of RUE_zp decreased with richness that was the most optimal variance covariate. In particular, an increase in 10 units of phytoplankton richness reduced RUE_zp variability by 18.8% (Supplementary Table 5). Finally, the RUE_zp model did not show any relevant remaining patterns in the residuals (Supplementary Figure 4).

FIGURE 3

Figure 3. Partial-effects plots from the RUE_zp GLS model depicted in Equation (2). Shown are the seasonal (A) and long-term (B) changes in RUE_zp, 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.

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, 5 compare the effects of taxonomic and functional diversity of phytoplankton on RUE_pp and RUE_zp, respectively. First, both taxonomic (Figure 4A) and functional (Figure 4C) richness had positive effects on RUE_pp, while the effects were negative on RUE_zp (Figures 5A,C). Second, both taxonomic (Figure 4B) and functional (Figure 4D) evenness had negative effects on RUE_pp, while the effects were positive on RUE_zp (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 R²) compared to FD metrics. Within TD, richness was a better predictor (48 and 47% for RUE_pp and RUE_zp, respectively), whereas within FD, FEve was a better predictor (29 and 20% for RUE_pp and RUE_zp, respectively). All slopes were statistically significant (p < 0.0001), and in all cases, elevated RUE_pp values occurred when diatom biomass was higher, while,when RUE_zp values were higher, phytoplankton biomass was less dominated by diatoms.

FIGURE 4

Figure 4. Bivariate plots showing the effects of phytoplankton taxonomic (A,B) and functional (C–E) diversity on RUE_pp. 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.

FIGURE 5

Figure 5. Bivariate plots showing the effects of phytoplankton taxonomic (A,B) and functional (C–E) diversity on RUE_zp. 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.

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 RUE_pp, the largest contributor was CWM_chain (Figure 6C), while CWM_motility (Figure 6B) and CWM_silica (Figure 6A) contributed less to explain phytoplankton resource use. For RUE_zp, all three CWMs contributed almost equally to zooplankton resource use (Figures 6D–F), though the contribution of CWM_motility was slightly higher (Supplementary Table 6). RUE_pp increased with CWM_silica and CWM_chain, and decreased with CWM_motility. However, the direction of the effects was the opposite for the case of RUE_zp.

FIGURE 6

Figure 6. Scatterplots showing the relationships between RUE_pp (A–C) and RUE_zp (D–F) with CWM_silica (A,D), CWM_motility (B,E), and CWM_chain (C,F). Lines show RMA fits, and corresponding numerical data are provided in Supplementary Table 6.

Complementary Analyses

Models using alternative calculations of RUE_pp, 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.

Discussion

Phytoplankton and Zooplankton Resource-Use Efficiency

The model for RUE_pp 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 RUE_pp, 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 RUE_zp 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 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 community-weighted 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 RUE_zp 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., Malej and Harris, 1993), or the reduced impact that copepod grazing has on phytoplankton biomass and production during production peaks compared to upwelling transition periods when phytoplankton biomass is lower, and dinoflagellates are more abundant (e.g., Bode et al., 2003). Indeed, when there was a predominance of traits in the community that characterize dinoflagellate species (e.g., ability to swim), zooplankton resource use was higher. Likewise, in lakes, when zooplankton preys upon phytoplankton communities dominated by cyanobacteria, RUE_zp decreases because cyanobacteria is a poor food source for zooplankton (Filstrup et al., 2014).

A simultaneous detection of a positive and negative effect of richness and evenness, respectively, on ecosystem function as the one observed here, was shown also by others either in fresh (Filstrup et al., 2019) or marine (Napoleón et al., 2014) waters. However, there are contrasting results on the relative contribution of the two effects, that is, the differences between the strength of the relationships. Overall, our models indicated that taxonomic richness consistently had much stronger effects on both phytoplankton and zooplankton RUE. This result contrasts with what has been observed in lakes, where evenness explained more variance in phytoplankton and zooplankton RUE than richness leading the authors to suggest that it might be the distribution of taxa rather than their number, the main driver of ecosystem function in lakes (Filstrup et al., 2019). Coastal upwelling ecosystems are highly heterogeneous and contain an enormous species richness rapidly changing with the environmental conditions. Such high number of species is probably necessary to maintain the overall productivity of the system and ensures the occurrence of highly productive species. This scenario would lead to both complementarity and selection effects operating simultaneously, with complementarity being stronger when functionally different species coexist (e.g., during the spring bloom) and take advantage of a heterogeneous resource supply, whereas selection effects would be important when a broad trait space coincides with a more homogeneous landscape (Hodapp et al., 2016). Nonetheless, taken altogether, our results suggest that the opposing effects of primary producer biodiversity on phytoplankton and zooplankton resource use efficiency would be of similar nature in the ocean and freshwater environment. However, there appear to be differences on the sensitivity of ecosystem functioning to richness and evenness that could be ascribed to intrinsic differences among systems, such as changes in hydrographic and biogeochemical characteristics, or the identity of the dominant taxa, or to other underlying mechanisms demanding further research.

Environmental Effects on Resource-Use Efficiency

Apart from the importance of the environmental conditions in determining the overall composition of the community outlined above, we detected direct abiotic effects on RUE and time-varying patterns. In particular, our statistical models showed an annual peak located in April for both phytoplankton and zooplankton RUE coinciding with the spring bloom (Casas et al., 1997). Moreover, there were somehow an interannual inverse pattern between RUE_pp and RUE_zp that would point to changes in the plankton community along years. In this regard, the studied system displayed phase shifts in plankton community structure around the turn of the 21st century following changes in climate and local hydrography (Bode et al., 2020). The models also revealed the influence of water temperature and upwelling intensity on RUE, albeit less pronounced than the effects of biodiversity. These effects were weak for RUE_zp and especially apparent for RUE_pp. Resource use efficiency had greater values when upwelling intensity was high and waters were warmer. This combination of high upwelling during the 15 days previous to the sampling date and high water temperature during the sampling date indicates the typical succession of intense upwelling followed by wind relaxation. The positive effect of upwelling contributes to fertilize the system, which promotes a proliferation of phytoplankton biomass and a more efficient use of resources during the subsequent upwelling relaxation (e.g., Huete-Ortega et al., 2010). The effect of temperature could also be explained because this variable favors cell growth. Besides that, RUE_pp decreased with depth likely associated to the decrease in phytoplankton biomass in deeper and less illuminated waters. Therefore, the abiotic environment was found to play a role in explaining part of the variability of RUE; however, its effect was minor compared to those of biodiversity. Recent meta-analyses showed that the effects of abiotic and biotic drivers in mediating ecosystem properties can overlap or be even stronger than the effects of species loss in controlled experiments (Godbold, 2012). However, Ptacnik et al. (2008) showed in the field that the relative importance of abiotic drivers vs. that of richness was minor in explaining phytoplankton RUE dynamics, thus agreeing with our results. Moreover, Hodapp et al. (2015) identified also a weaker effect of temperature and light compared to evenness on phytoplankton productivity.

Biodiversity Effects on Resource-Use Efficiency Variability

Our statistical models also showed that phytoplankton and zooplankton RUE were both less variable at higher levels of diversity, evenness in the case of RUE_pp, and richness in the case of RUE_zp. These results indicate that RUE_pp and RUE_zp would be stabilized when phytoplankton biomass is more evenly distributed, and when zooplankton preys upon richer phytoplankton communities, respectively. In general, theoretical, experimental, and field studies conclude that temporal stability of ecosystem processes would be greater at higher diversity, usually in the form of species numbers, and that those processes would be simultaneously enhanced also at higher levels of diversity (Tilman et al., 2014). Our variance models do not explicitly account for the temporal scale, but we can interpret the fact that biodiversity increases RUE stability as a consequence of that; first, more diverse assemblages are more likely to include different species maintaining the ecosystem function under different environmental conditions (i.e., the insurance hypothesis), and second, more diverse assemblages may show increased asynchrony in species' abundances (i.e., the portfolio effect) favoring the stability of RUE (Thibaut and Connolly, 2013). Furthermore, contrary to expectations, our results showed that the highest levels of RUE were not associated with the lowest levels of variability. This would be in line with recent synthesis of multiple experiments, which concluded that, while biodiversity of primary producers simultaneously increases both the production and stability of biomass in terrestrial and aquatic ecosystems, these effects are independent, suggesting that an ecosystem process and its variability would not necessarily reach the highest levels simultaneously with biodiversity (Cardinale et al., 2013).

Taxonomic vs. Functional Diversity Effects on Resource-Use Efficiency

Our analyses included the evaluation of the performance of three FD metrics resulting in the same patterns as their TD counterparts. This is in line with other works that pointed out the validity and potential of using phytoplankton functional traits to explain community assembly and ecosystem processes both in fresh (e.g., Leruste et al., 2018) and marine (e.g., Breton et al., 2017) waters. However, our results showed that FD metrics were slightly poorer predictors of RUE than TD metrics. This observation differs from what has been previously shown by other comparative studies undertaken in fresh (Abonyi et al., 2018) and marine (Ye et al., 2019) waters where FD usually outperforms TD in predicting ecosystem functioning. This discrepancy could be explained in part because the set of traits included in our analysis is insufficient and does not fully capture the fine variability along the species-specific ecological niches. For instance, quantitative light utilization traits and nutrient utilization traits are determinant for species-specific metabolism; thus, both functions, efficient use of light and nitrogen acquisition, might not be well represented reducing the percentage of variance of RUE explained by FD. Functional traits are typically sourced from published laboratory experiments, usually not available for all studied species (as would be the case here) but for the most common species, thus preventing a thorough examination of FD effects on RUE for this particular upwelling ecosystem. Apart from the comparison in explanatory power between TD and FD, it is worth mentioning that, in contrast to TD, FD evenness outperformed FD richness, which was a poor predictor, suggesting that FEve captures better the distribution of species across the functional trait space and its effect on RUE. Furthermore, the observed effects of FDis highlight that an increase in functional similarity of phytoplankton increases phytoplankton production per unit of nitrate, whereas at the same time, zooplankton biomass per phytoplankton biomass was higher when phytoplankton assemblages were more functionally dissimilar. This would reinforce the hypotheses that a selection effect would occur when trait dissimilarity was low in diatom-dominated assemblages inducing greater phytoplankton RUE and lower zooplankton RUE. In line with this, Cadotte (2017) using experimental plant assemblages found opposing relationships between selection (negative) and complementarity (positive) effects with functional dispersion.

Study Limitations

Our study has some limitations. On the one hand, we have not accounted for potential cell size changes through time that might bias biomass calculations as size is expected to decrease when temperature increases (Morán et al., 2010). On the other hand, and more importantly, we acknowledge that we have worked with an incomplete number of species. It is known that conventional sampling methods may underestimate phytoplankton diversity (Cermeño et al., 2014); however, we used a very specific assemblage that made the dataset more homogeneous along the time series. This assemblage was indeed representative as shown by the relationship between chlorophyll a and phytoplankton carbon, and covered the most abundant species occurring in this region in each recognized upwelling/relaxation/downwelling phase (Casas et al., 1997). Furthermore, this assemblage was able to explain the variability in phytoplankton and zooplankton RUE calculated in terms of chlorophyll a. Additionally, our definition of RUE was based on the quantification of biomass production using a measure of standing stock thus not accounting for gross production; nonetheless, using a direct measure of primary production (secondary production was not available), we also found comparable patterns in the dynamics of RUE_pp. These complementary analyses suggest that the patterns found are robust. Therefore, we believe that, given the consistency of the sampling methods along the three decades, and the comparable results obtained using different RUE calculations, the facts outlined above could have some influence on the effect size of our models, for instance, by minimizing the effect of rare species, though not on the direction of the effects. Finally, we assumed here a single resource limitation (nitrate); thus, further investigations are needed to explore the importance of colimitation or multiple resource limitation (Hodapp et al., 2019) following the dynamics of upwelling.

Conclusions

We found that taxonomic and functional diversity enhance planktonic resource use efficiency in a coastal upwelling system. The sign of this relationship differed with the components of diversity considered (richness vs. evenness) and had opposite patterns across trophic levels. These effects can be likely explained through a simultaneous complementarity effect, by which more diverse assemblages with more diverse functional traits use resources more efficiently, and a selection effect through the dominance of diatoms growing rapidly and exploiting typical upwelling-intermittent nutrient pulses more effectively than other taxa of similar cell size. Additionally, resource use efficiency was less variable at higher levels of diversity; however, greater stability was not associated with the highest levels of RUE. Our results contribute to the understanding of phytoplankton diversity effects on ecosystem processes underpinning the differing effects of two contrasting biodiversity metrics and how these relationships vary across adjacent trophic levels in a coastal upwelling area. Incorporation of this approach to the observation, experimentation, and modeling of plankton ecology in EBUEs will allow gaining understanding on the important processes occurring in these productive coastal ecosystems and how these processes are mediated by biodiversity.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.1594/PANGAEA.908815; https://doi.org/10.1594/PANGAEA.885413.

Author Contributions

JO designed the project. AB provided the data. JO analyzed the data with assists from AB and XÁ-S. JO wrote the manuscript with comments from co-authors. All authors contributed to the article and approved the submitted version.

Funding

This work is part of the time series project RADIALES conducted and funded by the Instituto Español de Oceanografía (http://www.seriestemporales-ieo.net) with additional support from project MarRisk (Interreg POCTEP Spain-Portugal) grant number 0262 MARRISK 1 E, from grants Contrato-Programa GAIN-IEO, and grant number IN607A2018/2 of the Axencia Galega de Innovación (GAIN, Xunta de Galicia, Spain). JO was supported by a Junta para la Ampliación de Estudios Fellowship (JAE-Doc programme 2011) from the CSIC and ESF.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This article is dedicated to the memory of Dr. Manuel Varela, respected friend and colleague, who passed away on October 24, 2019. We would like to acknowledge the dedication of a large number of technicians, crew members, and scientists who contributed to the observational time series project RADIALES. We are particularly indebted to M. Varela (phytoplankton) and M. T. Álvarez-Ossorio (zooplankton) for maintaining the plankton series for more than 20 years. J. Lorenzo continued the phytoplankton counts since 2011, and E. Rey and M. A. Louro ensured the continuity of the zooplankton series. Nutrient data were provided by N. González and R. Carballo (1989–2012) and by M. Álvarez and M. Castaño (2013–2016).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.592255/full#supplementary-material

References

Abonyi, A., Horváth, Z., and Ptacnik, R. (2018). Functional richness outperforms taxonomic richness in predicting ecosystem functioning in natural phytoplankton communities. Fresh. Biol. 63, 178–186. doi: 10.1111/fwb.13051