Predicting Microhabitat Suitability for an Endangered Small Mammal Using Sentinel-2 Data

: Accurate mapping is a main challenge for endangered small-sized terrestrial species. Freely available spatio-temporal data at high resolution from multispectral satellite o ﬀ er excellent opportunities for improving predictive distribution models of such species based on ﬁne-scale habitat features, thus making it easier to achieve comprehensive biodiversity conservation goals. However, there are still few examples showing the utility of remote-sensing-based products in mapping microhabitat suitability for small species of conservation concern. Here, we address this issue using Sentinel-2 sensor-derived habitat variables, used in combination with more commonly used explanatory variables (e.g., topography), to predict the distribution of the endangered Cabrera vole ( Microtus cabrerae ) in agrosilvopastorial systems. Based on vole surveys conducted in two di ﬀ erent seasons over a ~176,000 ha landscape in Southern Portugal, we assessed the signiﬁcance of each predictor in explaining Cabrera vole occurrence using the Boruta algorithm, a novel Random forest variant for dealing with high dimensionality of explanatory variables. Overall, results showed a strong contribution of Sentinel-2-derived variables for predicting microhabitat suitability of Cabrera voles. In particular, we found that photosynthetic activity (NDI45), speciﬁc spectral signal (SWIR1), and landscape heterogeneity (Rao’s Q) were good proxies of Cabrera voles’ microhabitat, mostly during temporally greener and wetter conditions. In addition to remote-sensing-based variables, the presence of road verges was also an important driver of voles’ distribution, highlighting their potential role as refuges and / or corridors. Overall, our study supports the use of remote-sensing data to predict microhabitat suitability for endangered small-sized species in marginal areas that potentially hold most of the biodiversity found in human-dominated landscapes. We believe our approach can be widely applied to other species, for which detailed habitat mapping over large spatial extents is di ﬃ cult to obtain using traditional descriptors. This would certainly contribute to improving conservation planning, thereby contributing to global conservation e ﬀ orts in landscapes that are managed for multiple purposes.


Introduction
Anthropogenic activities, concurrently with human population growth, are responsible for wiping out wildlife species at rates never experienced before [1]. In particular, agricultural intensification and infrastructure proliferation (roads, railways, etc.), which are considered among the main causes of habitat loss/fragmentation and populations declines, have been rapidly rising to an alarmingly level worldwide [2,3]. Traditionally, wildlife conservation priorities have been focused on megafauna, since species with a large body size have been associated with high extinction risks [4]. However, small body size can also be an important extinction driver [5], possibly exacerbated by species limiting ecological traits (e.g., short dispersal distances), restricted, and/or fragmented distribution and habitat specialization [6].
The Cabrera vole (Microtus cabrerae) is an Iberian-endemic small mammal, classified as "Vulnerable" in Portugal and Spain [7,8], and as "Near-threatened" by IUCN [9]. Within its restricted distribution range, the species presents a fragmented distribution [10], typically associated with marginal areas of agricultural systems, with local populations largely restricted [10][11][12] to sparse patches of tall and dense wet grasslands [11,13]. The major threats for this species include agriculture and grazing intensification [14], which destroy its preferred habitats, forcing individuals to disperse and occupy small habitat patches (often <500 m 2 [14]) like field margins or road verges [12,15,16]. The Cabrera vole often presents a metapopulation-like spatial structure, which together with the regular destruction and turnover of suitable habitat patches, makes the designation of special areas of conservation for this species a particularly challenging task. The designation of these conservation areas is however demanded by the European Union, as the species is listed in both Bern Convention (Appendix II; 82/72/CEE) and Habitats Directive (Annexes II and IV; Council Directive 92/43/EEC). The selection of those key areas should be supported by detailed and up-to-date species' distribution at multiple scales, and the use of efficient tools and frameworks able to appropriately identify them [17]. In this context, correlative species distribution models (SDMs), or habitat suitability/niche models [18], which provide probabilistic estimation of occurrence patterns over broad areas by relating species occurrences with environmental characteristics [18], have become a popular tool to develop potential species range maps.
Numerous studies have extensively reported the utility of SDMs for addressing a variety of ecological questions [19][20][21], related to biodiversity monitoring and conservation planning [17,22,23], including for the Cabrera vole [10,24]. Yet, SDMs applications on Cabrera voles, or other small and elusive species, at a local or regional scale are still challenging, likely due to their low detectability and/or narrow distribution, which may complicate data collection [25,26]. Moreover, the integration of fine grain habitat requirements for which suitability may change within short time periods makes SDM' building another challenging task, due to the lack of spatially explicit predictor variables able to capture habitat characteristics at small scales [27], as well as to account for species occupancy turnover and landscape dynamism [28], the latter being markedly pronounced in Mediterranean-type ecosystems [12,29,30]. Specifically, most available digital habitat proxy information (e.g., land cover/use maps) have low detail precision and have a static time nature (they are not expected to vary within the year) [28,31], and thus may fail to provide relevant ecological information for small species inhabiting dynamic habitat patch networks.
We used Cabrera vole as a model to create up-to-date spatially and temporally detailed habitat suitability maps for species with fine-scale habitat requirements occurring in dynamic landscapes. Opportunities to do this come from Earth Observation Satellites (EOS) due to their multispectral and systematic characteristics, which allows the identification of the vegetation composition and structure, as well as its physiological condition [32][33][34]. The usefulness of remote-sensing data for species habitat suitability mapping has been reported in numerous studies, as outlined in the review by He et al. [27]. In this review, the spatial-continuous nature and the reasonable time frequency of satellite-based data are highlighted as an added value to overcome SDMs limitations. By integrating this high-quality data into SDMs, their accuracy can be effectively increased as availability of resources may be better described [28,[35][36][37]. Moreover, remote-sensing data can be used for modelling changes in species distribution across time and understand how vegetation changes might affect patch quality and influence demographic parameters, including reproduction and dispersal movements [28].
While it may be straightforward to map habitat suitability areas, for example, for large mammal species [35,38], having broad-scale home range sizes (e.g., >1000 m 2 ), modelling species responding to fine-scale landscape requirements (e.g., small mammals or insects) is challenging from the remote-sensing perspective due to limitations associated to conventional imageries when identifying local resource patches [38,39]. Indeed, until recently, the available information from remote sensing (e.g., land-cover) was too coarse or too expensive to be properly applied on fine-scale modelling [28]. The Copernicus Program from the European Commission (EC) in partnership with the European Space Agency (ESA) has been developing several satellite missions under the scope of the Sentinel program [40]. Within this program, a constellation of two multi-spectral satellites called Sentinel-2A (launched on 23 June 2015) and Sentinel-2B (launched on 7 March 2017) are together collecting information at high spatial (up to 10 m), spectral (13 bands), radiometric (12 bits), and temporal (each five days) resolution [41]. Due to its technical features and the open data policy, Sentinel-2 brings new opportunities and capabilities for evaluating wildlife spatio-temporal response to habitat features [27] and dynamic processes [36], which may be of particular importance for SDMs developed for small species inhabiting dynamic systems (e.g., grasslands [42]) such as the Cabrera vole. To the best of our knowledge, modelling fine-scale habitat suitability for wildlife conservation, specifically with open-access remote-sensing data and with Sentinel-2 imagery, is still scarce in the literature. Besides, as Sentinel-2 derived-products mostly reflect biotic environmental attributes, the integration of these variables with abiotic descriptors (e.g., topography) into SDMs likely provide more realistic results than using each type of variables alone [28,36,43].
Therefore, by taking advantage from spectral, temporal, and spatial characteristics of Sentinel-2 sensors, the main goal of this study is to assess the usefulness of Sentinel-2 derived predictors for identifying suitable microhabitats for small and elusive species of conservation concern, using the Cabrera vole in a Mediterranean ecosystem as a model. In particular, we aimed to: i.
Quantify the importance of Sentinel-2 derived predictors relative to more conventional predictors (e.g., topographical and distance to landscape elements) in predicting vole microhabitat suitability; ii. Identify which Sentinel-2 derived predictors best explain vole distribution at fine spatial scales.
Overall, we predict that Sentinel-2-based variables should provide an important contribution for improving fine-scale habitat mapping of endangered small species, such as the Cabrera vole, thus supporting the view that remote-sensing products should greatly contribute for conserving biodiversity associated to small marginal areas in human-dominated landscapes. For this purpose, a methodological approach was devised for predicting suitable habitat areas for the Cabrera vole by using Boruta Random Forest algorithm [44] and different Sentinel-2-derived data (multispectral data, spectral indices, and textural and diversity indices), topographic variables, and distance to landscape key elements (roads, built-up areas, and water ponds).

Study Area
The study was conducted in a~176,000 ha area located in the Alentejo region, Southern Portugal (centroid: 586545 -4281192; EPSG: 32629-WGS 84/UTM 29N; Figure 1a). The area is characterized by an altitude ranging from 80 to 500 m a.s.l. with a gently undulating relief [and included within a bioclimatic zone commonly associated to the Cabrera vole, namely the meso-Mediterranean,29]. Climate is typically Mediterranean, with hot and dry summers (August: 31 • C Tmax), mild and wet winters (January: 6 • C Tmin), and medium annual rainfall (>600 mm) (Évora 1981-2010 [45]). The landscape is largely dominated by an agrosilvopastoral system called montado (or dehesa), an open woodland of cork (Quercus suber) and/or holm oak (Quercus rotundifolia) trees [46]. The system is characterized by high spatial variability in tree density and an understorey mosaic of annual crops, grasslands (intermixed perennial and annual herb communities), and shrublands [47]. While the montado is considered as one of the highest biodiversity-rich ecosystems of the western Mediterranean Basin [48,49] having been classified as a High Nature Value farming system (HNV) [50], it is also referred as one the most threatened in terms of conservation, mainly due to land use intensification [51]. grasslands (intermixed perennial and annual herb communities), and shrublands [47]. While the montado is considered as one of the highest biodiversity-rich ecosystems of the western Mediterranean Basin [48,49] having been classified as a High Nature Value farming system (HNV) [50], it is also referred as one the most threatened in terms of conservation, mainly due to land use intensification [51].

Cabrera Vole Field Surveys
Cabrera vole surveys were conducted through stratified random selection by initially identifying in the field suitable and unsuitable grass patches. A total of 146 patches with dense and tall perennial grasses and/or sedge/rush communities growing in high soil moisture conditions [13,14] were defined as locations of potential occurrence and 79 patches were considered not suitable for the species, due to very dry soil conditions and/or lower cover and height. Each of the selected patches was carefully surveyed by two observers for presence signs typical of this species (surface runways, grass clippings, and typical small, dark green faeces associated with latrines) to assess its presence [14]. These signs are easily recognizable, and together provide a reliable sampling method, at least when other species producing similar signs (e.g., M. agrestis) are absent in the area [11,13], as it is the case of our study region. Each surveyed habitat patch was classified according to the presence/absence of the species, and georeferenced with an accurate GPS device (Garmin eTrex 30x; Projected coordinate system: EPSG: 32629-WGS 84 / UTM 29N; precision up to 3 m). The absences were further classified as absences with and without suitable habitat conditions (as the first ones may correspond to patches potentially used by voles, but that were not occupied at the time of the survey); [12]. Although each patch was surveyed once, samplings were conducted in two sessions to account for habitat variation, namely soil humidity, vegetation dryness, and structure. The first session ran in Spring (February-April 2017), which is when Cabrera vole populations are typically close to their

Cabrera Vole Field Surveys
Cabrera vole surveys were conducted through stratified random selection by initially identifying in the field suitable and unsuitable grass patches. A total of 146 patches with dense and tall perennial grasses and/or sedge/rush communities growing in high soil moisture conditions [13,14] were defined as locations of potential occurrence and 79 patches were considered not suitable for the species, due to very dry soil conditions and/or lower cover and height. Each of the selected patches was carefully surveyed by two observers for presence signs typical of this species (surface runways, grass clippings, and typical small, dark green faeces associated with latrines) to assess its presence [14]. These signs are easily recognizable, and together provide a reliable sampling method, at least when other species producing similar signs (e.g., M. agrestis) are absent in the area [11,13], as it is the case of our study region. Each surveyed habitat patch was classified according to the presence/absence of the species, and georeferenced with an accurate GPS device (Garmin eTrex 30x; Projected coordinate system: EPSG: 32629-WGS 84 / UTM 29N; precision up to 3 m). The absences were further classified as absences with and without suitable habitat conditions (as the first ones may correspond to patches potentially used by voles, but that were not occupied at the time of the survey); [12]. Although each patch was surveyed once, samplings were conducted in two sessions to account for habitat variation, namely soil humidity, vegetation dryness, and structure. The first session ran in Spring (February-April 2017), which is when Cabrera vole populations are typically close to their peaks and breeding activity is presumably higher, due to increased soil humidity and vegetation growth (e.g., green grasses) [14]. The second session was conducted in autumn (October-early December 2018); when soil humidity was significantly lower due to the typical hot and dry summer conditions in the region, which were exceptionally hard and extended in 2018 (IPMA Évora 2018 [45]). This second session was also coincident with the period when more fallow areas can be found, being those of special interest for species' conservation [52]. A total of 97 and 128 herbaceous patches were surveyed in the first and second sessions, respectively. In order to lower model biases, all absences recorded in habitats identified as suitable were discarded from the dataset, as these may have resulted from possible low detectability [18,26]. We further applied a 500 m grid spatial filtering procedure, resulting in a roughly balanced dataset of 62 presences and 79 absences (Figure 1b).

Sentinel-2 Derived Predictor Variables
To better assess the capability of Sentinel-2 imagery in predicting Cabrera vole habitat suitability areas, three different types of Sentinel-2-derived variables were used: (1) Spectral bands, (2) spectral indices, and (3) textural and diversity indices.
Sentinel-2 multispectral images (Sentinel-2A MSI Level-1C) used in this study were downloaded from the Copernicus Science Data Hub portal (https://scihub.copernicus.eu/dhus/). For each of the study periods, the image with the lowest percentage of clouds was selected to represent environmental conditions at the time of vole surveys (5th April 2017 and 7th October 2018 in the case of the first and second period, respectively). The study area was entirely covered by the union of 4 multispectral images (0%-1% of clouds) for each selected period, which followed an atmospheric correction procedure using the Sen2Cor code implemented in the SNAP software [53].
In order to increase the spatial resolution of the 20 m spectral bands, a super-resolution enhancement method was applied, whereby high-resolution bands (10 m) were able to reconstruct coarser (20 m) at the given resolution while maintaining the associated spectral reflectance, as demonstrated by Brodu [54]. Super-resolved (SR) bands were computed using the Sen2res SNAP plugin (http://step.esa. int/main/third-party-plugins-2/sen2res/).
In order to capture different habitat features that are ecologically relevant to predict suitable areas for the Cabrera vole, three groups of spectral indices were computed: (1) Vegetation biomass indices (NDVI, NDRE1, NDRE2, NDRE3, NDI45, and SATVI), (2) senescent vegetation and soil surface indices (PSRI, SWIR32, and BI2), and (3) vegetation and landscape water content indices (NDII and NDWI) ( Table 1). These indices have been successfully used in retrieving different key biophysical vegetation information in semi-arid tree-grass ecosystems such as the one here addressed (montado) [55][56][57][58][59]. Calculated using the NDVI with a 3 × 3 pixels spatial moving window [70] To describe the montado vegetation and landscape structural and diversity properties, the grey-level co-occurrence matrix (GLCM) [69] and the Rao's Q index [36,70] were calculated, respectively. Prior to the textural calculation, the previously selected spectral bands underwent a Principal Component Analysis (PCA) fusion technique with the aim of obtaining a single Sentinel-2 image incorporating all bands' information [71]. The principal component image accounting for over the 90% of bands spectral variability was subsequently used to compute eight GLCM variables, namely, mean, correlation, contrast, Dissimilarity, entropy, homogeneity, second moment, and variance ( Table 1). The selected textural variables were calculated using the glcm package (v.1.6.1) [72] implemented in the R (v. 3.5.2) [73], and following the same parametrization settings described in Godinho et al. [57]. The Rao's Q diversity index, which accounts for both the abundance and the pairwise spectral distance among pixels [70], and thus is useful to assess spatial diversity, was calculated by using NDVI as input data and a moving window size of 3 × 3 pixels.

Distance to Landscape Elements
In order to quantify the potential influence of key landscape elements on Cabrera vole spatial distribution (e.g., [14,52]), distances to paved roads, built-up areas, and water bodies were calculated ( Table 2). A shapefile containing the information about paved roads was produced using OpenStreetMap data source [76]. Built-up areas and water bodies shapefiles were obtained from the imperviousness and the water and wetness high-resolution layers of the Copernicus Land Monitoring Service [77].

Habitat Suitability Model
The habitat suitability model was built using all previously described predictors using Cabrera vole presence/absence as response variable. The relationship between the predictors and the spatial distribution of Cabrera vole was evaluated in a three-step statistical approach. The first step consisted in selecting the relevant variables from a set of 67 candidate predictors using the Boruta algorithm [44,78,79]. Basically, Boruta algorithm relies on an extension of the random forest (RF) [80,81] method by introducing an iterative procedure to compare the relative importance of the original variables with the importance of their randomized copies [44]. After running iteratively a large number of random forest models, the Boruta algorithm computes the mean Z-score value to classify all the variables as confirmed, rejected, or tentative at a predefined threshold of statistical significance (p) and a maximum number of times the algorithm is run (maxRuns) [79]. In this study, the Boruta R package (v.6.0.0) [44] was used to execute the algorithm with maxRuns = 2000, ntree = 2000, and p value = 0.01. The second step consisted of running a Pearson's correlation analysis to determine pairwise correlations within the variables classified as confirmed in the previous step to remove highly correlated (r > |0.7|) ones. Finally, in the third step, and employing only the uncorrelated most important variables, an RF analysis was used to predict the spatial distribution of Cabrera vole in the study area. For the RF model, the number of trees (ntree) was fixed to 2000 and number of variables randomly tested on each split (mtry) to the square root of the number of variables. A 10-fold cross-validation resampling method was used to build the RF model. These analyses were done with the ggRandomForest R package (v.2.0.1) [82]. Each variable relative importance for the model was assessed and partial dependence plots [81] were used to explore interaction effects between variables on Cabrera vole presence probability. Model performance was verified using the area under the curve (AUC) of the Receiver Operator Characteristic (ROC), as well as the proportion of correctly predicted presences and absences [83].

Model Performance
The Boruta screening procedure resulted in a considerable reduction of possible explanatory variables, as only 26 predictors were confirmed (38.8% of all the candidate features set, Figure S1). From these, only 11 showed no strong correlation among them (r < |0.7|) and were retained for the multivariate analysis ( Figure S2; for more details regarding all pairwise correlation results, see Table S1). The results derived by the 10-fold cross-validation indicated that the RF model developed was robust given the low estimated error rate percentage, (19.15%), determining a high explanatory power of included predictors on the occurrence of the endangered Cabrera vole in our study area (about 80% of variance explained). Results also showed a 'high' AUC score (area under the curve) of 0.904, a sensitivity (true positive rate) of 0.73, and a specificity (true negative rate) of 0.778, therefore a higher performance for correctly predicted absences than presences was noticed.

Variable Importance
Following the multivariate analysis, the "Sentinel-2" variables group showed the highest contribution (65.7%) in explaining Cabrera vole habitat suitability, comprising 10 variables (Figure 2). The variables from the group "Distance to landscape elements" contributed to explain 34.22% of the variance, comprising only the distance to paved roads ( Figure 2). None of the "Topographic" variables were retained in the final model. Half of "Sentinel-2" variables concerned the Spring period and another half to the Autumn period ( Figure 2). The highest significant contributors from the "Sentinel-2" group were "NDI45 (Spring)" (14.9%), "SWIR1 (Autumn)" (10.4%), and "Rao's Q (Spring)" (9.9%) (Figure 2), meaning these variables incorporated most of the relevant habitat information from remote-sensing data. The habitat suitability for Cabrera vole increased when the spectral vegetation index NDI45 had low-medium values in Spring, and the spectral band SWIR1 and the metric Rao's Q showed intermediate values in Autumn and Spring, respectively (Figure 3c,d). Response curves for "Distance to paved roads" showed that suitability of Cabrera vole steeply decreased with the increase in distance from roads (Figure 3a). The habitat suitability map shows that the occurrence locations fell in high-probability areas in the final habitat suitability model (Figure 4).

Discussion
Results yielded evidence that fine-scale remote-sensing data may be useful to predict favorable habitats for the occurrence of small-sized species, with small home-ranges and specialized niches in spatially and temporally heterogeneous environments (e.g., [84]). Using ground-data from vole surveys across different periods, we are able to demonstrate that spectral, spatial, and temporal information from Sentinel-2 (Sentinel-2A MSI Level-1C) multispectral images analysis is significantly important to predict the Cabrera vole occurrence.
Results show that NDI45 vegetation index describing areas characterized with low-medium chlorophylls is the most important Sentinel-2-derived proxy for Cabrera vole habitat. High values of this index photosynthetically indicate higher biomass activity, i.e., dense canopies and crops linked to intensified agriculture practices, which are not suitable for the species. On the other hand, very low values of NDI45 indicate increasingly lower soil vegetation cover, which is also not suitable for the species occurrence. Reasons for a higher importance of this index during the 'Spring' should be related to increased wetness and mild temperature conditions during this period, which promotes annual grasses growth, ensuring higher vegetation cover, hence more available resources and improved habitat quality for Cabrera vole [13,14,52]. Multispectral satellite remote-sensing indices (e.g., NDVI) have been proven to successfully explain small mammal species distribution through the use of Landsat 7 [37] and Sentinel-2 data [42]. However, the present study showed that NDI45 is a better predictor than NDVI because it uses spectral information from the red-edge region, which has been recognized to provide more sensitive measurements of vegetation biophysical properties [60,85].
The SWIR1 spectral band obtained from the autumn season was ranked as the third most important variable in predicting Cabrera vole spatial distribution, and, during this season (in particular in 2018; see section 2.2), the grasslands over the study area were extremely dry due to the exceptional high temperatures and lack of rain. This is noteworthy because in the shortwave infrared region, the reflectance reduces as the amount of water content increases in vegetation [32] such that

Discussion
Results yielded evidence that fine-scale remote-sensing data may be useful to predict favorable habitats for the occurrence of small-sized species, with small home-ranges and specialized niches in spatially and temporally heterogeneous environments (e.g., [84]). Using ground-data from vole surveys across different periods, we are able to demonstrate that spectral, spatial, and temporal information from Sentinel-2 (Sentinel-2A MSI Level-1C) multispectral images analysis is significantly important to predict the Cabrera vole occurrence.
Results show that NDI45 vegetation index describing areas characterized with low-medium chlorophylls is the most important Sentinel-2-derived proxy for Cabrera vole habitat. High values of this index photosynthetically indicate higher biomass activity, i.e., dense canopies and crops linked to intensified agriculture practices, which are not suitable for the species. On the other hand, very low values of NDI45 indicate increasingly lower soil vegetation cover, which is also not suitable for the species occurrence. Reasons for a higher importance of this index during the 'Spring' should be related to increased wetness and mild temperature conditions during this period, which promotes annual grasses growth, ensuring higher vegetation cover, hence more available resources and improved habitat quality for Cabrera vole [13,14,52]. Multispectral satellite remote-sensing indices (e.g., NDVI) have been proven to successfully explain small mammal species distribution through the use of Landsat 7 [37] and Sentinel-2 data [42]. However, the present study showed that NDI45 is a better predictor than NDVI because it uses spectral information from the red-edge region, which has been recognized to provide more sensitive measurements of vegetation biophysical properties [60,85].
The SWIR1 spectral band obtained from the autumn season was ranked as the third most important variable in predicting Cabrera vole spatial distribution, and, during this season (in particular in 2018; see Section 2.2), the grasslands over the study area were extremely dry due to the exceptional high temperatures and lack of rain. This is noteworthy because in the shortwave infrared region, the reflectance reduces as the amount of water content increases in vegetation [32] such that SWIR1 can be sensitive to the existing senescent vegetation in the study area because it reaches a peak in terms of spectral reflectance [65]. Hence, it is reasonable to interpret grassy areas with some moisture conditions as associated with medium values of SWIR1. More specifically, a possible ecological explanation for the better support of SWIR1 during dryer periods is that the Cabrera vole might temporally respond to the leaf senescence spectral signals of perennial grasslands, which may help individuals' survival during most adverse environmental conditions (e.g., [30,52]).
Rao's Q metric is a measure of landscape beta diversity and can be a surrogate for landscape heterogeneity [70]. In the context of study area, the species occurs mainly in small marginal patches embedded in or surrounded by larger forest or agricultural areas, or on road verges [14]. The Rao's Q metric seems to be capturing this landscape diversity signal by showing that the species occurrence is favored in moderately heterogeneous landscapes. This pattern was particularly marked in Spring, when grasses become abundant, vegetation heterogeneity is higher, and vole populations increase given the higher availability of resources [12,16]. By contrast, a low suitability for homogeneous areas emerged from our analysis, suggesting vulnerability to habitat simplification, derived for instance from agricultural intensification or grazing pressure [86], which are known to have major impacts on small mammal habitat specialists [87] and for the Cabrera vole in particular [11]. Reasons for the slight decline in species probability of occurrence at the most heterogeneous areas (higher Rao's Q) are unclear, but may be related to the existence of shrubby areas where predation risk might be greater [11].
Apart from Sentinel-2, Cabrera vole occurrence probability peaks on close proximity to roads. This agrees with previous studies showing that the species often occurs on vegetated road verges, particularly in intensive agricultural or grazed areas [11,13,52]. This result does not necessarily suggest that the species is resilient to the negative effects that roads may exert on wildlife [88]. Instead, it emphasizes the compelling role of road verges in providing refuge habitats and corridors for small mammals, particularly where the surrounding matrix is mostly inhospitable [11,15,89,90]. Nevertheless, a major drawback of road verge habitats is that they may induce road-related mortality [91], which should be duly considered when the goal is to promote the use of verges as habitat and/or corridors for biodiversity.
Interestingly, along with the identification of suitable road verges, other semi-natural infrequently managed areas such as banks and field margins were identified in the habitat suitability model ( Figure 4). The conservation value of such areas is remarkable, as they usually support high levels of biodiversity, being key elements of High Nature Value farmland [92]. In addition, suitable areas for the Cabrera vole are often associated with Mediterranean temporary ponds [13], which are priority habitats under the EU Habitats Directive. Protecting such areas may be strategic for the conservation of the Cabrera vole, as well other species in human-dominated landscapes with limited availability of suitable habitats. Also, given the spatially limited and scattered distribution of those habitats, proper identification of priority conservation areas to ensure vole' populations viability, can potentially rely on landscape connectivity assessments (e.g., [93]). Once those areas are identified, conservation actions should consider the implementation of agri-environmental schemes, namely in the context of the European Union's Common Agricultural Policy, through which farmers are paid for restoring habitats, for instance by reducing the grazing pressure [11,15,94].
Earlier SDMs developed for Cabrera vole were carried out at broad scales and relied mostly on bioclimatic variables [10,24]. Despite the conservation value of macro ecological approaches for mapping environmental suitability at large scales [95], such models do not allow identifying, predicting, and mapping small key habitats [96], and thus are insufficient for defining concrete conservation actions. The use of fine-scale remote-sensing variables may thus provide a cost-effective tool to better support conservation planning with reduced survey costs [36], which may be crucial for rare and vulnerable species [97,98]. Higher mapping accuracy, especially when identifying grassland and linear land cover features, could be increased with images possessing very-high spectral and spatial resolutions, namely from data having a resolution spanning around 5m of detail, as suggested by Thornton et al. [99] and Rapinel et al. [100], possibly fulfilled through fusion of Sentinel 2 data [101]. Nevertheless, the use of very-high resolution data may be prohibitive for SDMs applications over larger areas due to its acquisitions costs. In this context, the use of Sentinel-2 data for habitat suitability mapping should be viewed as an effective compromise between spatial (10 m) and temporal resolution (5-6 days), as well as its open-data policy. Regarding the statistical methods inherent to SDMs, further studies are recommended in this research field in order to understand the best robustness of approaches able to handle high dimensional data [102], as well addressed to examine the predictive performances of multiple algorithms, especially when concomitantly integrated into an ensemble modeling framework [18,43]. This would be particularly interesting when evaluating how sub-sampled group of variables (remote-sensing products, topography, landscape variables) may singularly impact on the performance of species distribution models.
Our findings support the potential of remote sensing for mapping microhabitat suitability of rare small species, which until recently, was largely impracticable due to resource limitations [103]. Sentinel-2 is an open-access resource that provides spatial data at a resolution useful and necessary for this task, and, despite its relatively recent release, effective long-term ecosystem monitoring at local, regional, and national levels is planned to be continuously ensured by this satellite. As such, considering the increasing Sentinel-2 temporal span, future studies on conservation planning incorporating information for longer periods, as it is actually done with other satellites [104], may be valuable because they more likely minimize the common pitfall of assuming stable environmental suitability, and therefore populations persistence, over time [105,106].

Conclusions
Wildlife habitat selection is increasingly understood through the lens of earth observation remote-sensing instruments, either commercial or open-access. We demonstrated that the use of Sentinel-2-derived habitat variables, incorporating biophysical, spectral, and structural landscape information at fine-scales in different seasons, when integrated into RF machine learning methods, may support the identification of potential favorable areas for small and elusive species in dynamic landscapes. Overall, our study highlights that super-resolved remote-sensing data may provide an important tool for identifying linear habitat features (e.g., [99]). Sentinel-2 may provide high-quality and open-access data for fine-scale conservation planning and population monitoring, which may be particularly adequate when considering patchily distributed, small, rare, and elusive species. Finally, our study supports the view that the integration of detailed remote-sensing data into species distribution models is the next stage for linking species occurrences to environmental conditions at functionally relevant spatio-temporal scales, which is a central issue in ecology and conservation.