Remotely Sensed Variables of Ecosystem Functioning Support Robust Predictions of Abundance Patterns for Rare Species

remote sensing Article Remotely Sensed Variables of Ecosystem Functioning Support Robust Predictions of Abundance Patterns for Rare Species Salvador Arenas-Castro 1,2,* , Adrián Regos 1,3 , João F. Gonçalves 1, Domingo Alcaraz-Segura 4,5 and João Honrado 1,6 1CIBIO-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 2 CICGE—Centro de Investigaç ã o em Ci ê ncias Geo-Espaciais, Faculdade de Ci ê ncias, Universidade do Porto, Observatório Astronómico “Prof. Manuel de Barros”, Alameda do Monte da Virgem, 4430-146 Vila Nova de Gaia, Portugal 3Departamento de Zooloxía, Xenética e Antropoloxía Física, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain 4Department of Botany and Inter-University Institute for Earth System Research, University of Granada, 18071 Granada, Spain 5 Andalusian Center for the Assessment and Monitoring of Global Change (CAESCG), University of Almer í a, 04151 Almería, Spain 6Faculdade de Ciências, Universidade do Porto, 4169-007 Porto, Portugal *Correspondence: sar[email protected] Received: 7 August 2019; Accepted: 4 September 2019; Published: 6 September 2019   Abstract: Global environmental changes are affecting both the distribution and abundance of species at an unprecedented rate. To assess these effects, species distribution models (SDMs) have been greatly developed over the last decades, while species abundance models (SAMs) have generally received less attention even though these models provide essential information for conservation management. With population abundance defined as an essential biodiversity variable (EBV), SAMs could offer spatially explicit predictions of species abundance across space and time. Satellite-derived ecosystem functioning attributes (EFAs) are known to inform on processes controlling species distribution, but they have not been tested as predictors of species abundance. In this study, we assessed the usefulness of SAMs calibrated with EFAs (as process-related variables) to predict local abundance patterns for a rare and threatened species (the narrow Iberian endemic ‘Ger ê s lily’ Iris boissieri; protected under the European Union Habitats Directive), and to project inter-annual fluctuations of predicted abundance. We compared the predictive accuracy of SAMs calibrated with climate (CLI), topography (DEM), land cover (LCC), EFAs, and combinations of these. Models fitted only with EFAs explained the greatest variance in species abundance, compared to models based only on CLI, DEM, or LCC variables. The combination of EFAs and topography slightly increased model performance. Predictions of the inter-annual dynamics of species abundance were related to inter-annual fluctuations in climate, which holds important implications for tracking global change effects on species abundance. This study underlines the potential of EFAs as robust predictors of biodiversity change through population size trends. The combination of EFA-based SAMs and SDMs would provide an essential toolkit for species monitoring programs. Keywords: ecosystem functioning attributes (EFAs); essential biodiversity variables (EBVs); Iris boissieri; rare species; satellite remote sensing; species abundance models (SAMs); species distribution models (SDMs) Remote Sens. 2019,11, 2086; doi:10.3390/rs11182086 www.mdpi.com/journal/remotesensing Remote Sens. 2019,11, 2086 2 of 16 1. Introduction Human-induced global environmental change is threatening biodiversity globally at an unprecedented rate [ 1 ]. These environmental changes are expected to greatly impact both the distribution and the abundance of species. In that sense, both ‘species distribution’ and ‘population abundance’ are potential candidates to be considered as essential biodiversity variables (EBVs) in the EBV class ‘species populations’ [ 2 , 3 ]. Predictions from species abundance models (SAMs) can deliver relevant information both on species distributions and on population sizes for biodiversity monitoring and conservation management [ 4 ]. The integration of modeling tools and biodiversity monitoring protocols has increased our capacity to quantify and anticipate the response of species to environmental change across scales [ 5 – 7 ]. Model-assisted monitoring of distribution and abundances is therefore essential to assess the effectiveness of conservation policies globally [8–10]. However, while species distribution models (SDMs) have been greatly developed over the last decades [ 11 , 12 ] (and references therein), SAMs have generally received less attention [ 13 – 15 ], but see [ 16 ] and [ 17 , 18 ]. Moreover, these modeling techniques have been traditionally based on climate and static habitat features [ 19 ], failing to consider other environmental information that is ecologically relevant for species. Thus, the development and application of comprehensive environmental datasets and efficient modeling techniques are needed to predict and to monitor the geographic distribution of species and abundance, along with their related environmental drivers [ 20 ]. In this regard, satellite remote sensing time-series offer continuous and cost-effective measures of the environment across space and time [ 21 ], allowing near real-time assessment of species distribution and population dynamics [ 22 ]. This is especially relevant for rare and endangered species that often hold high habitat specificity, narrow niche breadth, limited dispersal ability, and increasing habitat fragmentation. As descriptors of exchange of matter and energy (i.e., ecosystem functioning [ 23 , 24 ]), and therefore processes that are key for species and their populations, satellite-derived ecosystem functioning attributes (EFAs) offer an integrative and earlier view of ecosystem responses to environmental changes than climate and landscape structural or compositional attributes [ 25 ]. Satellite-derived EFAs can, therefore, be used as process-related predictors in biodiversity models, since they provide information about ecosystem processes and properties related to the carbon (e.g., vegetation greenness) and water cycles (e.g., evapotranspiration) or with energy balance (e.g., land surface temperature and albedo). In fact, satellite-derived EFAs are being also tested as candidate EBVs related to the carbon cycle, energy, and radiation balance capable of informing about ecosystem components linked to species conservation status [ 26 , 27 ]. Thus, the incorporation of EFAs into SDMs was found to increase their predictive power and transferability [ 28 ], offering the opportunity to cost-effectively monitor multiple endangered species [ 27 , 29 ], at different spatial and temporal scales [ 30 ]. Still, despite these advantages, the predictive ability of EFAs in abundance models as well as to assess inter-annual population dynamics of rare species remains largely untested. The main goal of this study is to assess the potential of SAMs calibrated with satellite-derived EFAs (as process-related model predictors) to predict spatial patterns of abundance and to support predictions of the inter-annual abundance dynamics for rare plants, using as test species the narrow Iberian endemic ‘Ger ê s lily’ (Iris boissieri; protected under the EU Habitats Directive). Furthermore, we assess whether EFAs (as EBVs from the Ecosystem function class [ 2 ]) also hold the potential to inform on species abundance dynamics (Species populations EBV class) for biodiversity monitoring and reporting [ 4 , 26 ]. Sequentially, we: 1) assessed the ability of satellite-derived EFAs to support reliable predictions of Iris abundance; 2) analyzed the added-value of EFAs to replace or complement ‘traditional’ climate, land cover or topographic data in species abundance models; 3) projected (through back- and fore-casting) the best model to obtain yearly abundance predictions; and 4) compared those predicted abundances with disturbance factors (fire, climate anomalies) hypothesized to explain the predicted inter-annual abundance variations. The framework proposed here is flexible enough to be applicable to other species and socio-environmental contexts globally, and it is expected to assist conservation decision-making processes related to biodiversity monitoring and reporting schemes. Remote Sens. 2019,11, 2086 3 of 16 2. Materials and Methods 2.1. Test Species and Study Area We tested our approach with the ‘Ger ê s lily’ (Iris boissieri Henriq. =Xiphion boissieri (Henriq.) Rodion), a narrow-ranged endemic plant holding a ‘critically endangered’ conservation status [ 31 ], and protected under the European Habitats Directive (hereafter HD) for which EU member-states hold regular reporting obligations (under Article 17 of the HD/Annex II and IV) [ 32 ]. This species is restricted to mountainous areas of the northwest Iberian Peninsula (Figure 1a), where Portugal concentrates the largest populations of the species, especially in the Peneda-Ger ê s National Park (Figure 1b). This mountain-protected area represents a transitional area between the Mediterranean and Euro-Siberian and Alpine biogeographic regions and a mosaic of mixed vegetation types characterizes it. Average annual temperatures range from 17 ◦ C to 20 ◦ C, with some variation along the altitudinal gradient [ 33 ]. Highlands have an average temperature of about 10 ◦ C, ranging from 4 to 14 ◦ C, while areas in the valleys have milder climates with an average temperature of 14 ◦ C, ranging from 8 to 20 ◦ C. Total mean rainfall reaches 3000 mm per year (with more than 130 rainy days per year) and snowfall is frequent in the mountaintops. Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 18 2. Materials and Methods 2.1. Test Species and Study Area We tested our approach with the ‘Gerês lily’ (Iris boissieri Henriq. = Xiphion boissieri (Henriq.) Rodion), a narrow-ranged endemic plant holding a ‘critically endangered’ conservation status [31], and protected under the European Habitats Directive (hereafter HD) for which EU member-states hold regular reporting obligations (under Article 17 of the HD/Annex II and IV) [32]. This species is restricted to mountainous areas of the northwest Iberian Peninsula (Figure 1a), where Portugal concentrates the largest populations of the species, especially in the Peneda-Gerês National Park (Figure 1b). This mountain-protected area represents a transitional area between the Mediterranean and Euro-Siberian and Alpine biogeographic regions and a mosaic of mixed vegetation types characterizes it. Average annual temperatures range from 17 °C to 20 °C, with some variation along the altitudinal gradient [33]. Highlands have an average temperature of about 10 °C, ranging from 4 to 14 °C, while areas in the valleys have milder climates with an average temperature of 14 °C, ranging from 8 to 20 °C. Total mean rainfall reaches 3000 mm per year (with more than 130 rainy days per year) and snowfall is frequent in the mountaintops. Figure 1. Location of the study area (the Peneda-Gerês National Park) in the northwest Iberian Peninsula (a). Distribution of the sampled units of Iris boissieri in 2006 for each grid unit (250 × 250 m) in the study area (b). Figure 1. Location of the study area (the Peneda-Ger ê s National Park) in the northwest Iberian Peninsula ( a ). Distribution of the sampled units of Iris boissieri in 2006 for each grid unit (250 × 250 m) in the study area (b). The highest rocky outcrops with steep topography and low accessibility are the main habitat for the ‘Ger ê s lily’ [ 31 ]. The species typically occurs in open scrublands and crevices of granite outcrops at elevations between 500 and 1500 m above sea level. Many of these habitat areas were maintained by the traditional grazing system, which involved frequent small-scale burning for pasture maintenance Remote Sens. 2019,11, 2086 4 of 16 and was therefore an important control of vegetation encroachment. This long-term land management system has been declining for decades due to severe rural depopulation and other socioeconomic changes. As a result, widespread vegetation encroachment has changed habitat conditions as well as the local fire regime, potentially affecting the conservation of Iris boissieri and other plant species of high conservation concern. 2.2. Species Abundance Dataset The original dataset consisted of abundance estimates (counts of the number of individuals), performed in situ in the Peneda-Ger ê s National Park (Portugal) in 2006 by Park botanists. The field survey was conducted by visiting all previously known Iris occurrence sites as well as the neighboring areas. Iris clusters separated from each other by more than 100 meters were considered as distinct local populations to perform the counts. To match abundance data to the EFAs resolution (ca. 232 m MODIS pixel size; see below), abundance values were aggregated for each cell following an area-proportional allocation rule, i.e., the total number of plants within each sampling area (i.e., a local population defined by a polygon) was split between contiguous grid cells depending on the proportionally occupied area by the polygon in each cell. The final abundance dataset consisted of 61 cells (of the 232 m grid; Figure 1b), with the number of Iris plants per cell ranging from 1 to 285. 2.3. Modeling Framework To assess the predictive ability of satellite-derived EFAs in abundance models, we compared SAMs based on EFAs against SAMs based on climatic, topographic, or landscape variables. The best predictors were then combined to test whether more robust SAMs (and thus better abundance predictions) could be obtained. 2.3.1. Environmental Predictors We selected candidate predictive variables expressing environmental and ecological factors that are known to influence plant physiology, distribution, and local habitat preferences over a period (details below) as close as possible to the date of the in-situ abundance survey (2006). Based on a monthly time series of air temperature (mean maximum and mean minimum temperature) and total precipitation—which were calculated from meteorological stations in Spain and Portugal by using multiple regression techniques (see details in [ 34 ])—we derived 19 (bio-) climate variables at 200 m spatial resolution for each year from 2000 to 2010 by using the function ‘biovars’ available in the R package ‘dismo’, version 1.1. Land use/cover data was obtained by a hybrid classification procedure of Landsat-7 ETM +optical and thermal imagery (19 May 2010 and 30 July 2010; 30 m resolution), which combines unsupervised and supervised techniques (see details in [ 28 ]). From the classified land cover/use maps, we obtained the fractional cover of forest, shrubland, and rocky outcrop types to use as landscape composition variables at a spatial resolution of 232 m. From the digital elevation model (DEM) over Europe (EU-DEM based on SRTM and ASTER GDEM; https://www.eea.europa.eu/data-and-maps/data/eu-dem) produced at 30 m resolution, we derived topographic variables at 232 m grid size by using the function ‘terrain’ available in the R package ‘raster’, version 2.8–19. To compute EFAs, we used the MODIS Enhanced Vegetation Index (EVI) (MOD13Q1.v006; 232 m pixel every 16 days) for the 2001–2016 period. For that, we used Google Earth Engine [ 35 ] to derive the following eight summary metrics of the EVI seasonal dynamics: EVI annual mean (surrogate of annual primary production), EVI annual maximum and minimum (indicators of the annual extremes), EVI seasonal standard-deviation (descriptor of seasonality), and sine and cosine of the dates of maximum and minimum EVI [25,36] (indicators of phenology). Remote Sens. 2019,11, 2086 5 of 16 All environmental variables were processed to match the pixel size to the EFAs resolution (ca. 232 m grid size) at the geographical coordinate system WGS84, using the R package ‘raster’ (version 2.8–19). The initial dataset included 44 candidate predictors (Table S1—Supporting Information). 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). Based on these multicollinearity analyses, we finally retained three independent predictors for each competing model (with Pearson’s correlation coefficients of <0.7 and VIF of <4; Table S2; Figure S1). 2.3.2. Model Fitting For each model, the response variable (Iris abundance per grid cell; 61 records) was related to predictor variables using generalized linear models (GLMs) calibrated with three predictors. GLMs are flexible statistical methods for analyzing ecological relationships since they take into account the interactive behavior of variables that can be poorly represented by classical Gaussian (linear) distributions [ 37 ]. Furthermore, recent literature recommends using models based on Poisson and negative binomial distributions instead of log-transforming the data [ 38 ] with the main advantage that the parameters are estimated on the original scale and thus back-transformations are avoided [ 39 ]. Since this is an increasingly-used approach for count (abundance) data, we related Iris abundance values to predictors by computing a negative binomial distribution of errors (after checking for overdispersion effects in Poisson models and robust estimation via weighted likelihood), and applying a logarithmic link function by using the function ‘glm.nb’ available in the R package ‘MASS’, version 7.3–51.1. Based on a multi-model inference (MMI) framework [ 40 ], we tested and ranked the individual competing models fitted with three variables for abundance predictions performance (Table S2). We calibrated five GLMs as follows: (1) models based on process-related (i.e., ecosystem functioning) predictors (‘EFA’)—using only remotely sensed EFAs (their inter-annual mean for the 2001–2006 period); (2) climate-based models (‘CLI’)—using variables extracted from the monthly climate dataset (the inter-annual mean for the 2000–2010 period); (3) models based on land cover (landscape composition) predictors (‘LCC’)—from land cover maps (year 2010); (4) models based on topographic variables (‘DEM’); and (5) a best model, combining the most significant three predictors selected from all possible combinations of previous groups. The Akaike Information Criterion with a correction for finite sample size (AICc) was used to rank the competing models [ 41 ]. Complementarily, we used the explained deviance and Spearman’s correlation between observed and predicted values to assess model fitting. Models that were not significant, or with low performance in AICc or in explained deviance were excluded from the analysis. We also tested the predictive accuracy of each model by calculating the mean absolute scaled error (MASE) [ 42 ]. A MASE higher than 1 implies that the actual forecast is worse than a naive forecast as a reference method calculated as a sample, while a MASE lower than 1 implies that the current forecast behavior is better than the naive method. For that, we used the ´accuracy´ function through the R package ´forecast´ (version 8.9). 2.4. Inter-Annual Predictions of Iris Abundance The best model was backcasted and forecasted (i.e., projected) to the environmental conditions of each year to obtain a time-series of annual abundance predictions. Specifically, we used the most parsimonious and with greater predictive accuracy (‘best’) SAM (the one calibrated with the inter-annual (2001–2006) average values of EFAs plus topography; see Results) to obtain annual backward (from 2005 to 2001) and forward (from 2007 to 2016) yearly projections. We also calculated the inter-annual mean ( ¯ X ), the inter-annual standard deviation (SD), and the inter-annual coefficient of variation (CV; %) of predicted abundance at the grid cell level for the entire 2001–2016 time-period as a spatiotemporal proxy of the inter-annual variability of Iris populations. Finally, we explored possible relations between the inter-annual variability of species predicted abundance and the occurrence of disturbances, namely those related to fire and climate, by applying quantile regression analysis. Unlike ordinary least square regression, quantile regression is a Remote Sens. 2019,11, 2086 6 of 16 distribution-agnostic methodology, which makes no assumptions about the distribution of the residuals, and therefore not assume a parametric distribution for the response [ 43 ]. Thus, we explored different aspects of the relationship (as signals and effects) between the dependent variable (predicted abundance) and the independent variables (climate variables and anomalies). To do so, we calculated the total burned area per year from official wildfire maps based on Landsat-7 ETM+images [ 44 ] for the 2001–2016 period. Climatic means and annual anomalies were derived from TerraClimate [ 45 ] at 4 km resolution for the same period: total annual precipitation (mm · year −1 , PpT); average annual minimum temperature ( ◦ C; Tmin); and average annual maximum temperature ( ◦ C; Tmax). Since climate data, in general, have much more spatial autocorrelation, we resampled these climatic variables to ca. 232 m. The Iris abundance values showed a right-skewed distribution, and a prevalence of low abundance values makes it difficult to extract relations between population size and predictors in abundance trend analysis (heteroskedastic response [ 46 , 47 ]). To check the relationships and signals between the response variable (predicted abundance) and independent variables (climate and its anomalies), we applied quantile regression tests (‘quantreg’ R package) fitted to the median (quantile =0.5). To test the goodness-of-fit for quantile regression models, we used the pseudo-R 2 measure suggested by [ 48 ]. Furthermore, we tested the same associations between inter-annual predicted abundances and fire and climate variables, but for different levels of population size (quantiles). For that, we applied the same quantile regression test where abundances were separated by 5% intervals (i.e., in quantiles) in order to overcome limitations related to the right-skewed distribution and heteroskedastic response of abundance data. To check significant differences across quantiles, we used the test of equality of slopes by using the ´ANOVA´ function implemented by ‘quantreg’ developers. All analyses were conducted using the R environment for statistical computation [ 49 ]. QGIS version 2.18.11 and ArcGIS version 10.2 were used for managing and representing spatial data and projections. 3. Results 3.1. Ranking of Models and Predictors Overall, models that include EFAs as predictors showed the highest performance (explained deviance up to 0.3), whereas the LCC-based model attained the lowest performance (Table 1). Table 1. Results from the generalized linear models (GLM) model selection and multi-model inference explaining local abundance patterns of Iris boissieri. The competing models are listed in descending order from the best to least fitted model determined by their mean absolute scaled error (MASE) values. For comparison, a baseline ‘null model’ containing a single intercept term is used. Loglik: loglikelihood. AICc: Akaike Information Criterion value. ∆ AICc: represents the differences between the AICc values of the best model considered and other competing models. wi: represents the Akaike weights and indicates the probability that a particular model is the best among those considered. Explained deviance: the proportion of the explained deviance by the model. Spearman Correlation: Spearman correlation between observed and predicted values. MASE: Mean absolute scaled error as a measure of the predictive accuracy, and ranges from 0 =highest accuracy, to 1 =lowest accuracy. Competing Model Predictors LogLik AICc ∆AIC wi Explained Deviance Spearman Correlation MASE EFAs +DEM - EVImin - EVIdmic - SLO −239.68 487.95 0.00 0.79 0.33 0.44 0.73 EFAs +CLI - EVImin - EVIdmic - BIO5 −240.36 490.71 2.14 0.19 0.32 0.43 0.75 EFAs - EVImin - EVImn - EVIdmic −240.87 492.85 2.76 0.20 0.31 0.43 0.77 Remote Sens. 2019,11, 2086 7 of 16 Table 1. Cont. Competing Model Predictors LogLik AICc ∆AIC wi Explained Deviance Spearman Correlation MASE DEM - ELE - SLO - ASP −243.89 498.52 8.43 0.01 0.24 0.34 0.90 CLI - BIO5 - BIO15 - BIO19 −247.62 505.97 15.88 0.00 0.15 0.31 0.96 LCC - DeFo - OShr - RoAr −248.77 508.28 18.19 0.00 0.12 0.09 1 Null model - −253.38 510.97 20.88 0.00 0.00 - 1 The most parsimonious model (based on AICc) and with the greatest predictive ability (based on MASE) was obtained when combining EFAs and topographic variables (EFAs +DEM model; Table 1). The correlation between species abundance and predictor variables in this model (EVIdmic, EVImin, and SLO) was negative for the two EFAs and positive for SLO (Spearman correlations were − 0.13, − 0.37, and 0.10, respectively); however, only EVImin showed a significant correlation with Iris abundance (Figure 2). Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 18 Figure 2. Relation between observed abundance (number of individuals) of Iris boissieri and each of the predictive variables included in the most parsimonious species abundance model (SAM) (Ecosystem functioning attributes (EFAs) + digital elevation model (DEM)); (a) cosine of the date of minimum enhanced vegetation index (EVI) (EVIdmic as minimum vegetation greenness); (b) EVI annual minimum (EVImin as indicator of minimum photosynthetic activity); (c) Slope (SLO as an indicator of terrain flatness). A significant Spearman correlation (R) was only found for the relationship between observed abundance and EVImin. 3.2. Inter-Annual Variability of Iris Abundance On average (considering all pixels in the study area), the median and median absolute deviation (MAD) of abundance predictions based on EFAs projections over the full 2001–2016 time-period were of 17.43 (12 for observed values) and 12.04 (14.08) individuals, respectively, ranging from 12.79 (11.96 MAD) in year 2006 to 23.85 (11.84 MAD) in year 2007 (Figure S2). Overall, the inter-annual abundance predictions were spatially consistent with the known locations of Iris (Figure S3). The inter-annual mean of the predicted abundances for the 61 observed sampled cells based on the predicted spatial X , SD and CV were 25.02, 8.06 and 33.06 individuals, respectively. Spearman correlation analysis showed a high and positive significant correlation between the X  and SD (r = 0.83, p < 0.05) (Figure S4). There was also a high correlation of X with slope (r = 0.81, p < 0.05) and with EVImin (r = -0.33, p < 0.05), but no correlation with EVIdmic. The CV showed no correlation with slope (r = -0.07, p < 0.05), but a correlation with EVImin (r = 0.60, p < 0.05) and EVIdmic (r = - 0.23, p < 0.05). Those grid cells that combine higher values for the X , and lower values for the CV, therefore represent those areas where conditions are better to maintain large populations (Figures 3a, b). Figure 2. Relation between observed abundance (number of individuals) of Iris boissieri and each of the predictive variables included in the most parsimonious species abundance model (SAM) (Ecosystem functioning attributes (EFAs) +digital elevation model (DEM)); ( a ) cosine of the date of minimum enhanced vegetation index (EVI) (EVIdmic as minimum vegetation greenness); ( b ) EVI annual minimum (EVImin as indicator of minimum photosynthetic activity); ( c ) Slope (SLO as an indicator of terrain flatness). A significant Spearman correlation (R) was only found for the relationship between observed abundance and EVImin. 3.2. Inter-Annual Variability of Iris Abundance On average (considering all pixels in the study area), the median and median absolute deviation (MAD) of abundance predictions based on EFAs projections over the full 2001–2016 time-period were Remote Sens. 2019,11, 2086 8 of 16 of 17.43 (12 for observed values) and 12.04 (14.08) individuals, respectively, ranging from 12.79 (11.96 MAD) in year 2006 to 23.85 (11.84 MAD) in year 2007 (Figure S2). Overall, the inter-annual abundance predictions were spatially consistent with the known locations of Iris (Figure S3). The inter-annual mean of the predicted abundances for the 61 observed sampled cells based on the predicted spatial ¯ X , SD and CV were 25.02, 8.06 and 33.06 individuals, respectively. Spearman correlation analysis showed a high and positive significant correlation between the ¯ X and SD (r =0.83, p<0.05) (Figure S4). There was also a high correlation of ¯ X with slope (r =0.81, p <0.05) and with EVImin (r = − 0.33, p <0.05), but no correlation with EVIdmic. The CV showed no correlation with slope (r = − 0.07, p <0.05), but a correlation with EVImin (r =0.60, p <0.05) and EVIdmic (r = − 0.23, p<0.05). Those grid cells that combine higher values for the ¯ X , and lower values for the CV, therefore represent those areas where conditions are better to maintain large populations (Figure 3a,b). Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 18 Figure 3. Spatial projections of predicted Iris boissieri abundance representing the variability of the species populations based on: (a) inter-annual mean; and (b) inter-annual coefficient of variation (CV) for the whole 2001–2016 time-period based on the EFAs + SLO model. Blank cells represent values < 30. The overall mean for the study area of the annual predicted abundance matched the indicator of minimum photosynthetic activity (EVImin) (Figure 4a), and was higher after rainy years (total average annual precipitation; PpT), and mainly matched chilly years (lower values of maximum (Tmax) and minimum (Tmin) annual temperature) (Figure 4b). More specifically, climate conditions of the year with the lowest predicted abundance (16.79; 2002) were PpT = 1724.97 mm.year -1 , Tmax = 15.69 °C, and Tmin = 7.38 °C, while for the year with the largest predicted abundance (23.85; 2007) were PpT = 1087.49 mm year-1, Tmax = 15.88 °C, and Tmin = 7.21 °C. However, wildfires (total burnt area; TBA) did not have a significant effect on predicted Iris abundances (R 1 (pseudoR 2 ) = 0.025; Figure 4c; Figure S5). Figure 3. Spatial projections of predicted Iris boissieri abundance representing the variability of the species populations based on: ( a ) inter-annual mean; and ( b ) inter-annual coefficient of variation (CV) for the whole 2001–2016 time-period based on the EFAs +SLO model. Blank cells represent values <30. The overall mean for the study area of the annual predicted abundance matched the indicator of minimum photosynthetic activity (EVImin) (Figure 4a), and was higher after rainy years (total average annual precipitation; PpT), and mainly matched chilly years (lower values of maximum (Tmax) and minimum (Tmin) annual temperature) (Figure 4b). More specifically, climate conditions of the year with the lowest predicted abundance (16.79; 2002) were PpT =1724.97 mm · year −1 , Tmax =15.69 ◦ C, and Tmin =7.38 ◦ C, while for the year with the largest predicted abundance (23.85; 2007) were PpT =1087.49 mm year-1, Tmax =15.88 ◦ C, and Tmin =7.21 ◦ C. However, wildfires (total burnt area; TBA) did not have a significant effect on predicted Iris abundances (R 1 (pseudoR 2 )=0.025; Figure 4c; Figure S5). The quantile regression analysis between average annual predicted abundance and inter-annual climate fluctuations (in abundance levels across time), and anomalies (changes in climatic variables in relation to the mean), showed negative relationships for total precipitation (PpT; R 1 (pseudoR 2 )=0.30), and its anomalies (PpTa; R 1 =0.30), as well as minimum temperature (Tmin; R 1 =0.27) and its anomalies (Tmina; R1=0.28) (Figure S6). Remote Sens. 2019,11, 2086 9 of 16 Remote Sens. 2019, 11, x FOR PEER REVIEW 10 of 18 Figure 4. Yearly average values of predicted abundance of Iris boissieri (black line) in the study area for the time-period 2001–2016 against annual variables: (a) EFAs through EVI minimum (EVIm) and date of minimum EVI; (b) climate through total precipitation (Pp), maximum (Tmax) and minimum (Tmin) temperature; (c) total amount of burned area (TBA). Dark grey bars represent the calibration for the period 2001–2006. Red lines represent lag response of the species to fire events. Figure 4. Yearly average values of predicted abundance of Iris boissieri (black line) in the study area for the time-period 2001–2016 against annual variables: ( a ) EFAs through EVI minimum (EVIm) and date of minimum EVI; ( b ) climate through total precipitation (Pp), maximum (Tmax) and minimum (Tmin) temperature; ( c ) total amount of burned area (TBA). Dark grey bars represent the calibration for the period 2001–2006. Red lines represent lag response of the species to fire events. Remote Sens. 2019,11, 2086 16 of 16 54. Pearson, R.G.; Dawson, T.P. Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Glob. Ecol. Biogeogr. 2003,12, 361–371. [CrossRef] 55. Deblauwe, V.; Droissart, V.; Bose, R.; Sonk é , B.; Blach-Overgaard, A.; Svenning, J.-C.; Wieringa, J.J.; Ramesh, B.R.; St é vart, T.; Couvreur, T.L.P. Remotely sensed temperature and precipitation data improve species distribution modelling in the tropics. Glob. Ecol. Biogeogr. 2016,25, 443–454. [CrossRef] 56. Lassueur, T.; Joost, S.; Randin, C.F. Very high resolution digital elevation models: Do they improve models of plant species distribution? Ecol. Modell. 2006,198, 139–153. [CrossRef] 57. Cord, A.F.; Meentemeyer, R.K.; Leit ã o, P.J.; V á clav í k, T. Modelling species distributions with remote sensing data: bridging disciplinary perspectives. J. Biogeogr. 2013,40, 2226–2227. [CrossRef] 58. Cord, A.F.; Klein, D.; Gernandt, D.S.; de la Rosa, J.A.P.; Dech, S. Remote sensing data can improve predictions of species richness by stacked species distribution models: A case study for mexican pines. J. Biogeogr. 2014 , 41, 736–748. [CrossRef] 59. Alcaraz, D.; Paruelo, J.; Cabello, J. Identification of current ecosystem functional types in the iberian peninsula. Glob. Ecol. Biogeogr. 2006,15, 200–212. [CrossRef] 60. Hogrefe, K.R.; Patil, V.P.; Ruthrauff, D.R.; Meixell, B.W.; Budde, M.E.; Hupp, J.W.; Ward, D.H. Normalized difference vegetation index as an estimator for abundance and quality of avian herbivore forage in Arctic Alaska. Remote Sens. 2017,9, 1234. [CrossRef] 61. Ames, G.M.; Wall, W.A.; Hohmann, M.G.; Wright, J.P. Trait space of rare plants in a fire-dependent ecosystem. Conserv. Biol. 2017,31, 903–911. [CrossRef] 62. Pausas, J.G.; Keeley, J.E. Wildfires as an Ecosystem Service. Front. Ecol. Environ. 2019 ,17, 289–295. [CrossRef] 63. Renwick, A.R.; Massimino, D.; Newson, S.E.; Chamberlain, D.E.; Pearce-Higgins, J.W.; Johnston, A. Modelling changes in species’ abundance in response to projected climate change. Divers. Distrib. 2012 ,18, 121–132. [CrossRef] © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).