Model-Assisted Bird Monitoring Based on Remotely Sensed Ecosystem Functioning and Atlas Data

remote sensing Article Model-Assisted Bird Monitoring Based on Remotely Sensed Ecosystem Functioning and Atlas Data Adrián Regos 1,2,* , Pablo Gómez-Rodríguez 1, Salvador Arenas-Castro 2,3 , Luis Tapia 1, María Vidal 1and Jesús Domínguez 1 1Departamento de Zooloxía, Xenética e Antropoloxía Física, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain; [email protected] (P.G.-R.); [email protected] (L.T.); [email protected] (M.V.); [email protected] (J.D.) 2CIBIO-InBIO—Centro de Investigação em Biodiversidade e Recursos Genéticos, Laboratório Associado, Universidade do Porto, Campus Agrário Vairão, 4485-661 Vila do Conde, Portugal; salvadorarenascastr[email protected] 3Instituto Politécnico de Viana do Castelo, Escola Superior de Tecnologia e Gestão, 4900-347 Viana do Castelo, Portugal *Correspondence: adrian.r[email protected] Received: 22 June 2020; Accepted: 5 August 2020; Published: 7 August 2020   Abstract: Urgent action needs to be taken to halt global biodiversity crisis. To be effective in the implementation of such action, managers and policy-makers need updated information on the status and trends of biodiversity. Here, we test the ability of remotely sensed ecosystem functioning attributes (EFAs) to predict the distribution of 73 bird species with different life-history traits. We run ensemble species distribution models (SDMs) trained with bird atlas data and 12 EFAs describing different dimensions of carbon cycle and surface energy balance. Our ensemble SDMs—exclusively based on EFAs—hold a high predictive capacity across 71 target species (up to 0.94 and 0.79 of Area Under the ROC curve and true skill statistic (TSS)). Our results showed the life-history traits did not significantly affect SDM performance. Overall, minimum Enhanced Vegetation Index (EVI) and maximum Albedo values (descriptors of primary productivity and energy balance) were the most important predictors across our bird community. Our approach leverages the existing atlas data and provides an alternative method to monitor inter-annual bird habitat dynamics from space in the absence of long-term biodiversity monitoring schemes. This study illustrates the great potential that satellite remote sensing can contribute to the Aichi Biodiversity Targets and to the Essential Biodiversity Variables framework (EBV class “Species distribution”). Keywords: IUCN conservation status; land-surface temperature; life-history traits; remotely sensed essential biodiversity variables (RS-EBVs); surface energy balance and temperature; vegetation productivity 1. Introduction Urgent action needs to be taken to achieve the Aichi Biodiversity Targets of the Convention on Biological Diversity and thus halt the global biodiversity crisis [ 1 ]. The tenth meeting of the Conference of the Parties calls upon countries to implement the 2011–2020 Strategic Plan for Biodiversity, including the Aichi Biodiversity Targets, with attention given to the monitoring and review of national biodiversity strategies and action plans implementation in accordance with the Strategic Plan and national targets, making use of the set of indicators developed for the Strategic Plan as a flexible framework. To be effective in the implementation of the national biodiversity strategies and action plans, managers and policy-makers need to be supported by updated information on the status and trends of biodiversity. Remote Sens. 2020,12, 2549; doi:10.3390/rs12162549 www.mdpi.com/journal/remotesensing Remote Sens. 2020,12, 2549 2 of 14 In this context, the Group on Earth Observations Biodiversity Observation Network (GEO BON) has introduced the framework of Essential Biodiversity Variables (EBVs) [ 2 ]. The EBVs—defined as “derived measurements required to study, report, and manage biodiversity change”—aims to coordinate global monitoring initiatives and decision making [ 2 ]. Species abundance and distribution are intuitive biodiversity indicators widely used to quantify population trends and extinction risks, range shifts in response to climate or land-use change; among others [ 3 ]. In the EBV framework, “Species Distribution” has been proposed as an EBV candidate within the EBV class “Species Populations” [4]. Species distribution can be represented by occurrence (i.e., presence–absence) data based on geographical observations gathered from specific surveys (e.g., atlas or monitoring programs) or predicted from species distribution models (hereafter: SDMs) [ 3 , 5 , 6 ]. However, biodiversity monitoring programs based on repeated sampling designs are scarce. In the absence of such programs, model-assisted monitoring emerges as a cost-effective alternative to provide updated information on biodiversity status and trends [ 7 ]. SDMs trained with in situ data from atlas and inventories (i.e., presence–absence data) are a good compromise between presence-only models based on opportunistic data from citizen science initiatives (often associated with a large uncertainty in data quality) and long-term monitoring schemes (based on target-oriented protocols, supported by very expensive and time-consuming sampling designs) [3,4]. EBVs, by definition, should be able to be measured or modeled in standardized and cost-effective way from local to global scales. In addition, EBVs should ideally capitalize from integrating remote sensing with in situ observations [ 8 ]. Earth Observation (EO) and Remote Sensing (RS) technologies have a great potential to contribute to the Aichi Biodiversity Targets and to EBVs framework [ 9 ]. In fact, 14 of 22 candidate EBVs can be considered Remote Sensing Essential Biodiversity Variables (RS-EBVs), as they were found to have a fully or partly remotely sensed component [ 9 ]. Moreover, recent studies have demonstrated the great potential of EO data and RS technologies to improve SDMs [ 10 ]. The lack of spatial data on ecologically relevant factors driving species distribution constraints the full potential of SDMs for management and monitoring purposes [ 11 – 13 ]. In particular, RS can provide critical information on species habitat requirements such as soil moisture (measured through e.g., normalized difference of water indices [ 14 ]), microclimatic conditions (e.g., land-surface temperature [ 15 ]) or vegetation productivity (e.g., spectral vegetation indices [ 16 ]), thus improving species niche characterization and SDMs predictions [ 11 ]. In addition, the high spatiotemporal resolution of EO data obtained from satellite sensors such as Moderate Resolution Imaging Spectroradiometer (MODIS) (from TERRA/AQUA satellites) and MultiSpectral Instrument (MSI) (from Sentinel-2 mission) allows derivation of different descriptors of ecosystem functioning (also known as remotely sensed ecosystem functioning attributes; hereafter RS-EFAs) [ 17 , 18 ]. RS-EFAs have been recently tested as predictor variables in SDMs for both plants and animals [ 19 – 23 ] (being also considered to be EBV candidates for the “Ecosystem Function” class [ 2 ]). Some of these studies have already shown how RS-EFAs related to carbon cycle and water balance can be integrated into SDM frameworks to support monitoring for rare plant species [ 19 , 20 ] or waterbird conservation plans [ 14 ]. SDMs based on RS-EFAs could therefore offer an alternative and cost-effective approach for tracking biodiversity when long-term biodiversity monitoring schemes are not available. However, the ability of SDMs to predict the distribution of species could depend on different life-history traits of the modeled species —as already found in previous research (see e.g., [ 11 , 23 ]). In fact, the effect of species traits on the predictive ability of SDMs can even override any differences in modeling technique, as these traits may reflect the different responses of species to processes that control species distributions [ 24 ]. Recent studies highlighted that the predictive ability of SDMs differs regarding life-history characteristics such as range, migration, habitat and rarity of a species [ 25 , 26 ]. Despite the increasing acknowledgement on the role of species traits on SDM performance, the effects of life-history traits on SDMs based on a wide range of EFAs remains poorly studied (but see [11,23]). Protected areas are one of the main conservation strategies worldwide to cope with the ongoing biodiversity decline. In Europe, the Natura 2000 network is the cornerstone of conservation strategies, Remote Sens. 2020,12, 2549 3 of 14 comprising more than 25,000 sites representing around 18% of the area of the 27 EU Member States [ 27 ]. In the absence of long-term biodiversity monitoring schemes, model-assisted monitoring turns out to be an essential tool for managers of Natura 2000 sites. In this study, we specifically aims to: 1) assess the ability of SDMs based on EFAs to predict the distribution of 73 bird species with different life-history traits (habitat preference, biogeographic origin, migratory status, feeding and nesting strategies), 2) explore the relative importance of the different EFAs as predictor variables of the bird distribution, and the role of life-history traits on model performance and habitat suitability changes predicted by the SDMs. Finally, we also aim to 3) propose an alternative method to monitor inter-annual bird habitat dynamics from space in the absence of bird monitoring programs. To do so, we focus on a bird community breeding on a mountain range included in the Natura 2000 network—area located between the Mediterranean and Eurosiberian biogeographic regions and largely affected by wildfires, intensive forestry practices, and land-use changes. We expect that SDMs based on different components of ecosystem functioning and local atlas data could accurately predict the distribution of a large set of species, thus providing inter-annual predictions of habitat suitability dynamics from a time period of 17 years since the Nature 2000 network implementation. We also expect to find trait-predictor relationships since different descriptors of ecosystem functioning (e.g., variables related to energy balance) may be a priori important predictors for birds with specific life-history traits (e.g., birds nesting on rocky cavities). 2. Materials and Methods 2.1. Bird Data and Study Area We tested our approach with a bird community breeding in a mountain system of Galicia, northwestern Spain (Figure 1). The study area is a mountain range included in the Natura 2000 network: Special Areas of Conservation “Macizo Central” (c.a. 46,983 ha) and “Bidueiral de Montederramo” (c.a. 1984 ha). It represents a transitional area between the Mediterranean and Eurosiberian biogeographic regions [ 28 ]. Annual temperature is 10 ◦ C, with marked variations along the altitudinal (from − 12 ◦ C to 31 ◦ C) and latitudinal gradient [ 28 ]. Annual rainfall range between 300 and 1780 mm, and snowfall is frequent in the mountaintops [ 28 ]. Landscape is dominated by shrublands (63%) and forests (31%), whilefarmlandsrepresentless than5%.The study areahasbeen largelyaffected bywildfires (seeFigureS1), intensive forestry practices (clearcutting) and land-use changes (agricultural abandonment in the mountains and intensification in the lowlands) over recent decades, which have increased open shrublands and decreased closed shrublands and forest areas (Figure S2 and [29]). Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 [27]. In the absence of long-term biodiversity monitoring schemes, model-assisted monitoring turns out to be an essential tool for managers of Natura 2000 sites. In this study, we specifically aims to: 1) assess the ability of SDMs based on EFAs to predict the distribution of 73 bird species with different life-history traits (habitat preference, biogeographic origin, migratory status, feeding and nesting strategies), 2) explore the relative importance of the different EFAs as predictor variables of the bird distribution, and the role of life-history traits on model performance and habitat suitability changes predicted by the SDMs. Finally, we also aim to 3) propose an alternative method to monitor interannual bird habitat dynamics from space in the absence of bird monitoring programs. To do so, we focus on a bird community breeding on a mountain range included in the Natura 2000 network— area located between the Mediterranean and Eurosiberian biogeographic regions and largely affected by wildfires, intensive forestry practices, and land-use changes. We expect that SDMs based on different components of ecosystem functioning and local atlas data could accurately predict the distribution of a large set of species, thus providing inter-annual predictions of habitat suitability dynamics from a time period of 17 years since the Nature 2000 network implementation. We also expect to find trait-predictor relationships since different descriptors of ecosystem functioning (e.g., variables related to energy balance) may be a priori important predictors for birds with specific lifehistory traits (e.g., birds nesting on rocky cavities). 2. Materials and Methods 2.1. Bird Data and Study Area We tested our approach with a bird community breeding in a mountain system of Galicia, northwestern Spain (Figure 1). The study area is a mountain range included in the Natura 2000 network: Special Areas of Conservation “Macizo Central” (c.a. 46,983 ha) and “Bidueiral de Montederramo” (c.a. 1984 ha). It represents a transitional area between the Mediterranean and Eurosiberian biogeographic regions [28]. Annual temperature is 10◦C, with marked variations along the altitudinal (from −12 °C to 31 °C) and latitudinal gradient [28]. Annual rainfall range between 300 and 1780 mm, and snowfall is frequent in the mountaintops [28]. Landscape is dominated by shrublands (63%) and forests (31%), while farmlands represent less than 5%. The study area has been largely affected by wildfires (see Figure S1), intensive forestry practices (clearcutting) and land-use changes (agricultural abandonment in the mountains and intensification in the lowlands) over recent decades, which have increased open shrublands and decreased closed shrublands and forest areas (Figure S2 and [29]). Figure 1. Location of the study area in Europe. White and grey lines show the limit of the Natura 2000 network and the 2 km Universal Transverse Mercator (UTM) squares, respectively. Hillshade map shows the topographical shape of hills and mountains. Elevation is expressed in meters. Figure 1. Location of the study area in Europe. White and grey lines show the limit of the Natura 2000 network and the 2 km Universal Transverse Mercator (UTM) squares, respectively. Hillshade map shows the topographical shape of hills and mountains. Elevation is expressed in meters. Remote Sens. 2020,12, 2549 4 of 14 Data on bird species were obtained from breeding bird surveys conducted in the study area between 2002–2004 and 2008 [ 30 ]. The study area was surveyed using the Universal Transverse Mercator (UTM) grid system at 2 km resolution. All surveys were carried out during the breeding season (between Abril and July). Bird observations were georeferenced by a global positioning system (GPS) to the corresponding UTM grid. The final database comprises a total of 108 species, being represented by 37 taxonomic families and 13 orders [ 30 ]. More than 50% of the recorded species have a restricted distributional range within the Nature 2000 site (i.e., species present in less than 25% of the study area; see species maps in Figure S3). In addition, 20 species are included in Annex I of the Bird Directive and in the Galician Catalog of Threatened Species, which highlights the importance of the study area for bird conservation in Galicia [ 30 ]. The database was complemented with information on life-history traits of species (habitat preference, biogeographic origin, migratory status, feeding and nesting strategies) (see Table S1; species traits extracted from the SEO/BirdLife online bird encyclopedia, https://www.seo.org/listado-aves-2/). 2.2. Satellite-Derived Products and EFAs Three MODIS satellite products were selected to describe three dimensions of ecosystem functioning: the Enhanced Vegetation Index as a surrogate for the carbon cycle dynamics; the Land-Surface Temperature as a surrogate of sensible heat dynamics; and Albedo as a surrogate for the radiative balance [ 31 , 32 ]. All MODIS products were retrieved from the MODIS Collection V006 that includes improved aerosol retrieval and correction algorithms, new aerosol retrieval Look Up Tables updates, and refined internal snow, cloud, and cloud shadow detection algorithms—which represents a major improvement for environmental modeling and monitoring applications from previous versions. The three MODIS selected products were: (1) Enhanced Vegetation Index (EVI; MOD13Q1.006). EVI is an indicator of carbon gains since it is known to be more reliable in both low and high vegetation cover situations than Normalized Difference Vegetation Index (NDVI), and resistant to both soil influences, canopy background signals, and atmospheric effects on vegetation index values [ 33 ]. EVI values ranged from − 1 to 1, with healthy vegetation generally holding values between 0.20 and 0.80. (2) Land-Surface Temperature (LST; MOD11A2.006). LST is a good indicator of the energy balance at the Earth’s surface, and one of the key parameters in the physics of land-surface processes from regional to global scales. In addition, LST is directly linked to the primary environmental regimes and to habitat suitability attributes (e.g., productivity, vegetation structure, land-cover type; [ 34 ]). Temperatures (LST) ranged from −25 ◦C to 45 ◦C. (3) Albedo (ALB; MCD43A3.006). ALB is surrogate for surface properties such as the extent and nature of the vegetation cover, and it is affected by the change of the land-surface bio-physical factors such as vegetation, LST and soil moisture [ 35 ]. ALB values ranged from 0 to 1 (fresh snow and bare soil usually fall around 0.9). MODIS predictors were computed for the 2001–2017 time period at native resolution (~250 m) and then aggregated (mean of pixels) at the sampled square level (i.e., 2 km resolution) by using the cloud-based Google Earth Engine computational platform (https://earthengine.google.com) [ 36 ]. For each of these three dimensions of ecosystem functioning (EVI, LST and ALB), we derived the following 7 summary metrics (EFAs) of their seasonal dynamics: annual mean (surrogate of annual total amount); annual maximum and minimum (indicators of the annual extremes); annual standard deviation and coefficient of variation (descriptor of seasonality); and dates of maximum and minimum (indicators of phenology) (see description in Table 1). The complete EFAs dataset included 21 ecosystem functioning attributes (7 metrics ×3 dimensions) as competing candidate predictors. Remote Sens. 2020,12, 2549 5 of 14 Table 1. Ecosystem functioning attributes derived for each Moderate Resolution Imaging Spectroradiometer (MODIS) product. EFA Name Acronyms Description Annual mean Mean Surrogate of annual total amount. Maximum Max Annual maximum (indicator of the annual extremes). Minimum Min Annual minimum (indicator of the annual extremes). Annual coefficient of variation sCV A descriptor of the differences between seasons. Annual standard deviation sSD A descriptor of the variations between seasons. Date of maximum DMx Phenological indicator of the maximum growing season. Date of minimum DMi Phenological indicator of the minimum growing season. 2.3. Species Distribution Models Bird occurrence data consisted of all presence/absence records collected during the specific surveys between 2002–2004 and 2008 within each 2 km UTM grid. From the initial pool of 108 species, we only model those species with more than 20 presences (total of 73 species) to avoid the risk of overfitting the models with a low number of occurrences [ 37 – 39 ] (see the number of occurrence available for each modeled species in Table S1). For model calibration, EFA covariates were averaged for the years 2002–2004 and 2008 to match with the years in which bird data was collected. To avoid including highly correlated variables in model fitting, we conducted a multicollinearity analysis by testing Pearson pair-wise correlations and Variance Inflation Factors (VIF). The VIF analysis estimates how much the variance of a regression coefficient is inflated due to multicollinearity in the model. In general, a VIF above 10 indicates high correlation and is cause for concern. To be more conservative, we only selected variables with VIF values lower than 3. Based on these multicollinearity analyses, we finally retained 12 independent predictors for each SDM (i.e., variables with Pearson’s correlation coefficients <0.7 or >−0.7 and VIF <3; see Figure S4 for pair-wise correlations among predictors). Species distribution models were developed with the R-package “Biomod2” [ 40 ], where seven widely used modeling algorithms were applied: generalized linear models (GLM), generalized additive models (GAM), generalized boosted models (GBM, also known as boosted regression), flexible discriminant analysis (FDA), multivariate adaptive regression splines (MARS), Breiman and Cutler’s random forest for classification and regression (RF), and artificial neural networks (ANN) [ 40 ]. An ensemble approach of different modeling algorithms was used to deal with the inherent uncertainty of individual models and thus provide more informative and ecologically correct predictions [ 41 ]. SDMs were trained and evaluated by splitting data into calibration and validation subsets, including 70% and 30% of the data, respectively. We randomly repeated this procedure 30 times to yield predictions independent of the training data [ 42 ]. The predictive performance of these models was evaluated using the area under receiver operating characteristic curve (AUC) [ 42 ] and the true skill statistic (TSS) [ 43 ]. The final consensus predictions (i.e., ensemble models) were performed by using a weighted mean approach. This consensus method considers the weights proportional to the selected evaluation scores (i.e., the higher the AUC of the model, the greater the importance in the ensemble modeling [ 38 ]). Variables importance calculation was based on decreasing in accuracy. The value returned is based on the correlations of the model prediction with formal explanatory dataset and prediction with the same dataset where the columns of variable of interest has been shuffled. Greater is the correlation, less importance the variable has. We made the sum of variable importance across models and replicates to consider relative importance of different groups of variables. The effect of life-history traits and IUCN conservation status on model performance was tested with GLMs with a gamma distribution. AUC and TSS values were used as dependent variable and the different species traits as explanatory variables. No interactions were included in the model, as the objective was not to obtain the best model but to explore the effect of each single factor on model performance. Effects were considered significant at p-values <0.01. Remote Sens. 2020,12, 2549 6 of 14 2.4. Inter-Annual Predictions of Bird Communities The ensemble SDMs were projected to the environmental conditions of each year (from 2001 to 2017) to obtain a time series of annual habitat suitability predictions. Inter-annual variation in habitat suitability predicted by the SDMs for each bird species was estimated from the difference between the habitat suitability of each year in relation to year 2001. Predicted habitat suitability values at species level were then averaged at bird guild level (i.e., groups defined according to the life-history traits of each species, see Table S1) and represented with R-package “ggplot2” [ 44 ]. We also calculated the inter-annual mean and standard deviation (SD) of predicted habitat suitability for the entire study area and 2001–2017 time period as spatiotemporal descriptors of the inter-annual habitat dynamics. 3. Results The final set of RS-EFAs yielded ensemble SDMs with high predictive capacity (average AUC and TSS values of 0.94 ± 0.05 and 0.79 ± 0.12, respectively; see Table S1 for single-species SDMs). From the 73 modeled species only two of them showed AUC values lower than 0.7, the Black redstart (Phoenicurus ochruros) and the common buzzard (Buteo buteo). No significative effects of species traits on SDM performance were found (p-values >0.05). However, SDM performance was slightly better for Mediterranean and migratory birds (AUC and TSS values up to 0.96 and 0.84, respectively; Figure 2) than for Eurosiberian and sedentary species (AUC and TSS values up to 0.94 and 0.78, respectively; Figure 2). Regarding the remaining traits, SDMs for mountain birds nesting in rocky cavities with carnivorous diets showed the highest predictive accuracy (AUC and TSS values up to 0.96 and 0.86, respectively; Figure 2). The distribution was better predicted for vulnerable species (AUC and TSS values up to 0.97 and 0.86) than for species with the other protection status (Figure 2). Figure 2. Predictive accuracy of the ensemble species distribution models measured through the area under receiver operating characteristic curve (AUC) of the ROC curve and TSS metrics. Results grouped by bird guilds defined according to life-history traits and conservation status. Box plots show the median and interquartile range along with outliers, and violin plots show the data density (i.e., number of species) along the y-axis. IUCN conservation status: Vulnerable (VU), Near Threatened (NT), Least Concern (LC) and Data Deficient (DD). Remote Sens. 2020,12, 2549 7 of 14 Overall, minimum EVI and maximum Albedo values (descriptors of primary productivity and energy balance, respectively) were the most important predictors across the bird community (Figure 3). However, the SD of EVI values (descriptors of seasonality in primary productivity), mean and SD of LST and minimum of Albedo (descriptors of sensible heat dynamics and radiance balance, respectively) were also found to be important predictors of bird distribution in our study area (Figure 3). Minimum EVI is the most important predictor across species traits, with the exception of mountain birds nesting on rocky cavities whose distribution was found to be mainly predicted by descriptors of energy balance (minimum Albedo and mean LST; see nidification strategy and habitat preference groups in Figure 3). Figure 3. Variable importance for each bird guild defined according to life-history traits. Variable importance at guild level is a cumulative sum across the different species, modeling algorithms, and replicates. Ensemble SDMs predicted a large inter-annual variation in habitat suitability between 2001 and 2017 (Figures 4and 5). Overall, our models predicted sharp declines in habitat suitability for the bird community in years 2008 and 2014, with marked increases in years 2004–05, 2010–11 and 2017 (Figure 5). Despite the large variability across species, the predicted habitat dynamics were found to depend on the life-history traits and conservation status of the species (Figure 5). On average, migratory species ( − 0.62%) with preference for mountain habitats ( − 5.41%) nesting in rocky areas ( − 7.09%), Mediterranean climate conditions ( − 4.88%), carnivorous diets ( − 1.38%) and “Near Threatened” IUCN status ( − 8.64%) were predicted to experience an overall decline in habitat suitability for the time period 2001–2017 (Figure 4). Remote Sens. 2020,12, 2549 8 of 14 Figure 4. The spatiotemporal descriptors of the inter-annual habitat dynamics calculated from the predicted habitat suitability for the entire study area and 2001–2017 time period. The y-axis represents: (1) the inter-annual average of habitat suitability change over the 17-year period (HS_ch_mean; as indicator of total habitat suitability change), and (2) the inter-annual standard deviation (HS_ch_sd; as indicator of their inter-annual variability; the larger the standard deviation, the larger the inter-annual variability in habitat suitability). Results are grouped by bird guilds defined according to life-history traits and conservation status. Box plots show the median and interquartile range along with outliers, and violin plots show the data density (i.e., number of species) along the y-axis. IUCN conservation status: Vulnerable (VU), Near Threatened (NT), Least Concern (LC) and Data Deficient (DD). Figure 5. Inter-annual habitat dynamics between 2001 and 2017 grouped by bird guilds defined according to life-history traits and conservation status. The y-axis represents the habitat suitability change predicted from SDMs over the 17-year period (HS_ch), estimated from the difference between the predicted habitat suitability of each year in relation to year 2001. For all plots, colored lines indicate mean values while the transparent colored areas indicate the confidence limits for the population mean without assuming normality calculated through a nonparametric bootstrap. Remote Sens. 2020,12, 2549 9 of 14 4. Discussion The results showed that remotely sensed descriptors of ecosystem functioning are good predictors of bird distribution in our region (see Figure 2). Our ensemble species distribution models, exclusively based on EFAs, hold a high predictive capacity across 71 of our 73 target species (up to 0.94 and 0.79 of AUC and TSS, respectively; see Figure 2). This high model performance further supports to use of remotely sensed EFAs in species distribution models, as recently demonstrated for other taxa and spatiotemporal scales [ 19 , 20 , 23 ]. However, our SDMs were only validated within the calibration time frame (i.e., defined by bird atlas data) (see Figure 2). Although previous studies suggest that the higher model performance, the higher the transferability, a high model performance within the calibration time frame do not necessarily guarantee the ability of the model to predict changes in time [ 11 , 45 ]. Given the uncertainty when extrapolating the model beyond the range of data used for its calibration, further research will be needed to test the ability of SDMs based on RS-EFAs to predict species distribution in time [ 46 – 48 ], and not only in space (by, e.g., evaluating model projections over time with time series of bird distribution data). In addition, although the spatial heterogeneity is partially included in the SDMs through the spatial variation of the EFA across space (i.e., across the environmental gradient of the study area), some bird species could also benefit from the inclusion of other EFAs informing about the spatial heterogeneity within—and not only between—the sampled squares [ 49 – 51 ]. Our results showed that bird traits did not significantly affect SDM performance, even though the distribution was slightly better predicted for migratory and Mediterranean birds nesting in rocky areas, with carnivorous diets and vulnerable IUCN status (Figure 3). A recent study using remotely sensing descriptors of vegetation productivity also found that migrant birds were better predicted than sedentary species [ 11 ]. These results are supported by recent studies that demonstrated that migrant birds track intra-annual vegetation dynamics [ 52 ], i.e., the seasonally progressing green-up of vegetation [ 53 , 54 ]. Other organisms besides birds, such as rare plants or mammals, were also found to be well predicted by satellite descriptors of vegetation productivity [ 19 – 22 ]. For instance, a recent study found that the distribution of the European badger (Meles meles) was better predicted after including remotely sensed indicators of primary productivity (measured also from the Enhanced Vegetation Index)—which was related to food resources and habitat quality [ 21 ]. The overall high accuracy of our SDMs—which led to low variances in AUC and TSS metrics across species and, in turn, to no significative effect of species traits on model performance—suggests that our EFAs dataset was able to characterize well the niche of species with different life-history traits. This might be explained by the inclusion of not only EFAs related to carbon cycle (i.e., computed from time series of vegetation indices) but also other attributes related to heat dynamics and energy balance (i.e., computed from LST and Albedo) (see Figure 3). Our results underline the utility of remotely sensed descriptors of ecosystem functioning to better characterize the ecological niche of species (see Figure 3). For instance, forest species such as Eurasian nuthatch (Sitta europaea) or the European turtle dove (Streptopelia turtur) were mainly associated with vegetation productivity (minimum EVI values) (see Figure S5 for variable importance at species level)—which has been related to shelter and nesting habitats [ 55 , 56 ]. Birds nesting in rocky cavities such as the common raven (Corvus corax), the common rock thrush (Monticola saxatilis) or the Eurasian crag martin (Ptyonoprogne rupestris) were correlated with sensible heat dynamics and radiance balance (measured from minimum or maximum Albedo and mean or SD of LST) (see Figure S5). The intra-annual mean of LST time series was also critical for the niche characterization of Mediterranean species such as the Subalpine warbler (Sylvia inornata), the common nightingale (Luscinia megarhynchos) or the Melodious warbler (Hippolais polyglotta) (see Figure S5). The inclusion of remotely sensed indicators of energy balance allows us to accurately predict the distribution of Mediterranean, mountain birds nesting in rocky cavities (Figure 3). The SDM performance was even higher than for Eurosiberian forest-dwelling birds, traditionally benefit from the incorporation of descriptor of vegetation productivity (see e.g., [ 11 , 57 ]). These examples highlight the need for