Ecosystem Functional Attributes (EFAs)

We identified three EFAs that capture most of the variance in the time series of vegetation indices and with a meaningful biological significance (Paruelo et al., 2001; Alcaraz-Segura et al., 2006). These attributes were calculated from the EVI seasonal curve or annual dynamics: the annual mean (EVI_mean; an estimator of primary production), the EVI seasonal coefficient of variation (EVI_sCV; a descriptor of seasonality), and the date of maximum EVI (EVI_DMAX; an indicator of phenology). To summarize the EFAs of 2001-2018, we calculated the mean of the period for each attribute.

Ecosystem Functional Types (EFTs)

EFTs were identified following Alcaraz-Segura et al., (2013), from the three EFAs obtained. The range of values of each EFA was divided into four intervals, giving a potential number of 64 EFTs (4 × 4 × 4). For EVI_DMAX, the four intervals agreed with the four seasons of the year. For EVI_mean and EVI_sCV, we extracted the first, second, and third quartiles for each year and then calculated the quartiles means for the 16-year period. To name EFTs, we used two letters and a number: the first capital letter indicates net primary production (EVI_mean), increasing from A to D; the second small letter represents seasonality (EVI_sCV), decreasing from a to d; the numbers are a phenological indicator of the growing season (EVI_DMAX), with values 1-spring, 2-summer, 3-autumn, 4-winter. To summarize the EFTs of the 2001–2018 period, we calculated the dominant EFT (i.e., the mode value for each pixel) of the period.

Ecosystem Functional Diversity

In order to characterize ecosystem functional diversity we used EFT richness and EFT rarity as indicators:

• EFT richness map was calculated for each year by counting the number of different EFTs within a 4×4-pixel moving window (924 x 924 m; ~ 1 km^2) across the study area (Alcaraz-Segura et al., 2013). Then, an average richness map across all years was obtained.

• EFT rarity map was calculated for each year as the relative extension of each EFT compared to the most abundant EFT (Cabello et al., 2013). Then, we calculated an average rarity map across the period.

\[\begin{equation} \label{Eq:1} \textbf{Rarity of } EFT_{i} = \frac{\text{Area } EFT_{max} - \text{Area } EFT_{i}} {\text{Area }EFT_{max}} \end{equation}\]

where \(\text{Area } EFT_{max}\) is the area occupied by the most abundant EFT and \(\text{Area } EFT_{i}\) is the area of the \(i\) EFT being evaluated, with \(i\) ranging from 1 to 64.

Stability in ecosystem functioning

In order to show which are the most stable areas and with the greatest inter-annual variability (either due to directional changes or fluctuations) in ecosystem functioning, the number of different EFTs that occurred in the same pixel in the period 2001-2018.

To obtain a measure equal to the interannual variability but at landscape level, the dissimilarity, i.e., Jaccard’s index 1-coefficient, was calculated. Dissimilarity values range from 0 to 1, with 1 being the highest degree of dissimilarity in composition and relative abundance of EFTs and 0 being absent.