Integrating landscape resistance and multi-scale predictor of habitat selection for amphibian distribution modelling at large scale

Species distribution modelling is a common tool in conservation biology but two main criticisms remain: (1) the use of simplistic variables that do not account for species movements and/or connectivity and (2) poor consideration of multi-scale processes driving species distributions. We aimed to determine if including multi-scale and fine-scale movement processes in SDM predictors would improve accuracy of SDM for low-mobility amphibian species compared with species-level analysis. We tested and compared different SDMs for nine amphibian species with four different sets of predictors: (1) simple distance-based predictors; (2) single-scale compositional predictors; (3) multi-scale compositional predictors with a priori selection of scale based on knowledge of species mobility and scale-of-effect; and (4) multi-scale compositional predictors calculated using a friction-based functional grain to account for resource accessibility with landscape resistance to movement. Using friction-based functional grain predictors produced slight to moderate improvements of SDM performance at large scale. The multi-scale approach, with a priori scale selection, led to ambiguous results depending on the species studied, in particular for generalist species. We underline the potential of using a friction-based functional grain to improve SDM predictions for species-level analysis.


Introduction
Species distribution models (SDMs) are currently the most widely used models for making spatially explicit predictions of the suitable environment of a species in a given area (Araújo et al. 2019;Zurell et al. 2020) based on the relationships between species observations and environmental variables (Guisan and Zimmermann 2000;Franklin 2010). Careful selection of variables can improve model performance and limit multi-collinearity and over-fitting problems (Araújo and Guisan 2006;Austin and Van Niel 2011;Dormann et al. 2013;Fournier et al. 2017). It is, therefore, necessary to rigorously select predictors based on species ecology, geographical context and ecological processes, at the appropriate spatial extent and resolution (Guisan and Thuiller 2005;Austin and Van Niel 2011;Araújo et al. 2019). However, problems of variable selection (Tulloch et al. 2016;Fournier et al. 2017) underly two major criticisms of SDMs: (1) the use of simplistic variables that do not account for fundamental ecological processes such as biotic interactions and species movements or connectivity (Franklin 2010; Barve et al. 2011;Boulangeat et al. 2012;Decout et al. 2012;Foltête et al. 2012;Holloway and Miller 2017) and (2) a lack of consideration of the multi-scale processes driving species distributions (Mayor et al. 2009;McGarigal et al. 2016;Fournier et al. 2017;Hallman and Robinson 2020). Indeed, most SDM studies use simple distance layers and/or map all predictors at the same grain (Fournier et al. 2017), whereas the interactions between species and landscape operate at a variety of spatial and temporal scales (Wiens 1989;Levin 1992;Mayor et al. 2009;Jackson and Fahrig 2012;McGarigal et al. 2016;Miguet et al. 2016).
Variable scale-effect, or the scale at which species perceive and interact with the landscape (Miguet et al. 2016), is often related to species mobility and movement processes (Mayor et al. 2009;Fahrig 2012, 2015;Thornton and Fletcher 2013). Movements are essential for animal survival, especially in fragmented human-dominated landscapes (Wiens et al. 1993;Hanski 1998;Sinsch 2014). They can be defined as regular and local movements (e.g. daily foraging movement); migration (i.e. seasonal movement inside the species annual home range) or dispersion (e.g. larger scale occasional unidirectional movement) (Pittman et al. 2014;Sinsch 2014;Holloway and Miller 2017). These movements can involve different levels of biological organization, from the individual moving from a daily roosting site to a foraging site within its home range, to the metapopulation characterized by individuals dispersing between more or less isolated sub-populations (Sinsch 2014). Thus, movements influence distribution patterns through processes acting at different scales and hierarchical levels, such as resource accessibility, colonization of new suitable habitats, gene flow between populations and recruitment of new individuals (Levin 1992;Hanski 1998;Mayor et al. 2009;Pittman et al. 2014). It is therefore essential to consider movements in SDM (Soberon and Peterson 2005;Richard and Armstrong 2010;Boulangeat et al. 2012;Brotons et al. 2012) though they are not always included (Foltête et al. 2012;Holloway and Miller 2017).
Several approaches can be used to integrate movements and landscape connectivity at different scales into habitat suitability studies. Some examples are (1) including structural and/or functional connectivity predictors by using simple metrics of landscape structure or more complex metrics derived from circuit or graph theory (e.g. Foltête et al. 2012); (2) applying dispersal kernel layers (Summers et al. 2012); (3) using complex hybrid approaches (i.e. correlative SDM combined with spatially-explicit population models; (Brotons et al. 2012); and (4) using SDM outputs to define nodes used in connectivity analysis (Estrada-Peña 2005;Duflot et al. 2018;Foltête et al. 2020) or to apply spatial filters to constrain species dispersal from previously occupied areas (Cardador et al. 2014). However, all such methods mainly concern broad-scale processes such as migration or dispersal movements at population level (Franklin 2010;Barve et al. 2011;Holloway and Miller 2017) while finer scale movements are neglected (Holloway and Miller 2017).
The use of functional grain including landscape resistance might be relevant to consider finer-scale movements (i.e. within grain movements) in broader scale SDM (Baguette and Van Dyck 2007;Romero et al. 2009;Galpern and Manseau 2013a). Functional grain is defined as the perceptual range (i.e. the minimum distance from which a resource may be detected), which is related to species mobility, with or without landscape resistance (Lima and Zollner 1996;Baguette and Van Dyck 2007;Galpern and Manseau 2013a). It is related to the scale at which an organism interacts with landscape structure (Romero et al. 2009) according to the ecological processes involved (e.g. foraging behavior, migration or dispersion). Functional grain therefore determines a set of functionally connected resource patches in a continuous area within which an animal may move freely according to its mobility (Galpern and Manseau 2013b). Awareness of functional grain may greatly enhance species distribution predictions (Romero et al. 2009) especially for species like amphibians, with low movement capacities and which are highly affected by landscape features, composition, and configuration (Cushman 2006;Collins and Fahrig 2017).
Movements through heterogeneous landscapes are crucial for amphibian survival (Pittman et al. 2014;Sinsch 2014). Most amphibians have a biphasic life cycle (i.e. with aquatic and terrestrial stages) which means they seasonally travel relatively long distances despite their low locomotion capacity and sensitivity to desiccation (Janin et al. 2012). The distance between aquatic and terrestrial sites inside their annual home range can exceed several hundred meters (Marsh and Trenham 2001;Smith and Green 2005) and dispersal movements can reach a few kilometres (Pittman et al. 2014;Sinsch 2014). Consequently, amphibians are particularly sensitive to landscape permeability defined by landscape structure and composition at different spatial scales (Cushman 2006;Denoël and Lehmann 2006;Youngquist and Boone 2014;Collins and Fahrig 2017). Modelling their distribution over large scales (i.e. a scale with an effect of both climate and landscape variables) is particularly challenging and under-investigated although of high relevance for regional conservation, especially in a context of global change (Holloway and Miller 2017;Préau et al. 2019;Matutini et al. 2021). Hence, amphibians are an interesting focal taxon for testing whether functional grain analyses that integrate finer-scale movement processes improve SDM accuracy over large scales.
The aims of this study-at large scale and using low mobility amphibians as a focal taxon-were to compare the performance of (1) distance-based versus landscape composition predictors of species distribution, (2) single-scale versus multi-scale models based on empirical knowledge of maximal mobility, and (3) multi-scale models using landscape composition within functional grains defined by fine-scale movements constrained or not by landscape resistance. We tested three associated predictions: (1) models using landscape composition predictors perform better than models using only distance layers as predictors; (2) a multi-scale approach based on prior knowledge of movement processes and the scale-ofeffect improves model accuracy, and (3) integrating fine-scale movements using a friction based functional grain improves prediction of amphibian habitat distributions at large scale.

Study species and area
Our study was performed in Pays de la Loire, a region of western France covering 32,082 km 2 where landscapes have been greatly modified by intensive agricultural activities over the last century. These have notably led to a decrease in forest cover, a homogenization of landscape mosaics and a degradation of wetlands. Traditional hedgerow network landscapes with small pastures delimited by hedgerows, ponds networks and small woods associated with extensive livestock farming have high, local conservation value (Baudry et al. 2000), especially for amphibians.
The 21 species of amphibians present in the region are protected and 12 of them are priority species for conservation (See Table 1).

Biological data for model calibration
We accessed high-resolution presence-only occurrences of nine species of amphibian (from 1500 to 9000 presence data per species) from a regional database for the period from 2013 to 2020 (see Matutini et al. (2021) and Online Appendix 2 for data description and selection procedure). After filtering out data points \ 500 m apart to reduce spatial autocorrelation, we projected remaining presence data on a 100 m resolution grid. This presence data selection process was repeated for each modelling iteration.

Environmental variables
The environmental data used to establish predictor layers over the study area are described in Matutini et al. (2021) and in Online Appendix 3. Each land cover type was rasterized at 10 m or 20 m resolution for metric calculations (see Table 2) with Chloe Software 4.0 (Boussard and Baudry 2017). We tested four different models performed with different sets of predictors (see Fig. 1): (1) Distance-based predictors calculated as the Euclidean distances from the grid-cell centroid to the nearest patch of a given land cover type.
Predictors included in DIST model; (2) Composition-based and single-scale predictors calculated with a moving window analysis in a circular buffer around pixel centroids based on species migration ability, defining potential annual home range size. Predictors included in CIRC.SS model; (3) Composition-based and multi-scale predictors calculated with a moving window analysis in circular buffers corresponding to different scales of landscape effects on the focal species (i.e. home range size and population levels). Predictors included in CIRC.MS model; (4) Friction-based and multi-scale predictors calculated in friction-based buffers or weighted distances to account for resource accessibility and landscape resistance to movement. These predictors were calculated only to assess suitable resource accessibility. Predictors included in FRIC model (see Figs. 1 and 2).
Breeding sites (wetlands) and woodlands are strongly selected by amphibian species during different stages of their life-cycle (Laan and Verboom 1990;Guerry and Hunter 2002;Ficetola and De Bernardi 2004;Cushman 2006;Zanini et al. 2008;Hartel et al. 2010;Boissinot et al. 2019) and the true availability of these resources for individuals is conditioned by their RedL Region: regional red list of Pays-de-la-Loire; conservation status: vulnerable (VU) and least concern (LC); Cons. Priority: regional priority level for species conservation from the regional red list defined by Marchadour et al. (2009). Nb pres. data: number of presence-only data used for SDM calibration. References associated with migration and dispersion range are shown Online Appendix 1 (mainly from Smith and Green 2005;Matos et al. 2019;Sinsch 2014) and ND refers to ''no data available''. Scales 1 and 2 refer to the maximal distances used for metric calculation based on ecological processes closely related to migration or dispersion respectively accessibility. In FRIC, we calculated predictor layers of suitable resources with a moving functional window analysis using a cost-distance method based on landscape permeability to obtain a functional grain  (Ray et al. 2002;Joly et al. 2003;Janin et al. 2009). A functional window is defined as a buffer around a pixel centroid (referred to as available area Fig. 1) whose shape is determined by landscape resistance (i.e. a friction-based reachable zone, see Ray et al. (2002)). Given that ponds influence ecological possesses at different levels (home range and meta-population levels), we used two pond variables: distance to the nearest pond, as the main breeding habitat requirement, and pond density within a potentially accessible area for dispersing individuals, as pond density is known to improve landscape connectivity (Cushman 2006;Ribeiro et al. 2011;Arntzen et al. 2017). In FRIC, these two variables were calculated using the friction cost methods described below. A resistance map (20 m resolution) was computed using friction cost values to account for resistance or barrier effects of different categories of roads, railways and urban areas (see Table 3). The cost values of impermeable barriers (i.e. high-density urban area, highspeed train line, highway and dual carriageways) were set as the maximum number of pixels that the species is able to cross, as required for analysis with Chloe 4.0 software (Boussard and Baudry 2017). Permeable features facilitating movement across barriers (e.g., viaduct, wildlife pass, and others, see Table 3) were digitized from aerial photographs. The allocation of resistance values is strongly debated, can and FRIC models. Example for metrics associated with deciduous and mixed woodlands. Here, the available area or distance calculation is shown for pixel A only greatly affect model quality and scientific evidence is limited. In addition, the resolution of our resistance map was too coarse to consider fine linear elements of the landscape used by species for movements (e.g. hedges, field margins, ditches, road verges, channelized agricultural headwater streams). As a precaution, we allocated lower cost values than usually reported to permeable roads, railways and low settlement density areas, and other landcover classes were considered as permeable (i.e. with cost value equal to 1) (see Table 3).
Lastly, for CIRC.MS and FRIC, maximum distances from the buffer centroid were based on species migration (i.e. seasonal movements in the annual home-range-scale 1) and dispersal (i.e. long-distance occasional movements-scale 2) (Fig. 2). These distances have been reviewed by Smith and Green (2005) and Boissinot (2009) and variable scale effects have Fig. 2 Multi-scale functional grain with or without landscape resistance. For this example, friction-based functional grain was calculated using barrier elements only. Rough estimates of maximum window sizes at each scale (S1 or S2) are classified according to species ecology, distance and scale of effect referred in previous studies (see Table 1 and Online Appendix 1) Other landcover types 1 Maximum value max was defined for barrier elements (impermeable) according to maximal number of pixels that the species is able to cross in the 20 m resolution resistance map been highlighted by Boissinot et al. (2019) in similar geographical and biological contexts. Pond density, road density, urban area and crops mainly influence population-level processes such as dispersal (Fahrig et al. 1995;Cushman 2006;Hartel et al. 2010;Ribeiro et al. 2011;Arntzen et al. 2017) and consequently related metrics were calculated at scale 2. B. spinosus and H. arborea are two species with higher mobility than other species so maximum distance considered for scale 1 (home range) was 600 m and for scale 2 (dispersion) was 2000 m (see Table 1). For all other species, we used a maximum distance of 300 m for scale 1 and 1000 m for scale 2 (Table 1 and Online Appendix 1). All variables were mapped using R environment v.3.5.3 (R Core Team 2019), QGIS 3.10 (Quantum GIS Development Team 2019) and Chloe 4.0 Software (Boussard and Baudry 2017). Variables used are described in Table 2. Final predictor resolution was 100 m to match presence data resolution and was obtained by resampling original raster layers (value of the pixel centroid).
The spatial correlation between environmental predictors was investigated using the variance inflation factor (VIF) as a measure of multicollinearity and Pearson correlation tests with VIF \ 6 and |r| \ 0.6, lower than advised by O'Brien 2007 (see Online Appendix 3).

Statistical models
Detailed modelling methods are described in Matutini et al. (2021). We used one regression-based approach (Generalized Additive Models, GAM) and one machine learning algorithm (Random Forest, RF) to predict and assess habitat suitability within the studied region. We used random pseudo-absence selection constrained to take sampling effort into account as follows. We counted the number of dates with at least one amphibian observation on a 500 m grid over the studied area. 500 m-cells were considered as sampling units, i.e. if one species was detected in the cell, we considered the 500 m-cell as visited. Then, we proportionally sampled pseudo-absences according to sampling effort, with a minimum distance of 500 m from a presence point of the target species. In addition, we randomly sampled 5% of pseudo-absences in unsuitable habitat (e.g. at least 75% of the cell area covered by urban and/or crops area without any tree or water in a cell). The number of pseudo-absences was fixed equal to the number of presence data (Barbet-Massin et al. 2012;Liu et al. 2019) and we ran 20 replicates of the pseudo-absence generation processes.
Finally, we conducted ensemble modelling using the median value of all generated individual maps (i.e. 20 maps/algorithm) and extracted a standard deviation map for each species.

Model validation
For independent model validation, we used two external independent presence-absence datasets from citizen sciences and field complementation (Matutini et al. 2021). The first dataset, called EVAL, was filtered to increase independence from the training-set and robustness of evaluation as described in Matutini et al. (2021). We firstly applied a 500 m minimum distance between pixels containing data to reduce spatial autocorrelation (repeated for each of the 500 iterations). Secondly, if a pixel had been monitored but the focal species were not detected, we considered the cell as ''absence data'' if: (1) a minimum sampling effort had been applied (according to focal species and observer expertise level); (2) at least 1 site without fish had been sampled in the cell (fish are main predators of amphibians). Datasets used for evaluation poorly overlapped calibration sets in space and time (See Matutini et al. 2021 for details). For the second evaluation set, called EVAL_STRAT, we used a stratified sampling of EVAL on final predictive maps for each species (see Matutini et al. 2021). This method produced better data distribution along the suitability gradients, as recommended by Newbold et al. (2010), Guisan et al. (2017).
We calculated different evaluation metrics to improve discrimination between models with 500 bootstraps each. We calculated a Symmetric Extremal Dependence Index (SEDI, Wunderlich et al. 2019) using the sediWeighted function in R from Wunderlich et al. (2019) and AUC values, Kappa, specificity (true negative rate) and sensitivity (true positive rate) of ensemble models with EVAL and EVAL_STRAT for each species (Allouche et al. 2006) using PresenceAbsence package (Freeman 2015). The threshold value for presence-absence discrimination were calculated for each permutation using the threshold function from dismo package (Hijmans et al. 2020) with sensitivity and specificity values optimisation strategy (i.e. True Skill Statistics, TSS).
In addition, we compared outputs from this modelling with final suitability maps obtained by Matutini et al. (2021) for the same species with coarse resolution (i.e. 500 m) and single-scale SDM.

Results
The model accuracy varied among species, with AUC values from 0.63 to 0.89, Kappa from 0.29 to 0.63 and SEDI from 0.57 to 0.89. The accuracies were higher for specialist species models compared to those of generalist species. The models from this study performed better than the models from Matutini et al. (2021) for all species except L. helveticus (Table 4). In addition, the selected models are more similar and stable between evaluation methods for forest specialist species and Urodeles except for T. cristatus.

Distance versus compositional metrics (DIST vs CIRC.SS): Hypothesis 1 (H1)
Using more complex landscape composition metrics generally improved model accuracy in comparison with simple distances to the nearest habitat. The predictive power of models using landscape composition predictors (CIRC SS or MS and FRIC) was better than models using distance-based predictors for all species except R. dalmatina and R. temporaria, with more ambiguous results depending on the evaluation method used (see Table 4 and Online Appendix 4).

Single-scale versus multi-scale models (CIRC.SS vs CIRC.MS): Hypothesis 2 (H2)
Comparison between single-scale (SS) and multi-scale (MS) models shows more ambiguous results. MS models performed better for four species, T. cristatus, B. spinosus, R. dalmatina and P. punctatus while SS models performed better for T. marmoratus, L. helveticus, H. arborea and R. temporaria (see Table 4). Discrimination between the two models is not possible for S. salamandra.

Functional grain including fine-scale movements or not (CIRC.MS vs FRIC): Hypothesis 3 (H3)
Using functional grain accounting for landscape resistance to fine-scale movements (FRIC) slightly improved model accuracy in comparison with CIRC.MS including predictors calculated used circular sliding windows. Indeed, differences between evaluation metrics varied from 0 to 0.02 for AUC, from 0 to 0.05 for SEDI and 0 to 0.07 for Kappa. FRIC models performed better for three species, L. helveticus, H. arborea and R. temporaria while CIRC.MS models performed better for T. cristatus only (see Table 4). Model discrimination is more difficult for other species but FRIC slightly improved predictions for P. punctatus and S. salamandra. In addition, FRIC model improved sensitivity (proportion of observed presences correctly predicted with TSS optimization) for most species (see Table 4).
Differences in habitat suitability predictions resulting from the four models (DIST, CIRC.SS, CIRC.MS and FRIC) are shown in Fig. 3 and Online Appendix 4. Absolute differences between suitability index values obtained through DIST and CIRC.SS modelling varied from 0 to 0.6, especially near cities. Absolute differences between CIRC.MS and FRIC were generally low, though in a few small areas scattered across the region the difference in suitability could be as high as 0.5. In addition, general patterns of suitable habitat differed between final maps of DIST and other models but were similar between CIRC.SS, CIRC.MS and FRIC (Fig. 3).

Discussion
For most species, model accuracy was satisfactory (AUC from 0.61 to 0.91) and higher and more stable for specialist compared to generalist species, as found in numerous empirical studies (Hernandez et al. 2006;McPherson and Jetz 2007) and simulations (Connor et al. 2018). Indeed, modeling the distribution of generalist species is particularly challenging (Brotons et al. 2004) because generalists are eclectic feeders, adaptable to various ecological contexts, use a large proportion of available habitats and are usually more mobile than specialist species. Hence, the statistical ''signal'' of the relationship between species' occurrence and environmental conditions Table 4 Performance of discrimination capacity and accuracy of ensemble models for each studied species Results obtained with stratified data for external evaluation (EVAL_STRAT) with 500 permutations. DIST: model with distance predictors only; CIRC.MS: model with compositional predictors calculated in different circular buffer sizes according to species migration or dispersal ability (multi-scale); CIRC.SS: same as CIRC.MS but using only one circular buffer size (single-scale); FRIC: same predictors as CIRC.MS but suitable resources are considered in a friction-based buffer according to landscape permeability. Best models: H1: distance vs composition-based models; H2: single-scale vs multiscale models; H3: including vs not including finer-scale movements using a frictionbased functional grain. Brackets indicate high uncertainty for model selection (low discrimination). Bold values in ''Best model'' shows model selected using the two evaluation sets. Results for EVAL dataset are show Online Appendix 4. For each species, the bottom line shows the evaluation metrics calculated on the final maps obtain by Matutini et al. (2021) Fig. 3 Differences between predictions of habitat suitability with DIST, CIRC (SS and MS) and FRIC models for two species of high regional conservation concern: Triturus marmoratus and Pelodytes punctatus. DIST model with distance predictors only; CIRC model with compositional predictors calculated in circular buffers (SS single scale; MS multi-scale); FRIC same predictors as CIRC.MS but suitable resources are considered in a friction-based buffer according to landscape permeability. |model 1 -model 2| with associated maps are absolute differences between the two final suitability maps of the two models may be too weak to capture (Brotons et al. 2004;Connor et al. 2018). On the other hand, model accuracy of forest specialist species may be overestimated because forests are relatively scarce habitats and patchily distributed in the region (see Brotons et al. 2004).

Weakness of distance predictors
Models using landscape composition predictors performed better than models with only distance-based predictors. These results are consistent with numerous studies from Landscape Ecology (Moilanen and Nieminen 2002;Martin and Fahrig 2012). Distancebased layers only represent proximity to elements without considering the amounts of habitat which are important for species persistence (Moilanen and Nieminen 2002;Martin and Fahrig 2012). However, distance variables might complement compositional landscape predictors for some species and/or landscape elements. In our study, we maintained distance to breading sites and to permanent rivers in CIRC and FRIC models because of the importance of breading site proximity and the landscape gradient induced by permanent rivers (e.g. sandy soil, disturbances, probability of fish colonization and flood events, etc.). For other species (e.g. for bats or sea birds), certain rare landscape elements associated with suitable breeding or wintering sites (e.g. a cave with a bat colony), which are central in their life-cycle, distance variables may be much more informative. Indeed, habitat selection of bats (e.g. Dawo et al. 2021)) and some seabirds (e.g. Bolton et al. 2019), as well as others social species (e.g. social bees), are greatly influenced by distance to colony.
Combining multi-scale predictors for a functional approach Differences in accuracy between single-scale and multi-scale models were more ambiguous and varied among species, especially for generalist species. Hence, for B. spinosus and T. cristatus the multi-scale model performed best, while for H. arborea and L. helveticus the single-scale model was more accurate. For predictive approaches with SDM, the literature suggests that multi-scale or multi-level models are more accurate than single-scale models (Graf et al. 2005;Boscolo and Metzger 2009;Fournier et al. 2017;Bellamy et al. 2020;Hallman and Robinson 2020) but some authors have also found ambiguous results (Martin and Fahrig 2012). In our study, we used an a priori selection of appropriate scale for each predictor according to a species-specific known scale effects in similar geographical contexts but this might be misleading (Hallman and Robinson 2020). Although species mobility and some associated species traits as body size can be good indicators of potential variable scale of effect (Jackson and Fahrig 2012;Thornton and Fletcher 2013), scale-effects are much more complex and vary across large geographical extents. From simulation, Connor et al. (2018) demonstrated that scale-of-effect can vary between different landscapes (e.g. homogeneous vs heterogeneous), models and species and that statistical signal is difficult to capture for generalist species in homogeneous landscapes and with increasing grain size.
In addition, multi-scale amphibian movements are still poorly understood (Pittman et al. 2014) and dispersal distances may be under-estimated (Sinsch 2014). According to Miguet et al. (2016), the scale of effect may be influenced by five groups of factors: species traits, landscape variables, biological response, indirect influence (e.g. water quality or predator density) and regional context of the study. However, the determinants of the scale of effect are still only partly understood and empirical results tend to contradict theorical predictions making a priori scale selection particularly complex (Jackson and Fahrig 2015;Moll et al. 2020). Therefore, even with good prior knowledge on species ecology in neighboring and/or similar regions, transferring results of previous studies may be erroneous and defining scaleof-effect using species mobility capacity might be not sufficient (Mayor et al. 2009;Zanini et al. 2009). Furthermore, the model itself and the algorithm used might be sensitive to different scale-of-effect for the same variables (Connor et al. 2018;Rose et al. 2020;Hallman and Robinson 2020). Indeed, Rose et al. (2020) showed model-specific sensitivity to scale when predicting the distribution of a toad species in a large area in California using Random Forest, Point Process Models and Maxent. Some other algorithms such as Boosted regression trees (BRTs) are more robust to overfitting and collinearity then allow us to use several predictors at different scales (Hallman and Robinson 2020).

Functional grain and finer-scale movements in SDM
Including landscape resistance to movement produced only slight to moderate improvements in model accuracy for six species. This may have been a consequence of the evaluation-set. Indeed, the majority of sites used to evaluate the models were far from potential barriers to movements (cities and roads) which are the main elements mapped in the friction map. As a result, the evaluation sites are mainly located in areas with similar environmental characteristics and low variability between predicted suitability values from models including or not landscape resistance (i.e. CIRC.MS and FRIC).
Secondly, we mainly considered barrier elements in the friction map, which represent a small proportion of the total study area, also leading to low variation between predictors of CIRC.MS and FRIC. Many studies at local scale use several land-cover classes to generate resistance maps, in particular crops, forests, wetlands and edges (Ray et al. 2002;Joly et al. 2003;Compton et al. 2007;Decout et al. 2012;Janin et al. 2012;Youngquist and Boone 2014;Jeliazkov et al. 2019). However, the integration of other land-cover classes to develop a more precise permeability map must be carried out at higher resolutions (e.g. 5 m). This relies on sufficient computational power for metric calculation at regional extent and on knowledge for resistance value estimations (see Zeller et al. 2012;Keeley et al. 2017). Indeed, amphibians can use small or narrow elements of the landscape such as hedgerows, ditches, field margins, road verges and channeled agricultural headwater streams for movements (Pope et al. 2000;Mazerolle 2005).
Considering these elements, the slight improvements observed indicate that using functional grain with a friction-based approach might be relevant for predictions at large scales and merit further investigation, especially as crucial movement processes operate at fine-scale. Ray et al. (2002); Zanini et al. (2008) have underlined the interest of friction-based modeling approaches to improve the prediction of B. bufo presence or B. bufo and R. temporaria respectively compared to circular buffer zones around ponds. Other studies underline the importance of considering true accessibility of habitat to explain species distributions or population dynamics and avoid incorrect conclusions (Eigenbrod et al. 2008;Hamer 2018). Our approach, for a larger number of species, highlights the potential of this method at any point of the landscape (i.e. not only around ponds).

General limits of our study
Some ecological and methodological limits need to be mentioned. Firstly, the data we used are unstandardized presence-only data, from different sources and mainly collected by citizens. Data quality therefore limits the study of fine-scale processes and results need to be interpreted with caution (Guillera-Arroita et al. 2015); conversely, citizen science is a requisite to have data over large areas. Secondly, our large-scale study involved a diversity of landscapes that may present a high local level of variable collinearity undetectable over a large extent. These different landscapes are the expression of different agricultural practices (e.g. cattle farming vs cereal crops; extensive practices vs intensive agriculture), which significantly affect amphibian distributions at different scales (e.g. water pollution, decrease and contamination of food resources, hedge management more or less suitable, spatial heterogeneity). The associated environmental data were not available at regional scale. In addition, the wide extent of our study implies the presence of sub-specific variations in habitat selection and movements (i.e. local adaptations) reducing model accuracy in species-level analyses (McPherson and Jetz 2007).

Conclusion
We highlight the potential of using a friction-based functional grain to improve SDM prediction at the species level but further investigation is needed. Our results suggest that a multi-scale approach using a priori selection of scales, based on empirical data on species mobility and scale-of-effect, leads to ambiguous results depending on the species studied, in particular for generalists at species level. New methods are currently under investigation to propose a framework for multi-scale and multi-level SDM using specific scale-optimization or pseudo-optimization methods (Stevens and Conway 2019; Hallman and Robinson 2020) or post-modelling analysis (Miguet et al. 2016;Moll et al. 2020). Furthermore, spatial variability of scale-of-effect and species mobility through landscape matrices is a major limit for species-level analysis. Spatial structure of data according to landscape or sub-species variations induced by local adaptation is currently poorly included in SDM studies but should be considered explicitly, especially for high resolution studies over large extents (see Swanson et al. 2013;Crase et al. 2014). However, ecological relevance and final objectives of the study should remain the priority for model conceptualization and predictor selection. More complex models may be unreasonably computationally intensive for specieslevel analysis over large extents and could lead to difficulties in model interpretation, discrimination and evaluation (e.g. Bell and Schlaepfer 2016).
Decisions for biodiversity conservation often call for study of habitat suitability at large (regional) or global (national or European) scales. The ambiguity of certain results presented in this paper show that it can be difficult to include fundamental ecological processes linked to scales-of-effect and movements over large areas (i.e. with sub-specific variability, diversity of landscapes and agricultural practices, mainly heterogeneous data from citizen science). As a result, a balance must be found between research and application objectives, availability and quality of data (environmental and biological) and model complexity with a strong consideration of ecological justifications. Data availability Data sample and access procedure are available online: https://doi.org/10.5281/zenodo.4358147.

Declarations
Conflict of interest The authors of this preprint declare that they have no financial conflict of interest with the content of this article.
Ethical approval Access to citizen data is subject to a user agreement between our lab and the associations involved. We have done our best to involve volunteers and associations in this work with regular exchanges (attendance at meetings, reporting of results in different forms, organization of steering committees). In addition, all people involved in handling amphibians had specific ministerial authorization.
Consent to participate All participants in this study consented.
Consent for publication All authors contributed critically to the drafts and gave final approval for publication.