The Role of Remote Sensing Data in Habitat Suitability and Connectivity Modeling: Insights from the Cantabrian Brown Bear

: Ecological modeling requires sufﬁcient spatial resolution and a careful selection of environmental variables to achieve good predictive performance. Although national and international administrations offer ﬁne-scale environmental data, they usually have limited spatial coverage (country or continent). Alternatively, optical and radar satellite imagery is available with high resolutions, global coverage and frequent revisit intervals. Here, we compared the performance of ecological models trained with free satellite data with models ﬁtted using regionally restricted spatial datasets. We developed brown bear habitat suitability and connectivity models from three datasets with different spatial coverage and accessibility. These datasets comprised (1) a Sentinel-1 and 2 land cover map (global coverage); (2) pan-European vegetation and land cover layers (continental coverage); and (3) LiDAR data and the Forest Map of Spain (national coverage). Results show that Sentinel imagery and pan-European datasets are powerful sources to estimate vegetation variables for habitat and connectivity modeling. However, Sentinel data could be limited for understanding precise habitat–species associations if the derived discrete variables do not distinguish a wide range of vegetation types. Therefore, more effort should be taken to improving the thematic resolution of satellite-derived vegetation variables. Our ﬁndings support the application of ecological modeling worldwide and can help select spatial datasets according to their coverage and resolution for habitat suitability and connectivity modeling. A.G., T.G., M.C.M.-S., J.I.G.-V., M.M. and A.M.; visualization, P.C.-A., T.G., A.G., M.C.M.-S., J.I.G.-V. and M.M.; A.G., J.I.G.-V. and M.C.M.-S.; A.G., J.I.G.-V. and M.M.;


Introduction
Biodiversity loss is among the biggest challenges that nature conservation is currently facing. As a result, many conservation efforts are being continuously made by scientists, government administrations and conservationists. Computer-based approaches are acknowledged as one of the most useful and unbiased instruments to guide the design of nature management plans and policies, with habitat suitability and connectivity models being particularly helpful to produce spatially explicit information for conservation [1][2][3][4][5][6]. Nevertheless, sufficient spatial resolution of data and a careful selection of environmental variables are essential to ensure reliable model performance [7]. Accordingly, considerable effort has been taken to increase spatial resolution. The use of remote sensing imagery is gaining relevance because of its potential high spatial resolution and its broad temporal and spatial coverage and availability [8][9][10].
which includes the entire brown bear population of the Cantabrian Range. The mountain-ous zones have kept a considerable section of the study area with its original natural cover despite the human influence found in the region. The anthropogenic activity created a heterogeneous landscape with a mosaic of land covers where artificial and agricultural surfaces coexist with a diversity of grasslands, shrublands and forests. Natural areas are mainly composed of oak (Quercus sp.), beech (Fagus sylvatica) and chestnut (Castanea sativa) forests, heathlands (Erica sp.), shrublands (Cytisus sp.) and grasslands.

Study Area
The study area is located in the Cantabrian Range (NW Spain) and covers 35,700 km 2 , which includes the entire brown bear population of the Cantabrian Range. The mountainous zones have kept a considerable section of the study area with its original natural cover despite the human influence found in the region. The anthropogenic activity created a heterogeneous landscape with a mosaic of land covers where artificial and agricultural surfaces coexist with a diversity of grasslands, shrublands and forests. Natural areas are Remote Sens. 2021, 13, 1138 4 of 22 mainly composed of oak (Quercus sp.), beech (Fagus sylvatica) and chestnut (Castanea sativa) forests, heathlands (Erica sp.), shrublands (Cytisus sp.) and grasslands.

Species Occurrence Data
We used a comprehensive dataset of brown bear locations of the eastern and western subpopulations gathered from 2000 to 2010 by professional observers and rangers through direct observations and indirect evidence including excrements, footprints, beehive attacks and fur (details in [32]). This set of locations represents a time when both subpopulations were largely isolated and connectivity was identified as a major conservation concern for the species. This situation provides a strong opportunity to assess potential differences with far-reaching consequences in the definition of corridors and habitat areas where successful implementation of measures is fundamental. We unified the presence records and resampled them to 1 hectare pixel size to enhance the processing time. In total, we obtained 6207 brown bear locations to fit the habitat suitability models.

Environmental Variables and Spatial Datasets
We considered a set of six environmental variables to represent foraging resources, shelter (expressed by forest area and terrain ruggedness) and human pressure (expressed by building, highway and road density) to model the habitat suitability of the brown bear ( Figure 1, Table 1). Previous studies have shown this set of environmental variables to effectively predict brown bear presences in the Cantabrian Range [16,32]. Foraging resources and forest area variables were obtained from (1) a Sentinel-1 and 2-based land cover map of 2016-2017 (SLCM) [27], (2) the Forest Type Product 2015 [25] of CLMS, (3) the CORINE Land Cover 2018 map [26], (4) the Forest Map of Spain [29] (FMS;1998 and (5) LiDAR data (2009-2012) from the Spanish National Plan for Aerial Orthophotography [28]. See Appendix A for more details of these datasets. Terrain ruggedness and human pressure variables were constant in all the models and derived from the same data sources.  [33] We combined the spatial data from the previous sources into three datasets according to their spatial coverage-i.e., global, continental or national-to produce the environmental variables used as predictors of habitat suitability. However, the extent of the analysis was constant (i.e., the Cantabrian Range) and the terms "global", "continental" and "national" were only used to indicate dataset characteristics regarding spatial coverage, data availability and thematic, spatial and temporal resolution (see Appendix A for dataset details). Therefore, for the global level, we used the SLCM alone. Meanwhile, we combined CLC18 and FTY for the continental level, and LiDAR with FMS for the national level.
We resampled the environmental variables to 1 hectare pixel size to match species occurrence data resolution. We considered nine different spatial scales for each variable, each scale representing a particular brown bear's point of view and reaction to the environment. This environmental gradient is of particular relevance for its ecological requirements, as the brown bear may behave differently when selecting resources from patch to landscape scales [34]. We calculated mean values using circular moving windows with nine different radii (0. 25, 0.5, 1, 2, 4, 8, 16, 32 and 64 km). We standardized all variables by subtracting the mean and dividing by the standard deviation.

Foraging Resources
We calculated foraging resources as the abundance of plant species (trees, shrubs and herbs) weighted by the importance of each species in the diet of the brown bear (for details, see [32]). Importance was obtained from a previous report on the brown bear diet based on scat surveys [35]. Abundance was calculated in different ways depending on available information of each spatial dataset. For the national dataset, we used the information available in the FMS map, which includes tree species identities and total tree canopy cover and percentage of cover for each tree species. However, we replaced the canopy cover from FMS by the one derived from LiDAR data, which is more detailed due to its finest spatial resolution (25 m; Appendix A). We estimated the non-tree species abundance combining information from available floristic inventories, ecological niche models fitted through penalized logistic regression and based on climatic and lithological predictors (Appendix C) [36], and expert knowledge on the potential presence of specific plant species within plant communities identified by FMS. We used plant species that are basic for the brown bear diet (Appendix C). For the continental and global datasets, we used the abundances estimated in [35], adapting them to fit CLC18 and SCLM classes by averaging the original abundance values for each class. Additionally, for the continental dataset, we replaced the forest covers of CLC18 with the ones of the FTY due to its finer grain size.

Shelter
We used forest area and terrain ruggedness as proxies for brown bear shelter following the methods described in [16,32]. The definition of forest area varied slightly among datasets. For the global dataset, we obtained it from the forest class of the SLCM. We used the FTY layer for the continental dataset. For the national dataset, we applied a threshold of 10% of canopy cover of LiDAR data to determine forest pixels, which is the same threshold applied to the FTY (Appendix A; [25]). We calculated terrain ruggedness from a 25-m Digital Elevation Model [28], according to [37].

Human Pressure
To study human pressure, we used three variables that represent buildings and transport infrastructures: the density of buildings, highways and conventional roads. We obtained buildings data from a topographic map at 1:25.000 scale (MTN25 2015) [28] and transport infrastructures from Open Street Maps 2015 (downloaded from [33]). Transport infrastructures were divided into conventional roads and highways due to their different traffic volumes and physical limitations for the brown bear distribution. Both vector datasets were rasterized with 100 m pixel size and were then processed in the same way as other variables for the multi-scale analysis (see Section 2.3).

Habitat Suitability Models
First, we developed univariate models for every combination of environmental variable, scale (radii of the circular moving window; see Section 2.3) and dataset (global, continental and national) to find the operational scale, which is the scale that correlates best with the occurrence data ( Figures A1-A3) [32,38,39]. For this purpose, we applied lrm and pentrace functions of the rms R package [40] to fit penalized logistic regressions with data of each scale as a unique predictor and assuming a monotonic response. The outcome binary variable was built merging bear occurrences with 20,000 pseudo-absences generated randomly as background cells inside the study area. We selected the model with the best value of Akaike's information criterion (AIC; [41]). The scale of the best performing model was identified as the operational scale (i.e., the scale that best represents the habitat selection pattern of the species for each variable).
Second, we used every environmental variable at its operational scale to fit multivariate and multi-scale habitat suitability models for each of the three datasets. We calculated penalized logistic regression models using linear terms and no interactions among predictors, using lrm and prentrace functions of the rms package. To study the predictive ability of each of the three multivariate models, we relied on the area under the receiver operating characteristic curve (AUC; [42]), tested with ten-fold cross-validation. Next, we tested the difference between the AUC of the models by applying DeLong's test [43]. The predicted suitability values of each multi-scale model were mapped to obtain three habitat suitability maps (one for each dataset).

Connectivity Analysis
To assess brown bear connectivity, we used an individual-based approach and the graph theory, which is one of the most common methods to study ecological connectivity [44][45][46]. Graphs are mathematical representations of the landscape that allow effectively studying complex connectivity scenarios [47,48]. Graphs are composed of nodes and links between nodes, which ultimately represent habitat patches and movement corridors, respectively. The links are preferably characterized by effective distances which consider the species cost of movement through the landscape matrix [49]. Thus, higher distances indicate a high probability of isolation between source and destination habitat areas. To compute effective distances, it is therefore necessary to have a resistance surface (i.e., a raster layer) whose pixel values express the cost of movement, usually derived from the land covers present in each cell [50,51]. We derived resistance surfaces from each of the three bear habitat suitability maps by assuming an inverse relation of movement resistance with habitat suitability using the following expression: where R is the resistance value calculated for a given pixel, p is the value of habitat suitability in that pixel and max(p) is the maximum value of habitat suitability across all pixels in the study area. The transformation of habitat suitability values resulted in three resistance surfaces with values ranging from 1 to 100. The lowest values represent pixels of high habitat suitability and low resistance and vice versa. For each of the three resistance surfaces, we computed distances between nodes in two ways according to two of the three connectivity methodologies of this study. Nodes consisted of an array of 316 pixels resulting from integrating the brown bear presence records into a grid with 5 km cells (a size previously used for the species [38,52]) to reduce spatial autocorrelation and avoid overlaps.
We used three methodologies to assess connectivity, the first two to characterize links between nodes and a third one to determine the links' importance for connectivity. First, we used least-cost modeling to characterize links as the paths that imply the minimum effective distance between two nodes (i.e., least-cost paths), which are more likely to be used by animals due to lower dispersal efforts and mortality risk [53]. We used Linkage Mapper software to compute least-cost paths and their effective distances [54].
Second, we used circuit theory [55] to characterize links. Circuit theory allows the calculation of currents. Currents represent net probabilities of movement which are estimated from random walks between nodes and resistance surfaces. Corridors are characterized as a set of cells with high probabilities of being used-i.e., high current flows. As in least-cost modeling, circuit theory measures effective distances but incorporates multiple pathways connecting nodes (i.e., resistance distances). Therefore, several available dispersal routes among locations can be identified including optimal and suboptimal ones. We applied circuit theory with Circuitscape v4.0.5 [56].
For the third methodology, we used the probability of connectivity index (PC). This index is based on the habitat availability concept, which assumes that connectivity not only happens between habitat patches, but also inside them [48]. The probability of successful connection between two nodes is calculated as a decreasing exponential function of the effective distances between them. In our case, the effective distances between nodes were calculated through least-cost modeling. Then, these values are summed up to consider all the connections present in the landscape (for details, see [49]). Therefore, it is possible to calculate the contribution of each link to overall connectivity by calculating the difference in PC (dPC) when each link is removed, those links with a higher dPC being more important [57]. To calculate dPC values, we used the Link Removal function of Conefor software [58]. To estimate the probabilities, we used the median dispersal distance of the species, which represents a 0.5 probability of connection between two nodes separated at that distance. We considered a dispersal distance of 119 km for the brown bear [59], which we multiplied by the mean value of each resistance surface to obtain the effective dispersal distances.
We compared the results derived from the three resistance surfaces (global, continental and national) when applying the three methodologies (least-cost modeling, circuit theory and PC index). First, we obtained connectivity maps developed from each methodology and visually compared the distribution and location of movement pathways. Second, we compared the correlation of the connectivity metrics obtained from the three methodologies (effective and resistance distances and dPC values) between resistance surfaces. Last, we applied the three methodologies to a genetic resistance surface (genetic model) and compared the results to those obtained from the three previous resistance surfaces. The genetic resistance surface was obtained from brown bear genetic samples (17 polymorphic microsatellite loci) integrated with landscape variables following a genetic-multiplicative approach (for details, see [51]).
We analyzed the correlations of the three connectivity metrics for three different scenarios to characterize movements with strategic significance in terms of intra-and intersubpopulation movements. The three scenarios include movements within the eastern (EE) and western (WW) subpopulations and between subpopulation edges (WE)-i.e., border of each subpopulation closest to the other. We considered only the edges and not subpopulation centers to make a fair comparison between methodologies since Linkage Mapper and Conefor process the links of spatially adjacent nodes and only considered movements between subpopulation edges. Another reason is that movements inside subpopulations are already considered in two scenarios and happen in suitable habitat areas where these movements are feasible and frequent. In addition, connectivity conservation efforts are mostly focused on the low permeability gap between subpopulations [38,51].

Operational Scales and Predictors Effect
The operational scales of the vegetation variables selected by the univariate models varied among models. Foraging resources were similar between global (0.5 km) and continental (2 km) datasets but differed significantly from the national one (16 km), which was much broader. Forest area, on the contrary, was analogous for the three datasetsvariation between 0.5 and 1 km. Terrain ruggedness influenced at a relatively small scale of 2 km and all human pressure variables best performed at 16 km in the three models ( Table 2). Vegetation variables had a positive effect ( Table 2) on habitat suitability in most of the cases, especially the foraging resources of the national dataset. Together with foraging resources, terrain ruggedness was the environmental variable that most fostered bears' habitat suitability in every multivariate model. Forest area had a negative influence on the global model, while it was slightly positive in the continental and national models. In all cases, the relation of human pressure variables was unfavorable.

Model Performance and Habitat Suitability Spatial Pattern
All the models yielded a strong predictive performance with high AUC scores ( Table 2). The national model achieved the highest AUC (0.936) while the lowest belonged to the global model (0.900). Results of DeLong's test showed a significant difference between the three models (p-value < 0.05). Overall, the spatial pattern of the areas with high habitat suitability was very similar and almost identical between the global and the continental models, although with differences in their sharpness. High values appeared in two cores, which matches the brown bear distribution range in the study area. Both cores were separated by low-suitability areas where two main highways occur (AP-66 and N-630; see Figure 2). High suitability also appeared towards the southwestern and southeastern borders of the study area in the global and continental model predictions.

Corridor Comparison
Results from the connectivity analysis were variable according to the methodology used to compare the different datasets. Links obtained through least-cost modeling showed mostly similarities between datasets and the genetic model, especially inside both subpopulations (see Figures 3 and A4 in Appendix B for the genetic model). For the connections between subpopulations, the links derived from three datasets showed few differences as they mostly concurred, and when they did not, their routes were often parallel. Nevertheless, their differences were more evident for longer links in the east border of the eastern subpopulation and in the northwest.
The three datasets revealed moderately different current flow patterns when using circuit theory (see Figure 4). The global model drew a sharper and more scattered pattern of corridors, in accordance with its corresponding habitat suitability map. On the contrary, the spatial pattern of the continental model was more homogeneous and concentrated in fewer corridor areas. The national model showed a less dendritic pattern than the global model but with less condensed flow values than the continental model. In the three models, high current flow values appeared in similar areas but with more irregular patterns within the eastern subpopulation. Current flow values were low between subpopulations in the three models.

Corridor Comparison
Results from the connectivity analysis were variable according to the methodology showed mostly similarities between datasets and the genetic model, especially inside both subpopulations (see Error! Reference source not found. and Figure A4 in Appendix B for the genetic model). For the connections between subpopulations, the links derived from three datasets showed few differences as they mostly concurred, and when they did not, their routes were often parallel. Nevertheless, their differences were more evident for longer links in the east border of the eastern subpopulation and in the northwest. The three datasets revealed moderately different current flow patterns when using circuit theory (see Error! Reference source not found.). The global model drew a sharper and more scattered pattern of corridors, in accordance with its corresponding habitat suitability map. On the contrary, the spatial pattern of the continental model was more homogeneous and concentrated in fewer corridor areas. The national model showed a less dendritic pattern than the global model but with less condensed flow values than the continental model. In the three models, high current flow values appeared in similar areas but with more irregular patterns within the eastern subpopulation. Current flow values were low between subpopulations in the three models.

Comparison of Connectivity Metrics
We found high correlations when comparing connectivity metrics among the three resistance surfaces (see Table 3). Under least-cost modeling, Pearson's correlation estimates (r) of effective distances were remarkably strong, with coefficients higher than r > 0.95 between global, continental and national models. The highest correlation recorded was inside the western subpopulation between the global and continental models (r = 0.996). Strong correlations remained when comparing the effective distances of the three models with the genetic model (r > 0.9).
In comparison to least-cost modeling, the correlation patterns were similar when using circuit theory, although weaker correlations were observed. The best correlations remained between the global and continental models for every type of movement. The correlation coefficients were moderate (0.8 < r < 0.5) between each of the three models and the genetic connectivity model. The highest coefficients were recorded by the national model, especially for movements between subpopulation edges.
The dPC values showed less consistent correlations among models (see Table 3). Continental and national models showed moderate correlations. However, the correlations of the global model were weaker, with poor correlations between subpopulation edges. The global and continental models barely correlated with the genetic model. Conversely, the national and the genetic dPC values correlated moderately for movements between subpopulations.

Comparison of Connectivity Metrics
We found high correlations when comparing connectivity metrics among the three resistance surfaces (see Error! Reference source not found.). Under least-cost modeling, Pearson's correlation estimates (r) of effective distances were remarkably strong, with coefficients higher than r > 0.95 between global, continental and national models. The high-  Figure A5 in Appendix B for the current map of the genetic model). Table 3. Pearson's correlation of connectivity metrics: effective (least-cost modeling) and resistance (circuit theory) distances, and difference in probability of connectivity index (dPC) values (PC index) between the paired combinations of global, continental and national models and of each of these three models against the genetic model for each subpopulation (subp.) and between subpopulation edges.

Discussion
The choice of spatial data source is an important step in ecological modeling since the resolution of data is a key factor for landscape definition that deeply influences the nature of the results [60,61]. Accordingly, our results show similar outcomes but also meaningful divergences when modeling habitat suitability and connectivity using different data sources.

Habitat Suitability Models
Our results show similar predictive performances for the three habitat suitability models, which suggest that Sentinel-1 and 2 imagery and CLMS pan-European products are solid alternatives for modeling habitat suitability. Besides minor differences in AUC scores, habitat suitability maps showed analogous spatial patterns of high-suitable areas with the three models correctly predicting the distribution range of the bear in the Cantabrian Range.
The fine-spatial resolution data included in every model-i.e., Sentinel-1 and 2, FTY and LiDAR-may be the main cause behind the good performances, since finer grains make land cover classes more distinguishable for the model [62,63]. The high spatial resolution led to the selection of small operational scales (0.5-1 km) in the three models for the forest area variable and allowed the identification of isolated high-suitability areas (forest patches with abundance of food resources) in such a heterogeneous landscape [16]. However, the delineation of high-suitability areas did not have a detailed pattern in the continental and national maps. This could be due to CLC18 and FMS as their broader scale (MMU of 25 and 2.25 ha, respectively) may have affected the detail of the results.
To match the resolution of the bear presence data, we resampled the resolution of the environmental variables to 1 ha, coarsening the fine pixel size of some layers (e.g., from LiDAR, FTY and SLCM). Brown bear locations had to be resampled to 1 ha to avoid pseudoreplication since multiple traces could have been caused by the same individual.
GPS and tracking data could improve the accuracy of the locations, allowing using the abovementioned environmental variables with their best spatial resolution. This could give new insights into the influence of high spatial resolution on habitat suitability modeling. However, for most cases, brown bear habitat preferences are more influenced by coarser operational scales and not by the finest scale available of 250 m ( Table 2). Despite the operational scales, the high spatial resolution of the layers did influence how the configuration of the landscape was characterized, considering small habitat patches and landscape elements that would have been ignored with coarser spatial resolutions before resampling.
Thematic resolution may explain the slightly different predictive performances among models, as performance increased when their datasets classified more vegetation types [12,62]. For instance, the identification of multiple vegetation types included in the national model (due to the FMS) allowed characterizing important foraging resources for the brown bear-in comparison to the SLCM which only classified one class of forest. The richer thematic resolution of the national dataset may also explain the different operational scales selected for the foraging resources variable.
The inaccuracies of input layers could also produce lower performances in the continental and global models. Both the FTY (minimum overall accuracy of 0.90) [25] and the SLCM layers (kappa index of 0.89) [27] had high accuracies. However, some misclassifications are also visible, especially in the SCLM. The steep slopes negatively affected Sentinel-1 data, producing classification errors in mountainous areas where vital natural resources for the brown bear are abundant. More misclassifications appeared between forested areas and shrublands, both important for the definition of foraging resources [27]. These inaccuracies could explain the negative influence of the forest area variable in the global model (Table 2) as some relevant suitable areas for the bear were not correctly represented as forest in this model. Sentinel-1 allows mapping vegetation variables regardless of cloud conditions that limit optical imagery [64], but the potential issues derived from its use in hilly landscapes should be addressed [65].
Our findings support the capabilities of satellite imagery to predict species habitat preferences in heterogeneous landscapes [13,66]. However, this may not be appropriate for all circumstances. The considered global and continental datasets may be more suitable for cases with broad ecological requirements, as these may not attain enough detail in identifying a wide variety of natural resources, which can be important for understanding species' habitat requirements at small scales [32,67,68]. When trying to understand specific habitat-species relations, spatial analysts should use data that enrich the thematic resolution [69] and that complement the high spatial resolution of the input vegetation variables. For this purpose, free remotely sensed data sources, such as the recent space-borne NASA LiDAR, can be useful to inform about the structure of vegetation [22].

Connectivity Models
The similar trajectories, current flow patterns and high correlations of distances among models when applying least-cost and circuit theory methodologies indicate that varying data sources did not profoundly affect connectivity results. However, even the slight variations recorded in the results could have considerable implications for conservation [51,63]. For instance, lower correlations in resistance distances and slightly different corridor patterns suggest that the data source choice can have implications for corridor identification through circuit theory. Much of the difference among current maps may be derived from the effect of low-spatial resolution data (CLC18, FMS; see Section 4.1). Furthermore, the moderate to low correlations of dPC values (Table 3) indicate different contributions to connectivity for analogous links. Thus, each connectivity model would prioritize distinct corridors, which could have repercussions in the landscape planning of connectivity measures. In any case, predicted corridors at regional scales should be used as a first step in characterizing connectivity networks and in combination with studies conducted at different scales. More detailed data suitable for local scales should be used to ascertain if the predicted corridors are functional and better define the routes, habitat patches and landscape features before implementing specific measures of conservation projects [70,71].
Grain size can slightly but importantly affect the total cost of movement of a corridor [50,72], as finer grains capture landscape and vegetation patterns that should be represented to not misestimate corridor trajectories [73]. Considering the similar spatial grain of the vegetation variables, the differences among connectivity model outcomes may be due to variable thematic resolutions for estimating their resistance surfaces [72].
Genetic connectivity models are more representative of dispersal movements and successful reproduction between subpopulations, making them more appropriate for corridor planning [50,51,74]. Thus, the strong correlations (Table 3) with the genetic model when using least-cost modeling suggest that the spatial datasets considered in this study can offer robust connectivity results under this approach. It should be noted that leastcost modeling is deterministic and predicts only optimal paths [51,55,75]. Despite having similar effective distances for analogous paths, minor changes in their trajectories could focus the implementation of conservation resources on different places. This effect was more prominent as the Euclidean distance between nodes increased (Figure 3). Under circuit theory, the correlations with the genetic model were moderate but the strong coefficient between subpopulation edges of the national model highlights this model as the best option to predict corridors through degraded landscapes. A similar result was obtained for dPC values, supporting the national model also as the best alternative to assess the corridors' importance between habitat areas. These results are particularly relevant, considering that conservation efforts often focus on the identification of corridors between disconnected habitat patches to enhance gene flow, species dispersal and adaptation to current threats [51].

Other Considerations on the Choice of Spatial Data Source
Above, we discussed the effects that the choice of spatial data source can have on two types of ecological models (i.e., habitat suitability and connectivity). If multiple data sources are available, our results advocate for the data source with the highest spatial and thematic resolution. However, the regional availability of data can limit the choice of data source, which in most cases would be restricted to satellite data, such as Sentinel-1 and 2 imagery. Nevertheless, the temporal resolution of data can also be essential, especially for dynamic habitat and connectivity studies [76]. Satellite remote sensing vegetation variables offer a great advantage as these can be updated regularly with a high frequency-5 to 6 days of average revisit time intervals for Sentinel sensors at the equator [20,21]. Despite being derived from satellite data, CLMS products do not offer frequent updates-every 6 years for Corine Land Cover and every 3 years for the FTY. This period can be even longer for nationally restricted data; e.g., Spanish LiDAR data are updated approximately every 7 years, while the FMS has not yet been updated. Consequently, the temporal resolution would favor the choice of satellite imagery to model ecological processes when the time factor is important. Additionally, the combination of satellite data with other data sources can be recommendable to offset the caveats that these sensors can have [22]. Nonetheless, satellite imagery processing requires great expertise and effort to map environmental variables. This limitation supports the choice of already processed cartography, such as CLCS products, which allow conservationists to focus almost directly on the modeling process [77].

Conclusions
Our results support free satellite imagery as a global and powerful source to estimate vegetation variables for habitat suitability and connectivity modeling. Already processed land cover products (such as the ones from CLMS) are also an effective alternative for studies focusing on the European continent. Therefore, widely available data sources are suitable to produce meaningful results that guide the design of consequent conservation measures and nature management plans. However, some considerations depending on the required spatial, temporal and thematic resolution need to be taken into account. Globally available satellite data have the advantage of being updated frequently, but they may require sufficient processing skill and can be time-demanding and computationally intensive. Satellite data may be limited in many cases to understand precise habitat-species associations, which requires high thematic resolution besides fine spatial detail. Failing to map a sufficient diversity of natural resources for a species can have detrimental effects in some cases, such as in the identification and prioritization of corridors. For this reason, more effort should be applied to improve the thematic resolution of vegetation variables coming from broadly accessible sources, which can be achieved by combining sensors and other open source spatial data sources. The comparison presented by this study offers valuable insights to better understand the potential implications, biases and limitations that spatial data source selection has on the results of habitat suitability and connectivity modeling, which is key for the correct definition of conservation measures.

Data Availability Statement:
The data that support the findings of this study are available from Junta de Castilla y León, Gobierno de Cantabria, Principado de Asturias and Xunta de Galicia but restrictions apply to the availability of these data, which were used under license for the current study and are not publicly available. Data are, however, available upon request to Junta de Castilla y León, Gobierno de Cantabria, Principado de Asturias and Xunta de Galicia. The habitat suitability and connectivity maps are available from the authors upon request.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Spatial Datasets Details
Global dataset: • Sentinel Land Cover Map (SLCM): We used the land cover map of the study area produced by [27] (see for more details) using Sentinel satellite imagery of the Copernicus Programme of the European Space Agency (ESA). Specifically, the imagery used includes SAR data from Sentinel-1A Ground Range Detected (GRD) products (L1C) using VV and VH polarizations, together with optical data from Sentinel-2 (L2A) using ten spectral bands at 10 and 20 m of spatial resolution. • Forest Map of Spain (FMS) from the Spanish Ministry of Ecological Transition [29]: FMS is detailed cartography of the composition and structure of the forest stands in Spain at a 1:50,000 scale with an MMU of 2.25 ha. It was produced between 1998 and 2007 with photo interpreted aerial imagery, combined with pre-existing maps and field inventory data. FMS is only available for Spain, and it is free of cost. • LiDAR data obtained from the Spanish National Plan for Aerial Orthophotography between 2009 and 2012 [28]: This dataset has a mean density of 0.5 points/m 2 and a vertical root mean square error ≤0.15 m. It was processed with FUSION software [78] and then aggregated using 25 m of spatial resolution. We calculated forest canopy cover as the ratio between the number of first returns above 3.5 m-to filter out understory vegetation-and the total number of first returns (details in [32]). LiDAR data are only available for Spain, and they are free of cost.

Appendix B. Supplementary Figures
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 22 and then aggregated using 25 m of spatial resolution. We calculated forest canopy cover as the ratio between the number of first returns above 3.5 m-to filter out understory vegetation-and the total number of first returns (details in [32]). LiDAR data are only available for Spain, and they are free of cost. Figure A1. Operational scales of the foraging resources variable for the three datasets. Figure A1. Operational scales of the foraging resources variable for the three datasets.

Appendix B. Supplementary Figures
Remote Sens. 2021, 13, x FOR PEER REVIEW 18 of 22 Figure A2. Operational scales of the forest area variable for the three datasets. Figure A2. Operational scales of the forest area variable for the three datasets. Figure A2. Operational scales of the forest area variable for the three datasets. Figure A3. Operational scales of the ruggedness, buildings, highways and roads variable, which are common for the three datasets. Figure A3. Operational scales of the ruggedness, buildings, highways and roads variable, which are common for the three datasets.
Remote Sens. 2021, 13, x FOR PEER REVIEW 19 of 22 Figure A4. Links based on the least-cost modeling methodology for the genetic-multiplicative resistance surface. Figure A4. Links based on the least-cost modeling methodology for the genetic-multiplicative resistance surface. Figure A4. Links based on the least-cost modeling methodology for the genetic-multiplicative resistance surface. Figure A5. Current map based on the circuit theory methodology for the genetic-multiplicative resistance surface.

Appendix C. Further Details of the Datasets Used for the Estimation of the Foraging Resources Variable
For non-tree species, we used previously available habitat suitability model predictions for woody plant species [36]. Habitat suitability values were converted to abundance estimates using logistic regression models trained with habitat suitability as the predictor variable and observed species abundance from vegetation surveys as the response. Vegetation survey data were obtained from [79]. Figure A5. Current map based on the circuit theory methodology for the genetic-multiplicative resistance surface.

Appendix C. Further Details of the Datasets Used for the Estimation of the Foraging Resources Variable
For non-tree species, we used previously available habitat suitability model predictions for woody plant species [36]. Habitat suitability values were converted to abundance estimates using logistic regression models trained with habitat suitability as the predictor variable and observed species abundance from vegetation surveys as the response. Vegetation survey data were obtained from [79].