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

: Urgent action needs to be taken to halt global biodiversity crisis. To be e ﬀ ective 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 di ﬀ erent life-history traits. We run ensemble species distribution models (SDMs) trained with bird atlas data and 12 EFAs describing di ﬀ erent 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 signiﬁcantly a ﬀ ect 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”).


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.
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][12][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][20][21][22][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, 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).

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]).
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 networkarea 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).

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]).  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/).

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.

EFA Name Acronyms Description
Annual mean Mean Surrogate of annual total amount.

Max
Annual maximum (indicator of the annual extremes).

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.

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][38][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.

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.

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). 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).

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][47][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][50][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][20][21][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 incorporating ecosystem functioning through remotely sensed descriptors of nutrient cycling and energy balance to characterize the niche of species and thus accurately predict their distribution.
Despite the large variability in the inter-annual habitat suitability at species level ( Figure S6), our SDMs predicted sharp declines in habitat suitability at community level in 2008 and 2014, and marked increases in 2004-05, 2010-11 and 2017 (Figures 4 and 5). These inter-annual dynamics were mostly driven by inter-annual variability in minimum EVI and maximum Albedo ( Figure S7), since they were the most important predictor variables across our bird community (see Figure 3). For instance, inter-annual habitat dynamics for forest-dwelling species such as the great spotted woodpecker (Dendrocopos major) or the Western Bonelli's warbler (Phylloscopus bonelli) were predicted to be affected by changes in vegetation productivity ( Figure S5), as described by the inter-annual trends of minimum EVI ( Figure S7). Other birds that depend on sparsely vegetated areas for feeding or on rocky cavities for nesting were predicted to change according to the inter-annual patterns in land-surface Albedo and temperature. In our study area, Mediterranean and migratory birds with preference for mountain habitats and rocky cavities were predicted to be the most negatively affected by environmental changes between 2001 and 2017-trends that are in line with those reported for other European mountains [58]. For instance, inter-annual changes in habitat suitability for the Eurasian crag martin-a species strongly dependent on rocky cavities for establishing breeding colonies-were mostly driven by changes in mean LST ( Figures S5 and S7). In fact, in addition to the inter-annual climate fluctuations, other factors can be also driving the inter-annual habitat dynamics for these species. Previous research demonstrated that land-surface Albedo and temperature are affected by burn severity [59,60]. Wildfires cause large exchanges in matter and energy that largely influence surface energy balance and temperature [61]. In our study, area, mean LST and minimum Albedo correlated with the area burnt in our study area between 2001 and 2017 ( Figures S7 and S8). However, other local processes such as land-use change or intensive forest harvesting (clearcutting; see Figure S2) could have also been affecting vegetation productivity and surface energy balance. Our model-assisted monitoring approach based on remotely sensed EFAs offer the opportunity to track bird habitat dynamics jointly driven by anthropogenic climate warming, regional land-use change, and local disturbances.

Conclusions
Our results demonstrated that remotely sensed descriptors of ecosystem functioning are good predictors of bird distribution in our region. Our ensemble species distribution models-exclusively based on ecosystem functioning attributes-hold a high predictive capacity across 71 bird species. Our results showed the life-history traits did not significantly affect SDM performance-even though the distribution was slightly better predicted for migratory and Mediterranean birds nesting in rocky cavities with omnivorous diets and vulnerable IUCN status. Our results underline the utility of remotely sensed descriptors of ecosystem functioning to characterize the ecological niche of species.
This study illustrates the great potential that satellite RS can contribute to the Aichi Biodiversity Targets and to the Essential Biodiversity Variables framework (EBV class "Species distribution"). Our model-assisted bird monitoring approach underlines the utility of remotely sensed descriptors of ecosystem functioning to accurately predict bird distribution at regional scale-in a cost-effective and standardized way. Our approach also provides an alternative method to monitor inter-annual bird habitat dynamics from space in the absence of monitoring programs based on repeated sampling designs (i.e., time series of biodiversity data). In addition, this model-based monitoring approach leverages the existing biodiversity data based on standardized surveys (e.g., atlas and inventory data) that rapidly become obsoleted due to the lack of long-term and well-established monitoring schemes. However, the lack of time series of biodiversity data also constrains our capacity to test the model ability to predict changes in species distribution beyond the range used to fit the model (i.e., model transferability). Given the large uncertainty when extrapolating the model beyond the range of data used for calibration, further research will be required to test the ability of SDMs based on RS-EFAs to predict species distribution over time.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/16/2549/s1, Appendix A: Total burnt area and land-use/cover change in the study area. Figure S1: Percentage of area burnt between 2001 and 2017 in the study area calculated from the Collection 6 MODIS Burned Area product. Figure S2: Land-use/cover changes between 2000 and 2014 estimated from Landsat-derived maps: (A) circular plot illustrating the land-cover type transitions, in million hectares (Mha), (B) area covered by each land-cover type per year (in hectares). The size of the lines is proportional in width to the contribution of each land-cover type to the change. The colors refer to the land-cover types. Figure S3: Maps depicting the distribution and breeding status of each bird species in the study area according to the EOAC (European Ornithologist Atlas Committee) codes. Species acronyms available in Table S1. Table S1: Table with the target species, their life-history traits (habitat preference, biogeographic origin, migratory status, nesting and feeding strategies) and IUCN conservation status. Table S1 also shows the number of occurrences, the predictive accuracy (values of AUC of ROC and TSS) of ensemble SDMs at species level, and the mean and standard deviation of predicted habitat changes between 2001 and 2017 estimated from raw habitat suitability (HS_ch) and after binarization (Distr_ch). Figure S4: Heatmap depicting pair-wise correlations among the 21 candidate EFA predictors. Values represent the Pearson correlation coefficient. See EFAs acronyms in Table S1. Figure S5: Accumulative variable importance across individual models for each modeling algorithm for the 71 modeled species. Figure S6: Inter-annual habitat dynamics predicted by the ensemble SDMs between 2001 and 2017 for each species. Figure S7: Inter-annual variability between 2001 and 2017 for the most important predictors of our bird community (namely minimum and maximum Albedo, minimum, and standard deviation of EVI, mean, and standard deviation of LST). For all plots, colored lines indicate mean values across the study area while the transparent colored areas indicate the standard deviation. Bar plots indicates proportion of annual burnt area calculated from the Collection 6 MODIS Burned Area product. Figure S8: Linear relationships (and Pearson correlation coefficients) between annual burnt area (calculated from the Collection 6 MODIS Burned Area product) and the averaged values of the most important EFA predictors (minimum and maximum Albedo, minimum, and standard deviation of EVI, mean, and standard deviation of LST). Appendix D: Google Earth Engine Code snippet for remotely sensed ecosystem functioning attributes.

Conflicts of Interest:
The authors declare no conflict of interest.