Nutrient Addition and Drought Interact to Change the Structure and Decrease the Functional Diversity of a Mediterranean Grassland

Nogueira, Carla; Nunes, Alice; Bugalho, Miguel N.; Branquinho, Cristina; McCulley, Rebecca L.; Caldeira, Maria Conceição

doi:10.3389/fevo.2018.00155

ORIGINAL RESEARCH article

Front. Ecol. Evol., 05 October 2018

Sec. Agroecology

Volume 6 - 2018 | https://doi.org/10.3389/fevo.2018.00155

Nutrient Addition and Drought Interact to Change the Structure and Decrease the Functional Diversity of a Mediterranean Grassland

$\r\nCarla Nogueira*$ Carla Nogueira¹^*

Maria Conceição Caldeira¹

¹Forest Research Centre, School of Agriculture, University of Lisbon, Lisbon, Portugal
²Centre for Ecology, Evolution and Environmental Changes, Faculdade de Ciências, Universidade de Lisboa, Lisbon, Portugal
³Centre for Applied Ecology “Prof. Baeta Neves” (CEABN-InBIO), School of Agriculture, University of Lisbon, Lisbon, Portugal
⁴Department of Plant and Soil Science, University of Kentucky, Lexington, KY, United States

Anthropogenic activities are increasing nutrient availability and altering precipitation regimes. This may lead to critical changes in grasslands functioning. This is particularly important for grasslands in the Mediterranean Basin that have evolved in nutrient poor soils, and where more frequent and prolonged droughts are projected to occur. However, there is limited knowledge regarding the interacting effects of multiple nutrient inputs and rainfall variability on the plant functional structure and diversity of Mediterranean grasslands. We conducted a nutrient addition experiment in a Mediterranean grassland during four contrasting precipitation years. We established four treatments that varied in the number of added nutrients, from no added nutrients (control), to one added nutrient (Nitrogen-N, Phosphorus-P, or Potassium-K), two added nutrients, (NP, NK, or PK) and three added nutrients (NPK). We assessed the effect of increasing nutrient addition in wet, normal, dry, and very dry years on plant species functional traits at the community level. We determined the community functional structure (e.g., Community Weighted Mean, CWM) and functional diversity (e.g., Functional Dispersion, FDis) for eight key functional traits indicators of nutrient and water use strategies. We also assessed if CWM, FDis, and species richness were related to the aggregate grassland functioning property, i.e., productivity. We found that CWM was affected by nutrient addition and precipitation and, for some traits, by their interaction. However, FDis of most traits was affected by precipitation. The very dry year had a negative effect on FDis of most traits (e.g., dispersal modes, nutrient uptake strategies) and interacted with three added nutrients to decrease FDis of growth-forms. Conversely, FDis of reproductive traits decreased during the wet year. Species richness and FDis were not related to grassland productivity, whereas CWM was the main determinant of grassland productivity supporting the importance of species functional traits in determining ecosystem functioning. Our results highlight drought as a critical factor determining a decrease in the functional diversity of Mediterranean grasslands. Moreover, drought can also interact with nutrient addition changing the dominance of many traits and further decreasing functional diversity. This may have important implications for grasslands functioning in the context of global changes.

Introduction

Mediterranean grasslands are species-diverse ecosystems of high economic and ecological value (Hector and Bagchi, 2007; Bugalho and Abreu, 2008) likely to experience biodiversity changes due to interacting global change drivers (Sala et al., 2000).

Anthropogenic activities, such as fossil fuel combustion and fertilizer applications are increasing nutrient inputs, such as of nitrogen (N) and phosphorus (P), into the biosphere (Peñuelas et al., 2013). Nutrient addition often leads to decreases in plant species diversity and compositional shifts, which may potentially affect ecosystem functioning and resilience (Bobbink et al., 2010; Stevens et al., 2010; Flores-Moreno et al., 2016; Harpole et al., 2016). Climate change may also lead to reduction in species richness (e.g., Thuiller et al., 2005) and induce changes in community composition (e.g., Dornelas et al., 2014). Moreover, these two drivers of change, may interact, further affecting grassland diversity and functioning (e.g., Zavaleta et al., 2003; Harpole et al., 2007). This is particularly important in the Mediterranean region, a hot-spot for climate changes (Giorgi, 2006), where increased precipitation variability and frequency of droughts is forecasted to occur (Costa et al., 2012; Kovats et al., 2014).

Many studies have evaluated the responses of grassland plant species composition and diversity to alterations in environmental conditions, such as fertilization or drought (e.g., Tilman and El Haddi, 1992; Zavaleta et al., 2003; Harpole et al., 2007, 2016; Bobbink et al., 2010; Stevens et al., 2010; DeMalach et al., 2017). Species richness and diversity, per se, however do not reflect the role of plant species in ecosystem functioning (e.g., Tilman et al., 1997; McGill et al., 2006; Cadotte et al., 2011; Mouillot et al., 2011). Species functional traits, however, can provide a mechanistic and predictive understanding of ecosystem responses to environmental changes (e.g., Díaz et al., 2007; Mouillot et al., 2011). Functional traits are morphological, physiological, or phenological features of species that can be measured at the individual level (Díaz and Cabido, 2001). Functional traits determine the species functional role in ecosystems and/or their response to changing environmental conditions or disturbances (Díaz and Cabido, 2001). Ecosystem functioning and resilience depend on the functional traits of the dominant plant species in a community (i.e., community functional structure), in accordance with the mass-ratio hypothesis (Grime, 1998), and on the dissimilarity of functional trait values within that community (here defined as functional diversity sensu Laliberté and Legendre, 2010). Communities with a high functional diversity are expected to improve ecosystem functioning by increased complementarity in resource use among species (Tilman et al., 1997) and to provide high resilience to disturbances (Díaz et al., 2007; Cadotte et al., 2011; Mouillot et al., 2011; Volaire et al., 2014; Valencia et al., 2015).

Although the response of communities' functional structure and diversity to decreased nutrient limitation (e.g., Suding et al., 2005; Janeček et al., 2013; Helsen et al., 2014) or to drought (e.g., Polley et al., 2013; Carmona et al., 2015; Nunes et al., 2017) have been previously assessed, the concurrent effects of both remain unclear. Nutrient addition was shown to favor the abundance of species with high specific leaf area (SLA) and height, probably because these traits correlate better with competition for nutrients and light (Westoby et al., 2002; Suding et al., 2005; Ordoñez et al., 2009; Janeček et al., 2013; Eskelinen and Harrison, 2015). This may explain a positive response in the abundance of grasses (Zavaleta et al., 2003) to the detriment of the less competitive, rosette growth-forms (Pfestorf et al., 2013) when nutrients were added. The increased frequency of droughts predicted for the Mediterranean Basin, coupled with its high inter-annual precipitation variability, may also modify the functional structure and diversity of grassland communities, and affect ecosystem functioning (Pérez-Camacho et al., 2012; Carmona et al., 2015; Nunes et al., 2017). For example, drought generally favors small stress-tolerant and slow-growing species with low SLA, high leaf dry matter content (LDMC) (e.g., Westoby et al., 2002; Carmona et al., 2015), delayed flowering onset and shorter flowering duration (Crimmins et al., 2013), although an earlier flowering onset has also been observed with increasing drought (Pérez-Camacho et al., 2012; Shavrukov et al., 2017). In addition, Rota et al. (2017) observed a decrease in seed mass functional diversity due to drought, whereas functional diversity of SLA was unresponsive.

These findings indicate that nutrient addition and drought can select for different plant traits underlying opposite strategies, such as high and low SLA, respectively, thus drawing attention to the importance of understanding the response of grassland communities to interacting global change drivers (Bernard-Verdier et al., 2012; Pérez-Ramos et al., 2013; Eskelinen and Harrison, 2015; Rota et al., 2017).

In this 4-year study, we assessed how increasing nutrient availability and different precipitation years, ranging from wet to very dry years, interacted to affect community functional structure and functional trait diversity of a Mediterranean grassland. To achieve this, we focused on eight key plant functional traits known to be responsive to nutrient and water use strategies (Cornelissen et al., 2003; Garnier et al., 2007) and used the community weighted mean (CWM) and functional dispersion (FDis) as metrics of community functional structure and functional diversity, respectively. We hypothesized that (1) nutrient addition and precipitation are strong interacting filters of plant community functional structure; and (2) as precipitation plays a fundamental role in determining richness and community composition in Mediterranean ecosystems, functional diversity will be primarily influenced by precipitation.

Materials and Methods

Site Description

We conducted our study in a semi-natural Mediterranean grassland, located north-east of Lisbon, Portugal (38°49′45.13^′′N, 8°47′28.61^′′W). The grassland community, occasionally grazed by cattle until 2012, was dominantly composed by annual C3 species that grow in winter and early spring, and start senescing in late spring. Overall, dominant species are the forbs Tolpis barbata L., Plantago bellardii All., the graminoids Agrostis pourretii Willd., Avena barbata Link, and the nitrogen-fixing legumes Ornithopus compressus L. and Trifolium arvense L.

The climate is Mediterranean-type with long-term (1961–1990) mean annual rainfall of 709 mm, with considerable inter-annual variation. Mean annual temperature is 15.9°C, ranging from 10°C in January to 22.5°C in August [(Instituto Nacional de Meteorologia e Geofísica (INMG), 1991)]. Site topography is flat. The soil is a well-drained deep Haplic Arenosol (IUSS Working Groups WRB, 2006) with a low soil water retention capacity. Prior to treatment initiation, two soil cores (2.5 cm in diameter, 10 cm depth from soil surface) were collected from each plot. The plot subsamples were composited, homogenized, and air-dried. Soil C (%) and N (%) were determined by The Ecosystems Analysis Laboratory at the University of Nebraska using dry combustion GC analysis (COSTECH ESC 4010 Elemental Analyzer, Costech Analytical Technologies, Valencia, California, USA). Extractable soil P and K and soil pH were determined at A&L Analytical Laboratory (Memphis, TN). Soil pH was measured using a 1:1 soil to water slurry (Table 1).

TABLE 1

Table 1. Pre-treatment nutrient soil content and pH values.

Experimental Design

We added N, P, and potassium (K) combined with a micronutrient mixture to 5 × 5 m plots, following a full factorial combination in a complete randomized three block design. We considered four treatments that varied in the number of added macronutrients (0, 1, 2, and 3) following Harpole et al. (2016). Hereafter, we use the term “number of added nutrients” for simplicity. Specifically, the treatments were: Controls, i.e., no added nutrients, termed as 0 (6 replicates); additions of one single nutrient N, P, or K (9 replicates); addition of two join nutrients (NP, NK, or PK), named as 2 added nutrients (9 replicates); and plots fertilized with a combination of 3 nutrients, NPK that is 3 added nutrients treatment (6 replicates) (Table 2). The experiment is part of the global Nutrient Network experiment (www.nutnet.org; Borer et al., 2017)

TABLE 2

Table 2. Nutrient addition treatments grouped as number of added nutrients following Harpole et al. (2016).

The experiment started in 2012 and lasted until 2017. However, we excluded the year of 2015 as the cattle broke through the fence during February and heavily defoliated and trampled the plots. We added nutrients every fall, at the beginning of the growing season. We applied N, P, and K at a rate of 100 kg.ha⁻¹.yr⁻¹. N was added from 2012 to 2016 as slow-release urea (60–90 days), switching in 2017 to two separate additions of urea (3 months apart to mimic slow release) due to restricted availability of timed-release urea. P was added as triple-super phosphate, and K as potassium sulfate. A micronutrient mixture (6% Ca, 3% Mg, 12% S, 0.1% B, 1% Cu, 17% Fe, 2.5% Mn, 0.05% Mo, and 1% Zn) was added at a rate of 1,000 kg.ha⁻¹ with K and was only applied in 2012, the first year of the study, to avoid micronutrient toxicity.

The site was fenced to avoid intense cattle grazing and trampling. We used a brush cutter and removed all standing vegetation in the fenced experimental area at the end of each growing season after seed production and dispersal to eliminate potential effects of carried-over biomass (Dudney et al., 2017).

In situ precipitation was recorded using a tipping bucket rain gauge (RG2, Delta-T Devices, Cambridge, UK), connected to a CR1000 and AM16/32B multiplexer data logger (Campbell Scientific, Logan, USA).

Precipitation Years

During our study, the hydrological precipitation (October–September) of the four years was 483 mm (very dry year, 2017), 598 mm (dry year, 2016), 708 mm (normal year, 2013), and 814 mm (wet year, 2014) (Figure 1A). The year of 2017 was the driest hydrological year, with precipitation 32% below the long-term average (709 mm). April 2017 was unusually dry with 28 consecutive days without rain (Figure 1B) and unusually hot, from April to June (average for the 3 months: 3.2°C, above to long-term) (Figure 1C). Precipitation in 2016 was 16% below and in 2014, 15% above the long-term average (Figure 1A). The year of 2015, which was not included in our analysis, was also a very dry year, especially during winter and spring, with total precipitation 465 mm, 34% below the long-term average (Figure 1A).

FIGURE 1

Figure 1. (A) Precipitation difference from the long-term hydrological year average (709 mm, 1961–1990), for 2013–2017. 2013 (normal precipitation year—total precipitation, 708 mm), 2014 (wet year−814 mm), 2015 (very dry year−465 mm) not included in our study, 2016 (dry year−598 mm) and 2017 (very dry year−483 mm). Monthly values of total precipitation (B) and average air temperature (C) for 2013, 2014, 2016, 2017 and long-term (1961–1990) for comparison.

Vegetation Sampling

Species cover (%), to the nearest 1%, of each plant species was estimated annually by observation, using a modified Daubenmire method (Daubenmire, 1959), at peak vegetation growth, in a permanently marked 1 m² quadrat in each plot.

Species Richness and Aboveground Biomass

We have also examined the pattern of plant species richness and grassland productivity to increasing number of added nutrients (Nogueira et al., unpublished). Species richness was determined as the number of plant species per plot each year. Aboveground plant biomass was estimated by collecting, every year, all individual plants rooted within two 10 × 100 cm strips per plot, at the time of peak biomass accumulation. We considered aboveground biomass as a proxy for productivity. The locations of the clipping strips were moved within each plot every year to prevent clipping effects. Samples were dried to constant mass at 60°C prior to weighing.

Trait Selection

Selection of plant functional traits is a critical step to understand and predict community functioning and responses to environmental changes and disturbance (Petchey and Gaston, 2006). Therefore, we selected eight plant traits including continuous, ordinal and categorical traits, that reflect plant functional responses to water and soil nutrient availability (e.g., Cornelissen et al., 2003; Table 3). Trait data was obtained for 61 species identified in our grassland communities throughout the 4-year study. Trait information was obtained through direct observation or measurements in the field, following standard protocols (Cornelissen et al., 2003; Pérez-Harguindeguy et al., 2013), and/or derived from the literature (Franco, 1971, 1984; Castroviejo, 1986–2012; Valdés et al., 1987; Garnier et al., 1997; Paula and Pausas, 2013; Table 3). We used a single mean trait value per species to compute functional trait metrics. If different trait values for the same species could be found in the literature, we opted for plant trait values obtained in Mediterranean regions. Although intraspecific trait variability plays an important role in many ecosystem processes and contributes to the understanding of community response to environmental changes (Lepš et al., 2011), in this study we only consider differences in functional structure and functional diversity caused by changes in species composition (Bernard-Verdier et al., 2012; Carmona et al., 2015; Rota et al., 2017).

TABLE 3

Table 3. Functional traits studied, their typology/units and their main role in plant functioning.

Functional Trait Metrics

To assess community functional structure, we determined the Community Weighted Mean (CWM) for each trait, which corresponds to the average trait value in a community, weighted by the relative abundance of the species carrying each trait value (Garnier et al., 2007). CWM trait values are a quantitative translation of the biomass ratio hypothesis (Grime, 1998) by reflecting the dominant trait values in a community. For continuous traits CWM values represent the mean value of that trait in the community, for categorical and binary traits CWM values correspond to the proportion of each category in the community.

Functional trait diversity within the plant community was determined through Functional Dispersion (FDis), which is closely related to Rao's quadratic entropy (Laliberté and Legendre, 2010). It estimates the dispersion of species in functional trait space, weighted by their relative abundances. It is calculated as the weighted mean distance, in multidimensional trait space, of individual species from the weighted centroid of all species, where weights correspond to species relative abundances (Laliberté and Legendre, 2010). The minimum value FDis can take is 0 (e.g., communities composed by only one species) and it has no upper limit. FDis has several advantages over other functional diversity indices as it takes into account species relative abundances, it is unaffected by species richness, handles any number and type of traits (including more traits than species) and is not strongly influenced by outliers (Laliberté and Legendre, 2010). We calculated FDis for each functional trait individually and for the six traits combined, that responded significantly to nutrient addition and precipitation (multi-trait). The calculation of trait metrics was done using the FD package (Laliberté and Shipley, 2011), in software R (R Core Team, 2015).

Statistical Analyses

To examine the effect of number of added nutrients and precipitation (hydrological year) on community plant species composition, we performed a non-metric multidimensional scaling ordination (NMDS) based on species cover of each plot sampled in each nutrient addition treatment throughout our four-year study. We used meta MDS function from vegan package (Oksanen et al., 2013), with Bray-Curtis distance measure. Ordination stress statistic was used as a measure of goodness of fit. A permutational multivariate analysis of variance (PERMANOVA) was performed with 9,999 permutations on Bray-Curtis communities' dissimilarities, to test for the effect of number of added nutrients and precipitation year (Oksanen, 2011) on species composition.

The influence of number of added nutrients, precipitation year, and their interaction on single-trait CWM and FDis values and on multi-trait FDis computed for traits that were significantly affected by nutrient addition and precipitation years was examined by means of general linear mixed models. Number of added nutrients and precipitation year were used as fixed effects explanatory variables and plots were used as a random factor, to account for non-independent observations made along the four years of the study in each plot. Tuckey pairwise comparisons with Bonferroni corrections were performed to check for differences between number of added nutrients and between different precipitation years.

We used Spearman correlation to test for the significance of correlations between species richness, CWM, FDis, and productivity. We selected growth-form, maximum height, SLA, and LDMC for Spearman correlation as these traits may reflect grassland productivity (Table 3). Spearman correlation was used to account for possible nonlinear relationships.

All the analyses were performed under R statistical environment (R Core Team, 2015). Graphs, except Figure 2, were done with Sigmaplot 13.0 (Systac Software, Inc., San Jose, USA).

FIGURE 2

Figure 2. Non-metric multidimensional scaling (NMDS) based on species cover, with points representing plots sampled in each nutrient addition treatment in different hydrological precipitation years: normal year (2013−708 mm), wet year (2014−814 mm), dry year (2016−598 mm), and very dry year (2017−483 mm). Points representing plots are united according to (A) the number of added nutrients; (B) hydrological year. Final stress for the 2-dimensional configuration was 0.22.

Results

The Effect on Community Species Composition

Main gradients in community species composition are described by a 2-dimensional non-metric multidimensional scaling (NMDS) ordination with a final stress value of 0.22 (Figures 2A,B). The first axis (NMDS1) explained 35.5% and the second axis (NMDS2) explained 18.6% of the variation in species composition (Figures 2A,B). Plant community composition changed significantly with the increase in number of added nutrients (PERMANOVA test: F = 9.98, P < 0.001), and with the total hydrological precipitation, mainly in the very dry year (PERMANOVA test: F = 14.01, P < 0.001).

Functional Structure of Plant Communities

Ten out of sixteen CWM trait values were significantly related to number of added nutrients (Table 4), four of which had a significant interaction with precipitation year (graminoid and erect growth-forms, SLA, and LDMC; Table 4). Differences between nutrient addition treatments were, in most cases, only significant between three added nutrients and controls and one added nutrient treatment. In general, CWM values were significantly affected by the very dry year compared to other type of years (Table 4).

TABLE 4

Table 4. Statistical results of the linear mixed models for CWM for each functional trait.

Species Growth-Forms and Height

Overall, dominant species in our grassland were graminoids (proportional cover, average ± s.e.m.: 0.34 ± 0.02), erect (0.32 ± 0.02), and rosette (0.21 ± 0.02), while ascending (0.03 ± 0.003) and prostrate species (0.10 ± 0.01) were less abundant.

Increasing number of added nutrients had a significant positive effect on the cover of graminoids to the detriment of ascending, erect, and rosette growth-forms (Table 4). Prostrate species were significantly and positively affected by precipitation and unaffected by nutrient addition (Table 4). In the very dry year, prostrate species were residual (0.01 ± 0.003) as compared to other years (0.13 ± 0.01, average across treatments).

Graminoids were affected by a significant interaction between the number of added nutrients and precipitation (Table 4). During the wet year, abundance of graminoids was lowest (0.17 ± 0.02) compared to other years (0.40 ± 0.03); with no significant differences between nutrient addition treatments (Figure 3A). However, for other precipitation years, graminoids were significantly more abundant in the three added nutrients treatment as compared to one added nutrient (normal: t = −3.89, P = 0.004; dry: t = −3.48, P = 0.011; very dry: t = −5.73, P < 0.001) and controls (normal: t = −4.26, P = 0.001; dry: t = −4.09, P = 0.002; very dry: t = −4.89, P < 0.001; Table 4; Figure 3A). Abundance of rosette or ascending growth-forms were significantly and negatively affected by number of added nutrients and positively affected by precipitation (Table 4). Contrary to graminoids, rosette species were less abundant in three nutrients treatment (0.10 ± 0.021) compared to one nutrient (0.26 ± 0.03; z = −4.64, P < 0.001), two nutrients (0.19 ± 0.03; z = −2.76, P = 0.035), and marginally to controls (0.25 ± 0.02; z = −2.50, P = 0.076). Rosette species were also, marginally significantly more abundant during the wet year (0.30 ± 0.03) compared to the very dry year (0.19 ± 0.04) (z = 2.50, P = 0.074). During the very dry year, abundance of erect species was significantly lower in three nutrients treatment (0.05 ± 0.03) compared to one nutrient (0.29 ± 0.09; t = 3.73, P = 0.006) and controls (0.43 ± 0.16; t = 4.30, P = 0.001; Table 4).

FIGURE 3

Figure 3. CWM values for (A) graminoid growth-form; (B) height; (C) SLA; (D) onset of flowering; and (E) N-fixing ability, along a gradient of number of added nutrients, for a normal precipitation year (2013−708 mm), wet year (2014−814 mm), dry year (2016−598 mm), and very dry year (2017−483 mm). Symbols represent mean ± s.e.m, for number of added nutrients for each precipitation year. For one and two added nutrients n = 9; for controls and three added nutrients n = 6 for each precipitation year. Lines show mean response between CWM trait values and number of added nutrients for each precipitation year.

We found significantly higher maximum height in communities where three nutrients were added (average species height: 86.93 cm ± 3.19) compared to one nutrient treatment (64.32 cm ± 2.28; z = 3.27, P = 0.007) and controls (61.13 cm ± 3.73; z = 4.33, P < 0.001), with a significant positive precipitation effect (Table 4; Figure 3B).

Species Leaf Traits

SLA was the only trait with no significant association with precipitation (Table 4; Figure 3C). However, the nutrient addition treatments only significantly affected CWM of SLA during the very dry year, with the three added nutrients exhibiting higher SLA (29.45 mm².mg⁻¹ ± 2.46) compared to one nutrient treatment (24.11 mm².mg⁻¹ ± 1.35) (z = 3.30, P = 0.002) and controls (23.15 mm².mg⁻¹ ± 2.03) (z = 3.63, P = 0.006) (Table 4; Figure 3C). LDMC was also higher in three nutrients (301.83 mg.g⁻¹ ± 15.31) compared to one nutrient treatment (227.32 mg.g⁻¹ ± 20.86) (t = −4.13, P = 0.002) and to controls (236.62 mg.g⁻¹ ± 16.57) (t = −3.16, P = 0.024) during the very dry year. Nonetheless, precipitation had an overall significant negative effect on LDMC (Table 4).

Species Reproductive Strategies

Across all treatments and precipitation years, onset of flowering occurred between February and April, and was significantly affected by both the number of added nutrients and precipitation (Table 4; Figure 3D). Overall our grassland community responded with an earlier onset of flowering with increasing number of added nutrients. Communities had later flowering during the very dry year compared to normal (t = −2.84, P = 0.027) and wet years (t = −2.84, P = 0.061). Flowering duration, which varied between 3 and 6 months, was shorter during the very dry year compared to the wet year (z = 4.01, P < 0.001), normal year (z = 4.83, P < 0.001), and dry year (z = 3.72, P = 0.001). Also, only in the very dry year communities with three added nutrients had shorter flowering duration than one (t = 3.15, P = 0.025) and two nutrients treatments (t = 3.14, P = 0.025).

Species Dispersal Strategies

Across all nutrient treatments and precipitation years, our grassland community was dominated by species with wind dispersal (0.71 ± 0.02), followed by species with dispersal by gravity (0.26 ± 0.01); and rarely by internal (0.02 ± 0.003) or external animal transport (0.01 ± 0.002). Dispersal modes were unaffected by number of added nutrients (Table 4). Abundance of species with wind dispersal was highest during the very dry year (0.89 ± 0.02) compared to the, wet (0.63 ± 0.03; z = −3.92, P < 0.001), normal (0.71 ± 0.03; z = −4.16, P < 0.001), and the dry year (0.62 ± 0.03; z = −5.33, P < 0.001), while species with dispersal by gravity or internal animal transport were significantly less abundant (Table 4).

Species Nutrient Uptake Strategies

Most species in our grassland were non N-fixers (0.83 ± 0.015, average across treatments and years). Abundance of N-fixing species was not significantly affected by number of added nutrients (Table 4), and represented 0.14 ± 0.015 (average across treatments and years) of our grassland community. N-fixing species were significantly negatively affected by very low precipitation, compared to the wet year (z = 6.60; P < 0.001), normal year (z = 5.82, P < 0.001), and dry year (z = −7.64, P < 0.001) (Figure 3E). During the very dry year, N-fixing species decreased by 94%, from 0.18 ± 0.017 (3-year average) to 0.01 ± 0.003 (very dry year). Although, abundance of hemi-parasites species was low in our grassland, it was further decreased with the addition of three nutrients treatment (0.01 ± 0.005) compared to controls (0.09 ± 0.012), during the normal (t = 4.92, P < 0.001) and the dry years (t = 5.92, P < 0.001).

Functional Dispersion of Traits (FDis)

The functional dispersion of growth-form was significantly negatively affected by nutrient addition (Table 5), but only in the very dry year, with three nutrients having the lowest diversity in growth-form compared to the other treatments: two added nutrients (t = 5.91, P < 0.001); one added nutrients (t = 7.11, P < 0.001) and controls (t = 4.12, P = 0.002) (Figure 4A). FDis of growth-form was significantly lower during the very dry year compared to other years (wet: z = 4.27; P < 0.001, normal: z = 5.13, P < 0.001, and dry year: z = 4.21, P < 0.001).

TABLE 5

Table 5. Statistical results of the linear mixed models analyzing the effects of nutrient addition and precipitation on functional dispersion (FDis) for each functional trait and for a multi-trait combination of six relevant traits.

FIGURE 4

Figure 4. FDis values for (A) growth-form; (B) height; (C) onset of flowering; (D) nutrient uptake strategy; and (E) multi-trait, the combined six traits affected by nutrient addition and precipitation, along a gradient of number of added nutrients, for a normal precipitation year (2013−708 mm), wet year (2014−814 mm), dry year (2016−598 mm), and very dry year (2017−483 mm). Symbols represent mean ± s.e.m, for number of added nutrients for each precipitation year. For one and two added nutrients n = 9; for controls and three added nutrients n = 6 for each precipitation year. Line shows mean response between FDis values for each trait/ multi-trait and number of added nutrients, for each precipitation year.

FDis of most traits were not affected by increasing number of added nutrients but they were significantly affected by precipitation (Table 5). FDis of height, increased overall with precipitation, whereas, FDis of the onset of flowering and flowering duration decreased with precipitation (Table 5; Figures 4B,C). FDis of dispersal mode was significantly lower during the very dry year (0.15 ± 0.02) compared to the wet (0.30 ± 0.01; z = 4.41, P < 0.001), normal (0.27 ± 0.01; z = 4.48, P < 0.001), and the dry year (0.31 ± 0.01; z = 5.04, P < 0.001). FDis of nutrient uptake strategy followed the same trend and was significantly lower during the very dry year compared to the wet (z = 8.15, P < 0.001), normal (z = 9.29, P < 0.001), and the dry year (z = 10.07, P < 0.001) (Table 5; Figure 4D). FDis of SLA and LDMC were unaffected by nutrient addition and precipitation (Table 5). The FDis of the combined six traits (multi-trait FDis) decreased significantly during the very dry year in all nutrient treatments compared to the remaining years (wet: z = 5.99, P < 0.001; normal: z = 6.75, P < 0.001, and dry year: z = 7.05, P < 0.001) but it was not affected by the number of added nutrients (Table 5; Figure 4E).

Functional Traits and Productivity

Increasing number of added nutrients significantly decreased species richness and significantly increased productivity as our grassland is co-limited by multiple nutrients (Nogueira et al., unpublished). However, species richness had no significant relationship with productivity (ρ = −0.046, P = 0.621; n = 120), whereas most of the CWM trait values were significantly related to productivity (Table 6). Higher abundance of graminoid and prostrate species were significantly positively related to productivity, as were taller species, or/and higher LDMC. Contrarily, high abundance of species with an erect or rosette growth form were significantly negatively related to productivity (Table 6). FDis of the analyzed traits had no relation with productivity (Table 6).

TABLE 6

Table 6. Spearman correlations (ρ) between CWM (community weighted mean) and functional dispersion (FDis) of each trait and productivity (n = 120; for all treatments and for the four years).

Discussion

Overall, our results showed a combined effect of nutrients and precipitation on the community functional trait structure (CWM), while functional diversity (FDis) was mainly affected by precipitation. Nutrient addition only decreased functional diversity of growth-forms in the very dry year. A clear change in species composition was observed during the very dry year and in three nutrients addition treatment, which was reflected by marked changes in community functional structure and trait diversity.

Precipitation was the major environmental filter determining CWM and FDis. This is in agreement with the high inter-annual precipitation variability and low soil fertility context in which Mediterranean annual grasslands have coevolved (Noy-Meir, 1973; Peco et al., 1998; Pérez-Camacho et al., 2012); (Carmona et al., 2015; Rota et al., 2017). Functioning and resilience of the Mediterranean grasslands is probably mostly based on a high functional diversity of traits responsive to different levels of soil water availability combined with resource conservative traits (Pérez-Camacho et al., 2012; Carmona et al., 2015). The decrease in co-limitation by water and nutrients increased the dominance of traits associated with resource-acquisition such as maximum height (Westoby et al., 2002; Eskelinen and Harrison, 2015) but not SLA. Graminoids responded positively to nutrient addition as found in other studies (Hautier et al., 2009; Spasojevic and Suding, 2012; Helsen et al., 2014; Niu et al., 2014; Harpole et al., 2016) but negatively to precipitation (but see Harpole et al., 2007). Graminoids have higher root cation exchange capacity compared to most dicots (Woodward et al., 1984; DiTommaso and Aarssen, 1989) and higher root density (Craine et al., 2001) which makes them highly competitive for available nutrients. Interestingly, in the wet year, a low abundance year for graminoids, the functional group did not respond to nutrient addition. Indeed, Collins et al. (2012) observed that irrigation increased the abundance of forb species. On the contrary, during the very dry year, graminoids increased by two-fold comparing to other years, and were better competitors for nutrients. This was also supported by the increase in SLA and LDMC in the three added nutrient communities in the very dry year. In accordance, Wellstein et al. (2017) in a Mediterranean grassland, found that graminoids had a higher SLA under experimental drought and demonstrated better growth performance, which was most likely related to their strategy to allocate resources to roots. Additionally, in our study, the FDis of growth-forms was lowest in the very dry year, particularly in communities with lower co-limitation by nutrients. The niche dimension hypothesis (Harpole and Tilman, 2007) states that the decrease in the number of limiting resources should lead to loss of diversity. However, in our study it was during the very dry year, with a high limitation by soil water, that a decrease in FDis of growth-forms, dispersal modes, nutrient uptake strategies, as well as FDis of the main six functional traits combined, was observed. Contrastingly, during the wet year, with the highest soil water availability, only the FDis of onset and duration of flowering decreased. Our results indicate that the mechanisms underlying functional diversity loss can be dependent on the environmental and evolutionary context.

In our study, the increase in species height and SLA in three added nutrients communities probably resulted from the increase of graminoids which were taller and had higher SLA than rosette and erect growth-forms. This pattern was also observed by Ansquer et al. (2009) in grassland communities in the Pyrenees. Dominance of graminoids with increased nutrients may also help to explain the earlier flowering in nutrient rich communities, except in the very dry year. Later flowering during the very dry year is somewhat counterintuitive. If abundance of graminoids increased during the very dry year, onset of flowering should have been earlier. Also, early flowering is one of the many plant strategies to cope with increased water limitation (e.g., Shavrukov et al., 2017). However, in our study, Vulpia bromoides was the dominant graminoid species during the very dry year in all treatments. This species has a later onset and a shorter flowering duration which may explain main reproductive strategy observed during that year. During a wet year, we found a reduction in FDis of reproductive traits, which suggests that trait filtering, may also occur under more favorable conditions (Bernard-Verdier et al., 2012).

Species in our grassland were mostly dispersed by wind and gravity. Although nutrient addition had no effect on observed plant dispersal modes, drought affected wind dispersal positively and gravity dispersal negatively. The higher abundance of graminoids, all wind-dispersal species and the significant reduction of N-fixing species, most of them gravity-dispersal species, during the very dry year explains this result. In Mediterranean environments with low and highly fluctuating precipitation, small or large seeds can be both optimal (Volis and Bohrer, 2013). However, the larger seed production and greater longevity in the seed bank associated with small seeded species can be advantageous with increased aridity (Volis and Bohrer, 2013).

The very dry year had a negative effect on the abundance of N-fixing species as observed in other studies (e.g., Espigares and Peco, 1993, 1995; Del Pozo and Aronson, 2000). Also, N-fixing species were found to be highly vulnerable to drought during their reproductive phase, which may shorten the duration of reproductive development, reducing seed number and weight (Daryanto et al., 2015). As N-fixing species have, generally, a positive effect on the nitrogen budget of grasslands (Spehn et al., 2002), the projected increase in the frequency and intensity of droughts in the Mediterranean region, can negatively affect the nitrogen budget. Although in our study there was no association between added nutrients and abundance of N-fixing species, we found a negative correlation between N added treatments and abundance of N-fixing species, in accordance with other studies (DiTommaso and Aarssen, 1989; Suding et al., 2005; Silvertown et al., 2006).

Although species richness may be a useful indicator of ecosystem functioning (Tilman et al., 1997; Hector et al., 1999; Hector and Bagchi, 2007), in our study CWM but not species richness was correlated with productivity, an aggregate property reflecting ecosystem functioning. In line with the “mass ratio hypothesis” (Grime, 1998), we found that graminoids, an abundant growth-form in our grassland, contributed significantly to aboveground biomass. However, less abundant traits such as prostrate growth-forms were also positively related to productivity, suggesting that even the loss of low abundance traits may be of great importance to ecosystem functioning (McLaren and Turkington, 2010). Thus, global change drivers that have a significant impact on species composition (Dornelas et al., 2014) will most probably have a strong effect on ecosystem processes (Tilman et al., 1997). Beyond functional structure, functional diversity has also been found to be associated to ecosystem processes (Tilman et al., 1997; Valencia et al., 2015). In more productive systems, a few abundant species with similar traits (low functional diversity) usually dominate communities, whereas in low productive systems there is frequently a higher functional diversity (niche complementarity) (e.g., Tilman et al., 1997; Rolo et al., 2016). However, in our grassland we found no relation between functional diversity and productivity, suggesting that under strong environmental gradients (precipitation variability and nutrient availability), community structure was the main determinant of productivity as observed for other grasslands (Zhu et al., 2016).

Conclusions

Our results revealed an important shift in the mean community functional trait composition of our Mediterranean grassland, with decreased water and nutrient co-limitation. Precipitation was the main driver of functional diversity, affecting most of the studied traits. However, other factors such as air temperature can interact with total and distribution of precipitation to affect grassland composition (Peco et al., 1998; Zavaleta et al., 2003; Miranda et al., 2009; Nogueira et al., 2017) and thus the functional attributes of the community. We found, that nutrient addition only decreased growth-form FDis during the very dry year, highlighting the importance of the interacting effects of global change drivers on grassland functional diversity. These results will have strong implication to Mediterranean grassland dynamics mainly within global change contexts. Indeed, higher fertilization and increased drought frequency are likely to affect grassland functional diversity. We show that interacting global change drivers may further constraint trait variability in Mediterranean grasslands. Moreover, higher drought frequency will likely favor graminoids, ultimately affecting the diversity of the seed bank which may decrease, resulting in lower resilience and resistance of these ecosystems to global change (Volaire et al., 2014). This should be further tested using long-term grassland studies (e.g., Crawley et al., 2005) coupled to climate manipulation. Our work further supports that species functional traits may be more predictive of ecosystem processes than species richness making it an important tool to understand ecosystem responses to global changes.

Data Availability Statement

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

CN and MC conceived the idea and designed methodology. CN collected the data. AN and CN analyzed the data. CN and MC led the writing of the manuscript. MB, CB, and RM contributed critically to the draft. All authors gave final approval for publication.

Funding

The authors acknowledge the financial support of the Portuguese Science Foundation (Fundação para a Ciência e Tecnologia- FCT), through the project IF/740/2014 and PTDC/AAG-GLO/0045/2014 and fellowships to CN (SFRH/BD/88650/2012) and MB (IF/1171/2014). Centro de Estudos Florestais (CEF) a research unit funded by FCT, Portugal (UID/AGR/00239/2013). We are also grateful to FLAD/NSF for funding (A1/Proj 124/12).

Conflict of Interest Statement

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

The authors acknowledge the Nutrient Network (www.nutnet.org) and Companhia das Lezírias for providing access to the grassland site and permission to undertake research. We thank J. Mendes, A. Horta, P. Pais, L. Marçal, and M. Silva for field and laboratory assistance. We also thank the reviewers for comments on an earlier version of the manuscript.

References

Ansquer, P., Duru, M., Theau, J. P., and Cruz, P. (2009). Convergence in plant traits between species within grassland communities simplifies their monitoring. Ecol. Indic. 9, 1020–1029. doi: 10.1016/j.ecolind.2008.12.002