Understanding the mechanisms underlying community stability has become an urgent need in order to protect ecosystems from global change and resulting biodiversity loss. While community stability can be influenced by richness, synchrony in annual fluctuations of species, species stability and functional traits, the relative contributions of these drivers to stability are still unclear. In semi-natural grasslands, land-use changes such as fertilization might affect stability by decreasing richness and influencing year-to-year fluctuations. In addition, they can promote long-term directional trends, shifting community composition and influencing grassland maintenance. Thus, it is important to consider how species and community stability vary year-to-year but also in the long term.

Using a 14-year vegetation time series of a species-rich semi-natural Mediterranean grassland, we studied the relative importance of richness, synchrony, species stability and functional traits on community stability. To assess land-use change effects on stability, we applied a fertilization treatment. To distinguish stability patterns produced by year-to-year fluctuations from those caused by long-term trends, we compared the results obtained using a detrending approach from those without detrending.

Stability is influenced by richness, synchrony and functional traits. Fertilization decreases species and community stability by promoting long-term trends in species composition, favouring competitive species and decreasing richness. Studying stability at the community and species level, and accounting for the effect of trends is essential to understand stability and its drivers more comprehensively.

Study site and experimental design

In 2003, 12 plots of 15x5 m (hereafter called macro-plots) were established inside an area of 5500 m². Half of the macro-plots (six) were used as control plots and half were fertilized with sewage sludge in a single event in 2003, applying manually to the soil surface 5 kg m^-2 . The sludge came from a municipal urban wastewater treatment plant located in Tudela (Navarra, Spain), and it was sludge previously dried to 28% dry matter by centrifugation. To accurately assess vegetation changes, a 1x1 m permanent plot was placed in the centre of each macro-plot. Every year for 14 consecutive years (from 2004 to 2017), at the end of June, vegetation was sampled by R. Ibáñez, who identified and recorded every vascular plant species present in each of the permanent plots. The 1x1 m permanent plots were divided into 100 10x10 cm subplots to measure species abundance (frequency) by counting the number of 10x10 cm subplots in which the species was present (presence was recorded if shoots overlapped with the sampling unit/subplot, not according to rooted plants).

Richness, synchrony and stability measures

Species richness in each permanent plot was measured both as cumulative species richness, counting the number of species found at least once in a permanent plot during the 14 years of the study, and as mean species richness, averaging the number of species found in a permanent plot over the 14 years (Lepš et al., 2018). Community-level synchrony for each permanent plot was calculated using the log variance ratio index (“Logvar”), which is the log-transformation of the ratio of observed to expected variance (i.e. the ratio of variance of the total community abundance to the sum of variances of the abundance of each species; Lepš et al., 2018; Roscher et al., 2011). Stability at both the community and the species level was calculated as the inverse of the coefficient of variation (CV^-1) across years of cumulative or individual species abundances in each permanent plot. In order to distinguish the patterns produced by long-term trends from those caused by year-to-year fluctuations, we also used the three-term local quadrat variance detrending method (T3; Hill, 1973), which consists in calculating the variance in three year time periods, to remove the effect of long-term trends both on synchrony and stability indices (Lepš, Götzenberger, et al., 2019). Thus, we calculated synchrony (log variance ratio) and stability (CV^-1) indices using both the non-detrending (hereafter “long-term” synchrony and stability) and the T3-detrending approach (hereafter “year-to-year” synchrony and stability; for which the original variance used in the log variance ratio index or in the inverse of the coefficient of variation was replaced by the three-term local-quadrat variance) (Valencia, de Bello, Lepš, et al., 2020).  All synchrony and stability indices were calculated using the calc_sync function of the package “tempo” in R (Lepš, Götzenberger, et al., 2019).

Plant functional traits and indices

We obtained data for five functional traits (plant height, Leaf Dry Matter Content or LDMC, Specific Leaf Area or SLA, Leaf Area or LA and Seed Mass or SM) for most of the species present in the 1x1 m permanent plots (data available for 98%, 78%, 84%, 98% and 85% of the species, respectively). Trait data were collected in-situ. Although the number of individuals measured for each trait varied, traits were measured following the protocols provided by Cornelissen et al. (2003). Missing data were obtained from BROT, a trait database for Mediterranean Basin species (Tavşanoğlu & Pausas, 2018), or from TRY database (Kattge et al., 2020).

Mean trait values and abundance data for each species were used to calculate community functional composition and functional diversity for the five traits studied. Functional composition was measured as the Community Weighted Mean (CWM; Garnier et al., 2004) and functional diversity as the Rao Quadratic Index (Rao, 1982), which is the sum of pairwise functional distances between species weighted by relative abundance (Mouchet et al., 2010). To calculate these indices, we used the function dbFD of the “FD” package in R (Laliberté & Legendre, 2010; Laliberté et al., 2014).

Data analysis

To study how fertilization influenced stability, synchrony and species richness at the community level, and to test if the values displayed by these indices and their response to fertilization changed depending on the approach used to calculate them (e.g. long-term or year-to-year), we carried out three multiple linear regression models using as response variables the stability, synchrony or species richness index in each plot, respectively. As explanatory variables we used the treatment applied in each plot (i.e. control or fertilized), the approach used to calculate the index (i.e. long-term or year-to-year stability and synchrony, and cumulative or mean species richness), and the interaction between both.

In order to discover the main drivers of community stability, we used different linear regression models to test for relationships between community stability (long-term and year-to-year) and synchrony, richness, functional composition and diversity of the five traits studied and treatment. The explanatory variables were calculated in different ways depending on if we studied long-term or year-to-year stability, so that the variables would as well reflect “accumulated” or “yearly” values. For long-term stability we used long-term synchrony and cumulative species richness, and the functional composition and diversity indices were calculated using the cumulative abundance of each species across all years. By contrast, for year-to-year stability we used year-to-year synchrony and mean species richness across years, and the functional composition and diversity were calculated separately for each plot and year and then averaged across all years. This way, the first type of models was more focused on detecting processes acting across years and promoting trends, while the second type was focused on year-to-year processes. We build linear regression models by first running simple regression models for each explanatory variable and then applied multiple regression models with synchrony, richness, treatment and the functional indices selected as significant or marginally significant in the simple models.

We then studied stability at the species level. Species stability was averaged over control and over fertilized plots. As we found some extreme values corresponding to highly stable species (Brachypodium retusum (Pers.) P.Beauv. and Aphyllanthes monspeliensis L.), we log transformed average species stability. We tested for relationships between log-transformed species stability (long-term and year-to-year) and the five functional traits studied. We first used simple models for each trait and when a certain trait was significant we did multiple models using trait, treatment and the interaction between both as explanatory variables. On the other hand, we also studied how fertilization influenced species stability, checking if results changed when using the long-term or year-to-year approach. We calculated the Pearson’s correlation coefficient and applied a paired t-test to test for differences in species stability between control and fertilized plots and between the long-term and year-to-year approach. All analyses were carried out with the lm, cor.test and t.test functions in R software (v. 4.0.3; R Core Team, 2020).

datasheet.xlsx

Data used in the analyses, consists on two sheets:    
- boxplotindices: data with stability, synchrony and richness values for each plot or community, the approach used to calculate each index (long-term or year-to-year), and the treatment applied in each plot (control or fertilized).
    
- matrixalldatacom: data used to carry out the analyses at the community level. For each community it is given: the number of the plot or community, stability, synchrony, richness, functional composition and diversity of the five traits studied, and treatment (control or fertilized). For each of these variables there are two columns, one with the variable calculated using the long-term approach (named "longterm_" or "cum_") and the other one for the variable calculated using the year-to-year approach (named "yeartoyear_" or "mean_").

- matrixalldataspp: data used to carry out the analyses at the species level. For each species it is given: species name, long-term and year-to-year stability, species values for the five traits studied, and treatment (control or fertilized). Missing values are indicated by NAs.

indicescom.R: calculation of richness, synchrony and stability indices at the community level.

functionalindicescom.R: calculation of functional composition (Community Weighted Mean) and diversity (Rao) for the five traits studied (plant height, Leaf Dry Matter Content, Specific Leaf Area, Leaf Area, Seed Mass), at the community level.

spstabilityindices.R: calculation of stability indices at the species level. 

models_boxplots.R: R script for multiple linear regression models analysing the relationship between richness, synchrony or stability values and the approach used (long-term or year-to-year), treatment (control or fertilized) and the interaction between approach and treatment.
    
models_stability_communitylevel.R: R script for simple and multiple linear regression models analysing the drivers of long-term and year-to-year community stability.
    
SEMs.R: R script for Piecewise Structural Equation Models analysing the relationships between long-term and year-to-year community stability and their drivers.
    
models_stability_splevel.R: R script for simple and multiple linear regression models analysing the relationship between long-term and year-to-year species stability and functional traits. 
    
spstability.R: R script for Pearson's correlations and paired t-tests analysing the effects of treatment and long-term trends on species stability.

README.txt: text file describing the data.