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

20 Conservation Biology faces the challenge of safeguarding the ecological processes that sustain biodiversity. Characterization and evaluation of these processes can be carried out through attributes or functional traits related to the exchanges of matter and energy between vegetation and the atmosphere. Nowadays, the use of satellite imagery provides useful methods to produce a spatially continuous characterization of ecosystem functioning and processes at regional scales. Our dataset characterizes the patterns of ecosystem functioning in Sierra Nevada 25 (Spain) from the vegetation greenness dynamics captured through the spectral vegetation index EVI (Enhanced Vegetation Index) since 2001 to 2018 (product MOD13Q1.006 from MODIS sensor). First, we provided three Ecosystem Functional Attributes (EFAs) (i.e., descriptors of annual primary production, seasonality, and phenology of carbon gains), as well as their integration into a synthetic mapping of Ecosystem Functional Types (EFTs). Second, we provided two measures of functional diversity: EFT richness and EFT rarity. Finally, in 30 addition to the yearly maps, we calculated interannual summaries, i.e., means and inter-annual variabilities. Examples of research and management applications of these data sets are also included to highlight the value of EFAs and EFTs to improve the understanding and monitoring ecosystem processes across environmental gradients. The datasets are available in two open-source sites (PANGAEA: https://doi.pangaea.de/10.1594/PANGAEA.904575 (Cazorla et al. 2019) and 35 http://obsnev.es/apps/efts_SN.html), and bring to scientists, managers and the general public valuable information on the first characterization of the functional diversity at ecosystem level developed in a Mediterranean hotspot. Sierra Nevada represents an exceptional ecology laboratory of field conditions, where a long-term monitoring (LTER) program was established 10 years ago. The data availability on biodiversity, climate, ecosystem services, hydrology, land-use changes and management practices from Sierra Nevada, will allow to explore the 40 relationships between these other environmental data and ecosystem functional data that we provide in this work.

from the vegetation greenness dynamics captured through the spectral vegetation index EVI (Enhanced Vegetation Index) since 2001 to 2018 (product MOD13Q1.006 from MODIS sensor). First, we provided three Ecosystem Functional Attributes (EFAs) (i.e., descriptors of annual primary production, seasonality, and phenology of carbon gains), as well as their integration into a synthetic mapping of Ecosystem Functional Types (EFTs). Second, we provided two measures of functional diversity: EFT richness and EFT rarity. Finally, in 30 addition to the yearly maps, we calculated interannual summaries, i.e., means and inter-annual variabilities. Examples of research and management applications of these data sets are also included to highlight the value of EFAs and EFTs to improve the understanding and monitoring ecosystem processes across environmental gradients. functioning (Virginia and Wall, 2001;Pereira et al., 2013), since it has an important role in the carbon cycle (i.e., it is the energy input to the trophic web and therefore, the driving force behind many ecological processes). Moreover, primary production offers a comprehensive response to environmental changes, and constitutes a synthetic indicator of ecosystem health (Costanza et al., 1992;Skidmore et al., 2015).

70
To characterize ecosystem functioning through spectral vegetation indices, we can use the approach based on Ecosystem Functional Types (EFTs), defined as patches of the land surface that share similar dynamics in the exchanges of matter and energy 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: annual mean (a surrogate of annual primary production, the most essential 75 and integrative indicator of ecosystem functioning), annual standard deviation (a descriptor of seasonality or the differences between the growing and non-growing seasons), and the annual date of maximum (a phenological indicator of when in the year is the growing period centered). Since the concept appeared in 2001 (Paruelo et al., 2001), the EFT approach (or equivalent approaches) has exponentially grown to characterize functional heterogeneity from local to global scales (Alcaraz-Segura et al., 2006;Karlsen et al., 2006;Duro et al., 2007;80 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).
This article aims to illustrate how EFAs and EFTs can be used to assess the spatio-temporal heterogeneity and 85 inter-annual variability of ecosystem functioning in protected areas based on the vegetation dynamics captured through spectral vegetation indices (e.g. EVI). We introduce as a proof of concept the case of Sierra Nevada Biosphere reserve (SE Spain) from 2001 to 2018. First, for each year, we provide three Ecosystem Functional Attributes (EFAs) (i.e., annual primary production, seasonality and phenology of carbon gains), as well as their integration into a synthetic mapping of Ecosystem Functional Types (EFTs). Second, we present two measures of 90 functional diversity: EFT richness and EFT rarity. Finally, in addition to the yearly maps, we calculated interannual summaries, i.e., inter-annual means and inter-annual variability, to show the average conditions as well as the most stable and variable zones along the period (workflow in Fig. 2).

Site Description
Sierra Nevada (Andalusia, SE Spain) is a mountainous region covering more than 2,000 km 2 with an elevation range of between 860 and 3,482 m a.s.l (Fig. 1). It is considered one of the most important biodiversity hotspots

115
In Sierra Nevada, vegetation studies have mainly been developed considering a compositional perspective (phytosociological method) or successional perspective (vegetation series). These studies have been very useful for describing the vegetation heterogeneity at mesoscale (Loidi, 2017), for characterizing habitats of conservation importance (EU Directive 92/43/EEC), and for developing forest restoration policies (Valle et al., 2003).

120
However, these approaches are difficult to monitor the effects of environmental changes and management actions, to understand the environmental gradients at protected area scale that drive biodiversity patterns, and to evaluate the role of ecosystems as suppliers of benefits to society (Cabello et al., 2019).

Satellite images of Vegetation Indices (MOD13Q1 data product)
The characterization of ecosystem functioning in Sierra Nevada was based on the temporal dynamics of the Values of EVI*10,000 are given as real numbers between 0 and 10,000.

Calculating Ecosystem Functional Attributes (EFAs)
We identified three EFAs that are known to capture most of the variance in the time series of vegetation indices 135 and that are biologically meaningful (Paruelo et al., 2001;Alcaraz-Segura et al., 2006, 2009. These attributes were calculated from the EVI seasonal curve or annual dynamics. From the EVI seasonal curve of each year, we identified three functional attributes: 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 140 with maximum EVI) (Fig.3). To summarize the EFAs of the 2001-2018 period, we calculated the inter-annual mean and the inter-annual variability for each attribute.

Identifying Ecosystem Functional Types (EFTs)
EFTs were identified by synthesizing in a single map the variability contained in the three EFAs following a 145 similar approach to Alcaraz-Segura et al. (2013). 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_sSD, we extracted the first, second, and third quartiles for each year and then calculated the inter-annual mean of each quartile (means of the 18-year period) ( Table 1). These fixed limits between EFT classes were applied to each year. To summarize the EFTs of the 2001-2018 period, we calculated 150 the most frequent EFT of the period (i.e., the EFT mode for each pixel). 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_SD), 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.

Characterizing Ecosystem Functional Diversity
To characterize ecosystem functional diversity, we used EFT richness and EFT rarity. EFT richness was calculated for each year by counting the number of different EFTs in a 4×4-pixel moving window (924 x 924 m) around each pixel (top-left center pixel of the 4x4 Kernel) (modified from Alcaraz-Segura et al., 2013). Then, the average richness map of all years was obtained. EFT rarity was calculated for each year as the relative extension of each

160
EFT compared to the most abundant EFT (Equation 1) (Cabello et al., 2013). Then, the average rarity map of all years was obtained.

165
where Area_EFTmax is the area occupied by the most abundant EFT and Area_EFTi is the area of the i EFT being evaluated, with i ranging from 1 to 64.

170
To identify the most stable and variable areas (either due to inter-annual fluctuations or to directional trends) in ecosystem functioning, we provide three approaches. First, we calculated the inter-annual variability of each EFA (coefficient of variation for EVI_mean and EVI_sSD, and circular standard deviation for EVI_DMAX). Second, we recorded the number of different EFTs that occurred in the same pixel in the period 2001-2018. Third, to consider the changes not only at the pixel but also at the landscape level, the Jaccard similarity index (Jaccard, 175 1901) (Equation 2) was used in 4×4-pixel moving windows (924 x 924 m).

Jaccard Index = (the number in both sets) / (the number in either set) * 100
The same formula in notation is (Equation 3):

180
J(X,Y) = |X∩Y| / |X∪ Y| In Steps: 1) Count the number of EFTs which are shared between both windows; 2) Count the total number of EFTs in both windows (shared and unshared); 3) Divide the number of shared EFTs 1) by the total number of EFTs 2); 4) Multiply the number found in 3) by 100.

Dissimilarity = 1 -Jaccard Index
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.

210
The dataset is structured in three main subsets of variables: Ecosystem Functional Attributes, Ecosystem Functional Types, and Ecosystem Functional Diversity (see Table 2). For each variable there are two groups of data (two subfolders): one containing the yearly variables, and another one containing the summaries for the 18year period.

Ecosystem Functional Attributes patterns
Functional attributes of productivity, seasonality and phenology showed a clear altitudinal pattern. Productivity (EVI_mean) was much lower in the Crioro-and Oromediterranean bioclimatic belts than in the Supra-and Mesomediterranean belts. Productivity also decreased from west to east (Fig. 4a). Seasonality (EVI_sSD) was 220 high in the Supramediterranean, decreasing in Meso-, and Thermomediterranean belts, and in Crioro-and Oromediterranean (Fig. 4b). Phenology (EVI_DMAX) was characterized by a dominant summer peak in vegetation greenness in the Crioro-and Oromediterranean belts, and a late spring peak in the Supra-and Mesomediterranean belts. Dry and semi-arid thermomediterranean areas of the south and east showed greeness peaks in early autumn and winter months (Fig. 4c).

Ecosystem Functional Type patterns
As a result of the combination of the three functional attributes of the canopy, productivity, seasonality and phenology, represented in Fig. 4 a-c, we obtained the EFTs map ( Fig. 4d) that includes a synthetic characterization 230 of the spatial patterns of ecosystem functioning. A total of 64 classes were observed. The most abundant EFT presented the maximum greenness in spring, with productivity values from low to intermediate and with all possible combinations of seasonality: Aa1, Ba1, Cb1, Cd1, Ba1, and Cc1 accumulated 30% of the surface. On the contrary, the rarest EFTs were Ba4, Aa4 characterized by medium or low productivity, high seasonality and maximum greenness in winter.

235
Crioro and oromediterranean areas presented EFTs with low and intermediate productivity, high seasonality and moments of maximum greenness mainly in summer, but also in spring. Here, extreme conditions characterized by scarce soil (Peinado et al., 2019), high solar radiation, extreme temperatures, winds, snow and ice, give rise to a short vegetative period. This results in scarce vegetation cover, controlled by low temperatures, which can only 240 occur in summer, being the plant growth time, hence these areas have been referred to as "cold desert" (Blanca et al., 2019). The supra-and mesomediterranean levels had associated EFTs of intermediate-high productivity, medium-low seasonality and maximum green moment in spring and autumn (e.g., Cc1-3) (Fig. 4d). The supramediterranean is characterized by the presence of deciduous species, e.g., oak groves associated with the most productive and seasonal ecosystem functional type of the study area, with maximum in spring (EFT Da1).

245
In the dry and semi-arid thermomediterranean of the eastern end, characterized by thermophilic species, which hardly suffer from frost, a different functional behaviour of the ecosystems was detected. The functioning of this area showed low values of productivity, medium-low seasonality and maximum greenness of the vegetation in spring or winter (e.g., Ac1-4). Here, the main control of ecosystem functioning is water availability, presenting plant species with a fast response to scarce water inputs (Cabello et al., 2012).

Stability in ecosystem functioning
The interannual variability ranged from 1 to 17 different EFTs over the 18-year period in the same pixel (Fig. 5a). The number of EFTs observed in the same pixel over 18 years was higher in the supra-and mesomediterranean 255 levels, coinciding with the altitudinal range where interannual climate variability is most affected (e.g., they may present a lot of snow in cold years and be affected by drought in dry and warm years). The eastern end of the semi-arid thermomediterranean also highlighted with high inter-annual variability, where there exists a greater climate fluctuation and where small changes in precipitation produce large changes in the dynamics of primary production (Houérou et al., 1988;Cabello et al., 2012), as well as the area burned in 2005 near Lanjarón, where 260 the fire eliminated the vegetation that has been regenerating since then. On the other hand, the most stable vegetation types interannual, i.e., those that changed the least during the period, were located in the mesooromediterranean and crioromediterranean levels, specifically, the oak and borreguil vegetation types, ecosystems that are not subject to anthropic presence (e.g., low forest management and low presence of livestock).

265
The results of the inverse of the Jaccard coefficient to obtain the dissimilarity or functional changes between years in the composition of EFTs over the 2001-2018 period (Fig. 5b), showed an altitudinal pattern where the dissimilarity between EFTs was lower in the oro and cryoromediterranean levels, as well as in the mesomediterranean oak groves (functional stability already shown by other authors, i. e. Requena-Mullor et al, 2018). This pattern of dissimilarity increased towards lower levels, finding the highest values of dissimilarity (or 270 greater change) in areas where changes in land use and management are major (Zamora et al., 2016), such as autochthonous pine forests on dolomites, coniferous repopulations and meso-and thermomediterranean holm oaks. In addition, the eastern end of the Sierra Nevada had an area with low dissimilarity values, that is, where there were not many changes over the years and when they occurred they were towards very similar EFTs.

Functional diversity at ecosystem level
Richness oscillated between 1 and 13 EFTs. Highest EFT richness was observed in the supra-and mesomediterranean, particularly in the southern face of the Sierra (Fig. 5c), where the number of vegetation series is also greater than in the rest of the bioclimatic floors (Valle et al., 2003). The presence EFTs hotspots mainly in 280 the mid-mountain, and particularly in the southern face, could be related to two factors. On the one hand, many Mediterranean mountains show high values of beta diversity up to 1750-1800 m (Wilson and Schmida, 1984 (Mota et al., 2004). In the oromediterranean, EFT rarity decreased and reached its minimum, due to the great extension in the Sierra Nevada 305 of this bioclimatic level, which made its functioning not appear as rare, and increasing again in the supra-and mesomediterranean (Fig. 5d). The most rare supra-and mesomediterranean vegetation types corresponded to coniferous and holm oak repopulations (rarity 0.6). The high rarity of coniferous repopulations may be due to disturbances or management interventions that give rise to unique functions in the different masses of conifers. On the other hand, the rarity in holm oaks may be due to their exclusive functioning, i.e. they have very specific 310 associated EFTs (e.g., Cc1, Dc1). However, the rarity of the different vegetation types (between 0.45 and 0.64) was far from the maximum possible (1). (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. 2015). In fact, numerous studies have demonstrated the usefulness of satellite image time series to evaluate the functional changes in ecosystems at regional scale (Alcaraz-Segura et al., 2010) and at the protected area level (Alcaraz-Segura et al., 2009;Lourenço et al., 2018). Recently, the use of EFAs derived from spectral 320 indices of vegetation in species distribution models, has made it possible to evaluate with great spatial and temporal precision the suitability of habitat for plant species (Arenas-Castro et al., 2018) and animals (Requena-Mullor et al., 2017;Regos et al., 2019) and may even anticipate expected changes in the distribution of plant species threatened as a consequence of climate change (Alcaraz-Segura et al., 2017). In addition, based on the EFAs, a monitoring programme of the Spanish National Parks Network has been designed to identify changes 325 and anomalies in functioning, informing managers of the health and conservation status of ecosystems (Cabello et al., 2016). Furthermore, EFTs have been used to characterize spatial and temporal heterogeneity of ecosystem functioning at local and regional scales (Fernández et al., 2010;Cabello et al., 2013); to describe biogeographical patterns 330 (Alcaraz-Segura et al., 2006;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., 2019); to assess the representativeness environmental networks (Villarreal et al., 2018); to assess the effects of land-use changes on ecosystem functioning (Oki et al., 2013); or to improve weather forecast models (Lee et al., 2013;Müller et al., 2014).

335
The data sets that we are providing give to the scientific community valuable information of the first characterization of the functional diversity at ecosystem level developed in the entire protected area. We provided a detailed characterization of the functional diversity at ecosystem level for Sierra Nevada, that could be useful to monitor the response of ecosystems to global change and management actions, to understand the ecosystem 340 functioning and functional diversity across the environmental gradients at protected area scale, and to evaluate the role of ecosystems in providing ecosystem services (Cabello et al., 2019). Sierra Nevada is also a long-term ecological laboratory established 10 years ago (Zamora et al. 2016(Zamora et al. , 2017, that have available data on biodiversity, climate, ecosystem services, hydrology, land-use changes and management practices from Sierra Nevada. This will allow to explore the relationships between these other environmental data with the ecosystem functional data 345 that we provide.

350
The datasets described in this article are available in open-access sources. To broaden their use, first, we provide data in .tif format. Second, we have incorporated rendered versions of all layers as required by Google Earth Pro (called "filename..._forGoogleEarthVisualization.tif") for visualization. And third, we have also developed an adhoc visualization platform for all the layers. Datasets available for download in PANGAEA: The MODIS database used in this work are maintained by NASA (satellite Terra, sensor MODIS, product MOD13Q1.006), and the geospatial datasets of Sierra Nevada Park are included in public database of the Andalusian regional government (REDIAM).

6 Conclusion
This dataset provides a characterization of ecosystem functioning and ecosystem functional diversity in Sierra

365
Nevada Biosphere Reserve (SE Spain) through the analysis of time series of satellite images of spectral vegetation indices as surrogates of the carbon gains dynamics. First, three Ecosystem Functional Attributes (EFAs) describe the spatial and inter-annual variability in productivity, seasonality and phenology of vegetation photosynthetic activity. Second, the combination of these EFAs into a synthetic classification, i.e. Ecosystem Functional Types (EFTs), integrates in a single map the spatial heterogeneity of these descriptors of the seasonal dynamics of carbon 370 gains. Finally, by using EFTs as biological entities, the spatial patterns of ecosystem functional diversity were assessed by means of EFT richness and EFT rarity, as well as the inter-annual variability in ecosystem functioning through EFT inter-annual variability and EFT inter-annual dissimilarity.
Ecosystem Functional Types approach improve the understanding of ecosystem processes through environmental 375 gradients and provide both the scientific community with valuable information of the first characterization of the functional diversity at ecosystem level developed in the entire protected area.

380
DAS, JC, JP and BPC designed the study, and DAS, JC, JP coordinated it. BPC, AR and EG processed data and produced the associated data sets presented in this paper. BPC prepared the manuscript with contributions from all authors. BPC and EG prepared the final figures. AJPL design and made the application to visualize the data.
Competing interests. The authors declare that they have no conflict of interest.