A remote sensing-based dataset to characterize the ecosystem functioning and functional diversity in the Biosphere Reserve of Sierra Nevada (SE Spain)

Conservation Biology faces the challenge of safeguarding the ecosystem functions and ecological processes (the water cycle, nutrients, energy flow, and community dynamics) that sustain the multiple facets of biodiversity. Characterization and evaluation of these processes and functions can be carried out through functional attributes or traits related to the exchanges of matter and energy between vegetation and the 20 atmosphere. Based on this principle, satellite imagery can provide integrative spatiotemporal characterizations of ecosystem functions at local to global scales. Here, we provide a multi-temporal dataset at protected area level, that characterizes the spatial patterns and temporal dynamics of ecosystem functioning in the Biosphere Reserve of Sierra Nevada (Spain), captured through the spectral Enhanced Vegetation Index (EVI, using product MOD13Q1.006 from MODIS sensor) from 2001 to 2018. The 25 database contains, at the annual scale


Introduction
A better characterization of the functional dimension of biodiversity is required to develop management approaches that ensure nature contributions to human well-being (Jax, 2010;Cabello et al., 2012;Bennet et al., 2015).A set of essential variables for characterizing and monitoring ecosystem functioning is necessary to achieve this goal (Pereira et al., 2013).Such variables are critical to understanding the dynamics of ecological systems (Petchey and Gaston, 2006), the links between biological diversity and ecosystem services (Balvanera et al., 2006;Haines-Young and Potschin, 2010), and the mechanisms of ecological resilience (Mouchet et al., 2010).In addition, there have been calls for the use of ecosystem functioning variables to assess functional diversity at large scales to measure biosphere integrity (Mace et al., 2014;Steffen et al., 2015), one of the most challenging planetary boundaries to assess (Steffen et al., 2015).Despite the importance of ecosystem functioning variables and the conceptual frameworks developed to promote their use (Pettorelli et al., 2018), they have seldom been incorporated into ecosystem monitoring (see Alcaraz-Segura et al., 2009;Fernández et al., 2010;Cabello et al., 2016;Skidmore et al., 2021).
Characterization and evaluation of ecosystem functioning can be carried out through monitoring functional traits or attributes related to the matter and energy exchanges between biota and the atmosphere (Box et al., 1989;Running et al., 2000).Nowadays, the use of satellite imagery provides useful methods to derive such attributes, allowing for the spatially explicit characterization of ecosystem functioning and its heterogeneity (i.e., functional diversity) from local (Fernández et al., 2010) to regional (Alcaraz-Segura et al., 2006, 2013), and global scales (Ivits et al., 2013;Skidmore et al., 2021).Theoretical and empirical models enable essential functional variables of ecosystems, such as primary production, evapotranspiration, surface temperature, or albedo to be estimated based on spectral indices derived from satellite images (e.g., Enhanced Vegetation Index; EVI) (Pettorelli et al., 2005;Fernández et al., 2010;Lee et al., 2013).Among them, primary production is one of the most integrative and essential descriptor of ecosystem functioning (Virginia and Wall, 2001;Pereira et al., 2013) since it has a significant role in the carbon cycle (i.e., the energy input to the trophic web and the driving force behind many ecological processes) (Xiao et al., 2019).Moreover, primary production provides a holistic response to environmental changes and constitutes a synthetic indicator of ecosystem health (Costanza et al., 1992;Skidmore et al., 2015).
To characterize ecosystem functioning focusing on primary production, we can use the satellite-derived approach based on Ecosystem Functional Types (EFTs), defined as patches of the land surface that share similar dynamics in the matter and energy exchanges between the biota and the physical environment (Paruelo et al., 2001;Alcaraz-Segura et al., 2006).EFTs are derived from three Ecosystem Functional Attributes (EFAs) that describe the seasonal dynamics of carbon gains through time-series of spectral vegetation indices: 1) EVI annual mean (a surrogate of annual primary production), 2) EVI annual standard deviation (a descriptor of seasonality or the differences between the growing and non-growing seasons), and 3) the annual date of maximum EVI (a phenological indicator of the date within a year on which the growing period is centered).Biologically, these three metrics can be interpreted as surrogates of the total amount and timing of primary production (Paruelo et al., 2001;Pettorelli et al., 2005;Alcaraz-Segura et al., 2006), one of the most integrative indicators of ecosystem functioning (Virginia and Wall, 2001).Since the EFT concept appeared in 2001 (Paruelo et al., 2001), its applicability has grown to characterize functional heterogeneity from local to global scales (Alcaraz-Segura et al., 2006;Karlsen et al., 2006;Duro et al., 2007;Fernández et al., 2010;Geerken 2009;Alcaraz-Segura et al., 2013;Ivits et al., 2013;Cabello et al., 2013;Pérez-Hoyos et al., 2014;Müller et al., 2014;Wang and Huang, 2015;Villarreal et al., 2018;Coops et al., 2018;Mucina, 2019;Cazorla et al. 2020b).
Here, we present a dataset that describes the spatial heterogeneity and temporal variability of a key measure of ecosystem functioning and ecosystem functional diversity patterns in the Sierra Nevada Biosphere Reserve, a protected area in south-east Spain (Fig. 1).We derived the dataset from the analysis of the intraand inter-annual variation of vegetation greenness captured through the EVI spectral vegetation index, as a surrogate of primary production, during the 2001-2018 period.First, for each year, we provide maps of three EFAs: i) annual primary production, and ii) seasonality and iii) phenology of carbon gains, as well as their integration into a synthetic mapping of EFTs as discrete landscape functional units.Second, based on these units, we present two functional diversity metrics: EFT richness and EFT rarity.Finally, by considering the yearly maps, we calculated inter-annual summaries, i.e., inter-annual means and interannual variability, to show the average conditions and the stability of ecosystem functioning of the period (the workflow is provided in Fig. 2).The abundance of long-term datasets from multiple disciplines in Sierra Nevada constitutes an opportunity to explore the role of ecosystem functioning and functional diversity on ecohydrological and species distribution modeling, climate change mitigation and adaptation, ecological resilience, adaptive management, and mountain ecosystem services supply since humankind depends on freshwater of mountain regions, are hotspots of biodiversity and a key destination for recreation activities (Grêt-Regamey et al., 2012).

Data acquisition and processing
We based ecosystem functioning characterization on analyzing the temporal dynamics of the EVI from 2001 to 2018.We chose EVI instead of any other vegetation index (such as SAVI, ARVI, or NDVI) as an indicator of carbon gains since it is more reliable in both low and high vegetation cover situations (Huete et al., 1997).Furthermore, EVI reduces the influence of atmospheric conditions on vegetation index values and corrects for canopy background signals.
EVI is computed following this equation (Equation 1): Equation 1, where NIR/RED/Blue are atmospherically corrected (Rayleigh and ozone absorption) surface reflectance for near infrared, red and blue wavelengths, respectively; L is the canopy background adjustment that addresses the non-linear and differential transfer through a canopy of the NIR and red radiations; and C1, C2 are the coefficients of the aerosol resistance term, which uses the blue band to correct for aerosol influences in the red band.The coefficients adopted in the MODIS-EVI algorithm are; L=1, C1 = 6, C2 = 7.5, and G (gain factor) = 2.5.The EVI values range from -1 to +1, where negative values correspond to snow, ice, or water, and values closer to +1 represent a high density of green leaves (Huete et al., 2002).
We obtained EVI from the MOD13Q1.006product of the Moderate Resolution Imaging Spectroradiometer (MODIS sensor) onboard NASA's Terra satellite (Didan, 2015).MOD13Q1.006EVI product is computed from atmospherically corrected bi-directional surface reflectance by choosing automatically the best available pixel value from all the acquisitions (4 per day) in a 16-day period based on quality, cloud presence, and viewing geometry (Huete et al., 1999;Didan et al., 2015).To further remove the potential remaining effect of snow, ice, and water in our dataset, we transformed negative EVI values into zeros.Thus, we obtained a maximum-value composite image every 16 days (23 images per year).Despite its moderate spatial resolution (232 meters spatial resolution at the equator), we chose the MOD13Q1.006product as the basis for our data since it offers a long time series (almost 20 years) every 16 days, which allows for the characterization of the temporal dynamics of ecosystem functioning (Anderson et al., 2018).

MOD13Q1.006 images
were available in Google into the Earth Engine (GEE) servers https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13Q1),where we processed them.GEE combines a multi-petabyte catalog of satellite imagery and geospatial datasets with planetary-scale analysis capabilities (Gorelick et al., 2017).We used the main Javascript programming interface to build the algorithms and requests to GEE servers.EVI values are multiplied by 10,000 to store them as real numbers to occupy less disk space (both in the original MOD13Q1.006product and in our dataset).

Calculating Ecosystem Functional Attributes (EFAs)
We calculated three EFAs maps known to capture most of the variance of the seasonal curve or annual dynamics of vegetation indices (Paruelo et al., 2001;Alcaraz-Segura et al., 2006, 2009): the EVI annual mean (EVI_mean; an estimator of primary production), the EVI seasonal Standard Deviation (EVI_sSD; a descriptor of seasonality, i.e., the differences between the growing and non-growing seasons), and the date of maximum EVI (EVI_DMAX; a phenological indicator of the month with maximum EVI) (Fig. 3).To summarize the EFAs of the 2001-2018 period, we calculated the interannual mean for each attribute (using linear statistics for EVI_mean and EVI_sSD and circular statistics for EVI_DMAX).
To explore the robustness of these metrics in our study area, we examined their correlation with the first axes of a Principal Component Analysis run on the EVI annual curve of the average year (i.e., 12 EVI values corresponding to the inter-annual means of the 18 (one per year) maximum EVI values of each month.These three metrics were highly correlated with the first two PCA axes (which accumulated 96.5% of variance), as previously reported for many regions (Townshend et al., 1985;Paruelo andLauenroth, 1998, Paruelo et al., 2001;Alcaraz-Segura et al., 2006, 2009;Ivits et al., 2013).This component of the analysis is presented in full in Supplement A.

Identifying Ecosystem Functional Types (EFTs)
EFTs were identified by dividing the EFAs into four intervals and combining them into a potential number of 64 EFT classes (4 × 4 × 4), following a similar approach to Alcaraz-Segura et al., (2013) (Figure 2).For EVI_DMAX, the four intervals agreed with the four seasons of the year: January to March = Winter, April to June = Spring, July to September = Summer, October to December = Autumn.For EVI_mean and EVI_sSD, we extracted the first, second, and third quartiles for each year.We verified that an 18-year period was long enough to stabilize the quartiles by running a sensitivity analysis (see Supplement B).Then, for each quartile, we calculated the inter-annual mean of the 18-year period and used them as breaks between classes (Supplement B, Table S4).These breaks were applied back to each year as the thresholds for EVI_Mean and EVI_sSD to set EFT classes (Table 1).Finally, to summarize the EFTs, we calculated the most frequent EFT (i.e., the EFT mode for each pixel) of the 2001-2018 period.To name EFTs, we used two letters and a number: the first capital letter indicates primary production (EVI_mean), increasing from A to D; the second lower case letter represents seasonality (EVI_sSD), decreasing from a to d; finally, the numbers are a phenological indicator of the growing season (EVI_DMAX), with values 1-spring, 2-summer, 3-autumn, 4-winter (Table 1).The EFT alphanumeric code (Aa1 to Dd4) corresponds to the numeric code (1 to 64) in the .TIF files, which is shown in the legend of Figure 4d.

Deriving Ecosystem Functional Diversity metrics
We derived two diversity metrics of ecosystem functional diversity from the EFT map: EFT richness and EFT rarity.EFT richness was calculated for each year by counting the number of different EFTs in a 44pixel moving window around each pixel (top-left center pixel of a 44 kernel) (modified from Alcaraz-Segura et al., 2013).Each pixel received a richness value calculated by counting how many different EFTs there were in the 44 pixels search window.We chose a 44-pixel window because it offered the most acceptable spatial resolution without saturating the number of EFT classes per kernel (i.e., smaller kernel sizes resulted in a high proportion of moving windows saturated) (see sensitivity analysis on kernel size in Section 3.2 and Supplement C).
We calculated EFT rarity as the extension of each EFT compared to the most abundant EFT in the study area (Equation 1) (Cabello et al., 2013).Then, the average rarity map of all years was obtained.
where Area_EFTmax is the area occupied by the most abundant EFT, and Area_EFTi is the i EFT area being evaluated, with i ranging from 1 to 64.
Once we had obtanied the rarity value of each EFT (using Equation 2), we assigned to each pixel in the EFT map its rarity value.Hence, the EFT rarity map spatial resolution was the same as the resolution as the EFT map (232 m).

Characterizing inter-annual stability in ecosystem functioning
We followed two approaches to characterize inter-annual stability of ecosystem functioning (either due to inter-annual fluctuations or directional trends).First, as an estimate of inter-annual variability of EFT occurrence at pixel level, we recorded the number of different EFTs that were observed at each pixel throughout the 2001-2018 period.Second, as an estimate of inter-annual dissimilarity of EFT composition at 4x4-pixel level, we started by calculating the dissimilarity (Equation 2) in the EFT composition of each 44-pixel window (924 x 924 m) between all combinations of two years.
Then, we obtained the mean of the indices obtained from all combinations of years.We calculated dissimilarity as follows (Equation 3): Dissimilarity = 1 -Jaccard Index (Equation 3) where Jaccard similarity index is calculated as (Equations 4 and 5; Jaccard, 1901): Jaccard Index = (the number of shared classes between two sets) / (the total number of classes in either set) * 100 (Equation 4).
The same formula in notation form is (Equation 5): 5) where the Jaccard index for two data sets (X = set 1; Y =set 2) is equal to the size of the intersection divided by the size of the union of the data sets.This way, to calculate the Jaccard Index, we proceeded as follows: for each 4x4-pixel window, we first counted the number of EFTs shared between the two sets, i.e., the two years; second we counted the total number of EFTs that occurred in either set (shared and unshared between the two years); then, we divided the number of shared EFTs by the total number of EFTs; and finally, we multiplied the result by 100.Dissimilarity values ranged from 0 to 1, being 1 the highest degree of dissimilarity in EFT composition throughout years, and 0 full similarity in EFT composition throughout years.We processed inter-annual variability and inter-annual dissimilarity with IDL® software (Interactive Data Language).IDL is commonly used for interactive processing of large amounts of data, including image processing.

Effect of the EVI inter-annual variability on the boundaries set among EFT classes
To assess how inter-annual variability of EVI dynamics could affect the quartiles of EVI_Mean and EVI_sSD (which set the boundaries among EFT classes), we determined the minimum number of years needed in a study period to get stability in all quartiles (see Supplement B).For each quartile, we plotted (Figure B1) the maximum inter-annual coefficient of variation observed across all combinations of consecutive years from 2001 to 2018 (i.e., from the 17 separate combinations of two consecutive years to a single combination of all 18 years together) against the number of years considered.Starting with two consecutive years, we plotted the maximum of 17 coefficients of variation (i.e., 2001-2002, 2002-2003, … 2017-2018); for three consecutive years, the maximum of 16 coefficients of variation (i.e., 2001-2002-2003, … 2016-2017-2018), and so on.
The inter-annual Coefficient of Variation (CV) of the 2001-2018 period was around 5% for the EVI_mean quartiles and around 10% for the EVI_sSD quartiles (Table B1, Supplement B).The quartiles of EVI_Mean (our surrogate for productivity) required at least 14 years to stabilize around 5% of CV.The quartiles of EVI_sSD (our surrogate for seasonality) required at least 17 years to stabilize around 10% of CV (Figure B1, Supplement B).Despite the variation observed in the quartile values across years, we did not change the limits among EFT classes from one year to the next but instead applied the same limits to all years so we could compare the classification output across years.That is, we followed a fixed-classification approach with fixed limits among EFT classes for the entire period to make the classification able to detect such inter-annual changes.For example, if a mega wild fire burns the entire protected area in the future (e.g., in 2022), our use of fixed limits among classes for the 2001-2018 period will allow the detection of such disturbance (most pixels would be classified as low productivity "A EFT class").On the contrary, if the limits among EFT classes were adapted to the data distribution of each year, the classification would not detect the effect of wildfire on EVI dynamics and would impede the 2022 classification to be compared to previous years.

Kernel size and borderline effect on EFT richness
To assess the effect of the size of the sliding window kernel on EFT richness and rarity, we calculated EFT richness for different kernel sizes (22, 33, 44, 55 pixels) and compared the outputs (see analysis in Supplement C).The 44-pixel kernel offered the most satisfactory spatial resolution of the EFT richness map without saturation of this variable (Figure C1, Supplement C).When the size of the sliding window kernels was 22 or 33 pixels, there was a high proportion of kernels that reached the highest possible richness value for that kernel size (i.e., 4 and 9 EFT classes per kernel, respectively), so the EFT richness variable was highly saturated.The 55-pixels kernel never reached the maximum number of pixels in a kernel but resulted in too coarse outputs (grain size of 55 MO13Q1 pixels, i.e., 11501150).Hence, the 44-pixel kernel offered a balance between output resolution and variable saturation, since we observed a maximum EFT richness of 13, while the maximum potential richness in a 44-pixel kernel was 16.
Pixels with NoData values were not considered a distinct class to compute EFT richness along the study area borderline.For this reason, it is important to note that the sliding windows along the borderline of the study area could systematically show lower EFT richness in our dataset than reality.

Data structure and availability
Overall, the collection of datasets that we present here provides a characterization of ecosystem functioning, ecosystem functional diversity and inter-annual dynamics in Sierra Nevada Biosphere Reserve (SE Spain) through a key measure based on primary productivity derived from remote sensing.To broaden the use of data, we provided the datasets in two institutional scientific repositories.Datasets are available for downloading in PANGAEA: https://doi.pangaea.de/10.1594/PANGAEA.924792?format=html#download (Cazorla et al., 2020a), and we have also developed an ad-hoc visualization for the inter-annual summaries at Sierra Nevada Global Change Observatory-LTER site (http://obsnev.es/apps/efts_SN.html).To increase the transferability of the work, we also provide open source to calculate the EFTs in GEE: https://doi.org/10.5281/zenodo.7524973.
The dataset is structured in three main subsets of variables: Ecosystem Functional Attributes (EFAs), Ecosystem Functional Types (EFTs), and Ecosystem Functional Diversity (see Table 2).For each subset of variables, there are two groups of data (two subfolders): one containing the yearly variables, and another one containing the summaries for the 18-year period.A Data Management Plan with our dataset formal metadata is also available in PANGAEA data repository.

Data attribution
The MODIS database used in this work is maintained by NASA (satellite Terra, sensor MODIS, product MOD13Q1.006)and copied by Google into the Earth Engine servers (https://developers.google.com/earthengine/datasets/catalog/MODIS_006_MOD13Q1).The Sierra Nevada Biosphere Reserve boundaries shapefile is included in the public official geodatabase of the Andalusian regional government (REDIAM:https://descargasrediam.cica.es/7_PATRIMONIO_NATURAL/01_ESPACIOS_PROTEGIDOS.

Data usage in Ecology and Conservation
Ecological research based on time series of spectral vegetation indices plays an essential role in biodiversity conservation (Cabello et al., 2012;Pettorelli, 2016Pettorelli, , 2018) ) and management (Pelkey et al., 2003;Cabello et al., 2016), and for the study of biodiversity and ecosystems responses to environmental changes (Alcaraz-Segura et al., 2017;Pérez-Luque et al., 2015a, Skidmore et al., 2021).Recently, the use of EFAs derived from spectral vegetation indices in species distribution models has made possible to evaluate with high spatial and temporal precision habitat suitability for plant (Arenas-Castro et al., 2018) and animal species (Requena-Mullor et al., 2017;Regos et al., 2019) and may even anticipate expected changes in the distribution of threatened species because of climate change (Alcaraz-Segura et al., 2017).EFAs are also the basis of the monitoring program of the Spanish National Parks Network to identify changes and anomalies in ecosystem functioning, and to inform managers of ecosystem health and conservation status (Cabello et al., 2016).
Datasets based on an EFT approach can be useful for different purposes: to characterize spatial and temporal heterogeneity of ecosystem functioning from local to global scales (Alcaraz-Segura et al., 2006, Fernández et al., 2010, Cabello et al., 2013;Ivits et al., 2013); to evaluate the environmental and human controls of ecosystem functional diversity (Alcaraz-Segura et al., 2013); to identify priorities for Biodiversity Conservation (Cazorla et al., 2020b); to assess the representativeness of environmental observatory networks (Villarreal et al., 2018); to assess the effects of land-use changes on ecosystem functioning (Oki et al., 2013); and to improve weather forecast models (Lee et al., 2013;Müller et al., 2014).

Case study
Sierra Nevada is a mountain range between 860 and 3482 m a.s.l covering more than 2000 km2 in SE Spain (Fig. 1).It is one of the most critical biodiversity hotspots in the Mediterranean region (Blanca et al., 1998;Cañadas et al., 2014), hosting 105 endemic plant species and a total of 2353 taxa of vascular plants (33% and 20% of Spanish and European flora, respectively; Lorite 2016).Forest cover in Sierra Nevada is dominated by pine plantations covering approximately 40000 ha.The primary native forests are dominated by the evergreen holm oak Quercus ilex subsp.ballota (Desf.)Samp.occupying low and medium mountain areas (8800 ha) and by the deciduous Pyrenean oak Quercus pyrenaica Willd.ranging from 1100 to 2000 m a.s.l.(about 2000 ha).Above the tree line, plant communities of the Oromediterranean and Crioromediterranean belts (above 1800-2000 m a.s.l.), dominated by chamaephytes and hemicryptophytes (scrublands, grasslands, and cliff and scree communities), are the habitat to many endemisms.Sierra Nevada receives legal protection and international recognition in multiple ways: UNESCO Biosphere Reserve (1986), Natural Park (1989), National Park (1999), Important Bird Area (2003), Special Area of Conservation in Natura 2000 network ( 2012), and it is in the IUCN Green List of Protected Areas (2014), a global standard of best practice for area-based conservation.Sierra Nevada is also a site within the European Long Term Ecological Research (LTER) network, with many available ecological data records from multiple disciplines (Zamora et al., 2017, LTER_EU_ES_010), DEIMS link https://deims.org/e51cee43-dc12-4545-8e5b-dad35431e3f7.
The dataset presented here provides the first characterization of functional diversity at ecosystem level for the Sierra Nevada Biosphere Reserve.Our dataset provides a baseline of ecosystem functioning for Sierra Nevada ecosystems, which opens the possibility to assess responses of ecosystem functioning to global change and management actions, to understand the drivers of ecosystem functioning and functional diversity, and to assess the supply of ecosystem services.This dataset provides valuable information for the Global Change Observatory of Sierra Nevada, a long-term ecological research site (name: ES-SNE, code: LTER_EU_ES_010) established more than a decade ago (Zamora et al., 2017), which is also now a mountain node of the LifeWatch ERIC (European Research Infrastructure Consortium).Thus, our dataset provides information at the level of ecosystem functioning that can be combined with other available datasets on this LTER site concerning on species distributions and dynamics, climate, ecosystem services, hydrology, land-use changes, and management practices (Pérez-Luque et al., 2014, 2015b, 2015c, 2016;Ros-Candeira et al., 2019, 2020;Lorite et al., 2020).

Description of the data and insights for its use Ecosystem Functional Attributes spatial patterns
Functional attributes of productivity, seasonality, and phenology showed a clear altitudinal pattern.Productivity (EVI_mean) was much lower above the treeline, i.e. in the high mountain bioclimatic belts (Cryoro-and Oromediterranean belts) than in lower belts (Supra-and Mesomediterranean belts).Productivity also decreased from west to east (Fig. 4a).Seasonality (EVI_sSD) was the highest in the Supramediterranean belt (Fig. 4b).Phenology (EVI_DMAX) was characterized by a dominant summer peak in vegetation greenness in the Cryoro-and Oromediterranean belts, and by a late spring peak in the Supra-and Mesomediterranean belts.Dry and semi-arid Thermomediterranean areas of the south and east showed greenness peaks in early autumn and winter months (Fig. 4c).

Ecosystem Functional Types map
As a result of the combination of the three Ecosystem Functional Attributes, productivity, seasonality, and phenology (Fig. 4 a-c), we obtained the EFT map (Fig. 4d).This map represents a synthetic characterization of ecosystem functioning spatial heterogeneity based on the primary production dynamics.A total of 64 EFTs classes were observed.High-mountain areas (Cryoro-and Oromediterranean bioclimatic belts) showed EFTs with low and intermediate productivity, high seasonality, and maximum greenness in summer and spring.The extreme conditions of these areas, characterized by poorsoils (Peinado et al., 2019), high solar radiation, extreme temperatures, winds, snow, and ice, constrain so much the vegetative period that they are known as "high-altitude cold desert" (Blanca et al., 2019).Mid-mountain areas (Supra-and Mesomediterranean belts) were associated with EFTs of intermediate-high productivity, medium-low seasonality, and maximum greenness in spring and autumn (e.g., Cc1-3) (Fig. 4d).The low dry and semiarid areas (Thermomediterranean belts) of the south and east, characterized by thermophilic and xerophytic species, displayed different EFTsto the rest of the park, with very low productivity, mediumlow seasonality, and maximum greenness in spring or winter (e.g., Ac1-4).

Functional diversity at the ecosystem level
The highest EFT richness was observed in the Supra-and Mesomediterranean belts (Fig. 5c).Such ecosystem functional diversity hotspots (i.e., EFT hotspots) could be related to two facts.First, many Mediterranean mountains show high beta diversity (in terms of species turnover) around 1700-1900 m a.s.l.(Wilson and Schmida, 1984), where there is a great structural and compositional replacement of vegetation.Second, for the case of Sierra Nevada, in the mid-mountain and especially in its southern face, there is a very diverse land cover mosaic, composed by different types of natural vegetation (Valle et al., 2003), mixed with different types of pine afforestations, traditional croplands, and land-uses (Camacho et al., 2002).The areas with the lowest EFT richness were located in Oro-and Crioromediterranean belts, and in the eastern semi-arid Thermomediterranean extreme, where the harsh soil and climatic conditions (Peinado et al., 2019) diminish floristic diversity (Fernández Calzado et al., 2012).EFT rarity was highest in the peaks (above 2800 m a.s.l., Cryoromediterranean belt) and the lowest areas of the Eastern side of Sierra Nevada (semi-arid Thermomediterranean belt), both areas characterized by a high concentration of narrowly endemic species (Mota et al., 2004;Cañadas et al., 2014;Peñas et al., 2019.The high mountain areas (Oromediterranean belt) showed the lowest EFT rarity, since their EFT composition was the most abundant and extensive in the Biosphere Reserve.Mid-mountain areas (Supra-and Mesomediterranean belts) (Fig. 5d) showed medium to high EFT rarity values.The relatively higher rarity of ecosystem functioning in these belts was associated with the particular phenology of coniferous forests with autumnwinter maxima, also identified in other areas of the Iberian Peninsula (Aragones et al., 2019).

Stability in the ecosystem functioning
The inter-annual variability ranged from 1 to 17 different EFTs over the 18-years period in the same pixel (Fig. 5a).The inter-annual variability and inter-annual dissimilarity (1 -Jaccard index) (Fig. 5b) observed was higher in the Supra-and Mesomediterranean levels, coinciding with the altitudinal range where interannual climate variability is also higher (e.g., areas above the treeline are subjected to both cold-snowy years and warm-dry years, Zamora et al., 2016).The eastern end of the semi-arid Thermomediterranean areas also displayed a high inter-annual variability.There also exists significant climate fluctuation in these areas, where small changes in the seasonal pattern of precipitation can produce large changes in primary production (Houérou et al., 1988;Cabello et al., 2012b).On the other hand, the most inter-annually stable areas (i.e.those that changed the least during the period) were located in the Meso-Oromediterranean and Crioromediterranean belts, specifically, in the oak forests and high-mountain meadows, ecosystems that are subjected to low human intervention.

Conclusion
We introduce a new dataset based on the processing and analysis of the temporal dynamics of the Enhanced Vegetation Index (EVI) data from the MODIS sensor (MOD13Q1.006).The dataset contains EFAs, EFTs, EFT richness, and EFT rarity across the Sierra Nevada Biosphere Reserve (SE Spain).EFAs represent functional attributes at the ecosystem level related to primary production, seasonality, and phenology of carbon gains.EFTs are groups of pixels of the land surface that share similar dynamics in the exchanges of matter and energy between the biota and the physical environment, derived from the combination of the EFAs.EFT richness and EFT rarity are two metrics that provide information about the spatial heterogeneity of primary production as a focal ecosystem function.We also provide an estimation of the inter-annual stability of the former functional variables both at the pixel and at landscape levels throughout the 2001-2018 period.For this, we estimated the EFT inter-annual variability at each pixel and the inter-annual dissimilarity in the composition of EFTs at 4x4 sliding-window resolution.Overall, these data characterize the spatial and temporal patterns of a key measure of ecosystem functioning and ecosystem functional diversity.The EFT approach adopted here improves our understanding of ecosystem processes and ecosystem response through environmental gradients.It provides scientists and managers valuable information to identify conservation priorities, assess ecosystem response to environmental changes, estimate ecosystem services provision, and to model species distributions and ecological and hydrological processes.

FiguresFigure
Figures Figure 1.Study area: Sierra Nevada Biosphere Reserve.a) Location in the context of the Iberian Peninsula; b) remote view of Sierra Nevada mountain region (image from the International Space 715 Station took in December 2014; courtesy of "Earth Science and Remote Sensing Unit, 615 NASA Johnson Space Center"); c) delimitation of the Biosphere Reserve and the distribution of the main ecosystems (Pérez-Luque et al., 2019) and thermotype bioclimatic belts (Molero-Mesa and Marfil, 2015).

Figure 2 .
Figure 2. Workflow to characterize the ecosystem functioning and functional diversity of Sierra Nevada.MODIS (Moderate Resolution Imaging Spectroradiometer) sensor product MOD13Q1 onboard NASA's Terra satellite was used.This product contains maximum value composite images with 16-day temporal resolution (23 images per year) and 232 m spatial resolution for the Enhanced Vegetation Index (EVI) (1).The study period was from 2001 to 2018.Three ecosystem functional attributes (EFAs) describing ecosystem functioning were calculated from the EVI seasonal curve for each year (2 and 3).The range of values for each attribute was divided into four intervals, resulting in a potential number of 64 ecosystem functional types (EFTs; 444=64) (4).From EFTs, we derived four metrics related to ecosystem functional diversity (i.e., EFT richness and rarity) and ecosystem functional stability (i.e., inter-annual variability and dissimilarity) (5).

Figure 3 :
Figure 3: Seasonal dynamics of Enhanced Vegetation Index (EVI) and EVI derived metrics or Ecosystem Functional Attributes (EFAs).The "X" axis corresponds to the months of the year and the "Y" axis to the EVI values.EFAs were: the annual mean of EVI, an estimator of annual primary production (EVI_mean); the EVI seasonal coefficient of variation, i.e. a descriptor of seasonality or the differences between the growing and non growing seasons (EVI_sSD), and the date of maximum EVI, a phenological indicator of the growing season (EVI_DMAX).We chose these three EVI metrics or EFAs since they capture most of the variance (96.5%) of the EVI seasonal dynamics in a Principal Component Analysis.

Figure 4 .Figure 5 .
Figure 4. Ecosystem Functional Attributes (EFAs; a-c) and Ecosystem Functional Types (d) describing the functioning of vegetation canopy based on the Enhanced Vegetation Index (EVI), derived from MOD13Q1-TERRA (pixel 232 m) for the period 2001-2018.EFAs were (ac): the annual mean of EVI, an estimator of annual primary production (EVI_mean); the EVI seasonal coefficient of variation, i.e. a descriptor of seasonality or the differences between the growing and non growing seasons (EVI_sSD), and the date of maximum EVI, a phenological indicator of the growing season (EVI_DMAX).To name EFTs (d), we used two letters and a number: the first capital letter 745