The inﬂuence of road networks on brown bear spatial distribution and habitat suitability in a human-modiﬁed landscape

Roads are human infrastructure that heavily affect wildlife, often with marked impacts on carnivores, including brown bears Ursus arctos . Here, we assessed the potential impact of road networks on the distribution of brown bears in the small, isolated and endangered Cantabrian population of north-western Spain. To ascertain whether local road networks affect brown bear spatial distribution, we ﬁ rst assessed potential in ﬂ uences on the distance of bear locations to roads using candidate models which included topographic variables, landcover types, bear age and reproductive status, traf ﬁ c volume and road visibility. Then, we built two sets of habitat suitability models, both with and without roads, to discern the possible loss of habitat suitability caused by roads. The mean distance of bear locations to the nearest road was 968 (cid:1) 804 m and the closest road was a low traf ﬁ c road in 72.5% of cases. Candidate models showed little in ﬂ uence of our variables on bear distance to the nearest road, with the exception of elevation. Habitat suitability models revealed that road networks in our study area seem to have almost no effect on brown bear habitat suitability, except for females with yearlings during the denning season. However, this result may also be a consequence of the fact that only a small proportion (16.5%) of the cells classi ﬁ ed as suitable bear habitats were crossed by roads, that is, most of the roads are primarily located in unsuitable bear habitats in the Cantabrian Mountains. Compared to previous studies conducted in other populations, mainly North American ones, our ﬁ ndings might suggest a different response of Eurasian brown bears to roads due to a longer bear-human coexistence in Europe versus North America. However, the indirect approach used in our study does not exclude other detrimental effects, for example, road mortality, increased stress and movement pattern disruption, only detectable by more direct approaches such

Wide-ranging mammals with low reproductive rates and low densities (e.g. large carnivores) are particularly vulnerable to the various effects of roads and vehicle traffic (Alexander et al., 2005;Kautz et al., 2021;Rytwinski & Fahrig, 2015). For example, roads can negatively impact tiger Panthera tigris survival and reproductive rates (Kerley et al., 2002). They can also influence, through road avoidance, the movement behaviour of jaguars Panthera onca (Colchero et al., 2011) and wolves Canis lupus (Ciucci et al., 2018;Dennehy et al., 2021). Indeed, the fragmentation of carnivore populations caused by roads (Fahrig & Rytwinski, 2009) can affect the functioning of entire ecosystems due to the ecological role large carnivores play as apex predators (Ordiz et al., 2013(Ordiz et al., , 2014. Among large carnivores, the multiple, non-exclusive effects of roads on brown bears Ursus arctos have been extensively studied, mainly in North America. For example, avoidance patterns (Proctor et al., 2020;Støen et al., 2020) and fast displacement rates (Kite et al., 2016;Roever et al., 2010) have been described in areas surrounding roads, with bears selecting higher elevations and steeper slopes because they are further away from roads and, consequently, less accessible to humans (Goldstein et al., 2010;Nams et al., 2006). The impact of road networks also seems to depend on landcover types surrounding roads and periods of the bear life cycle, with higher probabilities of road crossings where roads intersect areas with dense vegetation cover offering shelter (Find'o et al., 2019;Lyons et al., 2018;Roever et al., 2010) or during hyperphagia (Fra z ckowiak et al., 2014; Stewart et al., 2013), due to the need to consume a large amount of food before hibernation. Similarly, roads have also been described as occasional attractants to brown bears as they facilitate communication (Gonz alez-Bernardo et al., 2021;Sato et al., 2014) or movement (Hill et al., 2021;Roever et al., 2010). Along the same line, some age classes such as females with cubs and subadults may use roads to avoid risky encounters with adult males (Graham et al., 2010;Penteriani et al., 2018), as also occurs in other ursids (Gantchoff et al., 2019). Moreover, traffic volume associated with the type of road has also been shown to determine the severity of road impacts on bears (Elfstr€ om et al., 2008;Northrup et al., 2012). In fact, increasing traffic intensities are often associated with stronger avoidance (Jacobson et al., 2016;Mace et al., 1999) and decreasing permeability (Find'o et al., 2019;Skuban et al., 2017), affecting movement rates (Proctor et al., 2012;Roever et al., 2010) and altering rhythms of activity (Støen et al., 2020;Waller & Servheen, 2005). Finally, one aspect of roads that might affect bears and has never been explored as a factor of disturbance (but see the analyses of traffic kills sites in Huber et al., 1998) is road visibility (i.e. a viewshed, the part of the environment assumed to be visible to an animal in a given location, Tandy, 1967). In recent years, it has been advocated that the quantification of an animal's potential visual space may improve our understanding of animal use of an environment (Aben et al., 2018), but still few examples exist on the link between visibility and large carnivores (Davies et al., 2016;Grant et al., 2005).
Here, we assessed the impact of the road network on the spatial distribution of brown bears in the small, isolated and endangered population inhabiting the human-modified landscape of the Cantabrian Mountains in north-western Spain. We predicted that bear locations would be closer to roads when: (1) individuals were in areas with dense vegetal cover compared to open habitats; (2) roads were not visible; (3) traffic level was low; (4) females had cubs during the mating season to decrease the risk of infanticide from encounters with adult males; and (5) during hyperphagia, as bears might use food resources associated with road edges. Predictions (1), (2) and (3) are related to a lower expected level of disturbance in these cases and therefore to less avoidance of roads by bears, while predictions (4) and (5) are related to specific seasonal circumstances. Furthermore, we tested the possibility that road networks may reduce habitat suitability for bears, an effect that we expect might be different depending on bear classes and seasonal cycles Penteriani et al., 2018).

Study area and bear population
The study area encompasses the distribution range of the Cantabrian brown bear population in the provinces of Le on, Palencia and Asturias, north-western Spain (Fig. 1), which occupies an area of 4474 km 2 along the Cantabrian Mountains at an average elevation of 1100 m a.s.l. (maximum 2650 m; Penteriani et al., 2019). The climate of the region is oceanic, more continental and drier along southern slopes and temperate and more humid on northern slopes (Ortega & Morales, 2015). The region is mainly covered by forests, shrublands and farmland. The forests of the southern slopes are mainly composed of semi-deciduous and evergreen oaks (Quercus sp.), whereas the northern slopes host mostly deciduous forests (Fagus sylvatica; Q. robur, Q. petraea; Betula sp.; Garc ıa de Celis et al., 2004;Mateo-S anchez et al., 2016). Non-forested areas are covered with shrubs such as broom (Cytisus sp.) and heather (Erica sp., Calluna sp.), while above the tree line, berry shrubs such as bilberries appear (Vaccinium myrtillus; Mateo-S anchez   Pato & Obeso, 2012). In general, most suitable bear habitat patches are embedded in a matrix of urbanized and cultivated areas with a high density of transport routes and human settlements, where the main economic activities are livestock breeding, recreational activities, mining and timber harvesting . Although the Cantabrian brown bear population trend in recent years has been positive (Gonzalez et al., 2016), its size is still small (c. 250 individuals 10 years ago [P erez et al., 2014], and probably larger now due to this trend), and it continues to be considered as threatened (Boitani et al., 2010;MITECO, 2022). The road network in the study area consists of a few high-and medium-capacity roads that connect the main urban areas and a majority of minor local roads that connect small mountain towns with each other or with the main roads. This road network has hardly changed in the last decades in the two regions studied, Asturias and Castilla y Le on (increases of 0.56% and 1.89%, respectively, between 2005 and 2019; MITMA, https://www.mitma. gob.es/recursos_mfom). Viewing points used by rangers and ourselves, as well as camera trap locations, were evenly distributed over: (a) the entire bear range in the study area Zarzo-Arias et al., 2019); and (b) the whole year, winter included, as bear hibernation is facultative in the Cantabrian Mountains (Gonz alez-Bernardo, Bombieri, et al., 2020;Uzal et al., 2022). Each observation was accompanied by information on the date and the age class of the observed bear, that is, adult or subadult solitary bear, female with cubs, female with yearlings. We classified the bear locations into three seasons, according to the annual cycle of this population: denning (January 1 to April 15), mating (April 16 to June 30) and hyperphagia (July 1 to December 31; Zarzo-Arias et al., 2021). Due to the mild winters of the area, food is also available during the winter months and not all bears in this population hibernate (Gonz alez-Bernardo, Bombieri, et al., 2020;Gonz alez-Bernardo, Russo, et al., 2020).

Model covariates
On the basis of the information obtained from previous studies on this bear population (Lamamy et al., 2019; Zarzo-Arias et al., 2019), we selected variables that may influence the distance of bear locations to the nearest road, including: topographic variables (elevation, slope, terrain ruggedness), landcover type, season, bear class, traffic volume and road visibility (Tables 1 and 2). Landcover information was extracted from the CNIG (Environmental Thematic Cartography of the Principality of Asturias; Sheets of the Map of Vegetation, Lithology, Rocks and Habitat of the Bear, 2011. Scale 1:25 000. © Principado de Asturias, Spain) and reclassified into seven landcover types: (1) forests, (2) shrubs, (3) pastures, (4) farmland, (5) rocky areas, (6) bare ground/water and (7) urban areas. The first six landcover types were used in the construction of the distance to road models (see 2.4.1 Road network effects on brown bear spatial distribution), while for the habitat suitability analyses we also included urban areas (see 2.4.2 Potential impact of road networks on brown bear habitat suitability).
The road layer of the study area was obtained from the CNIG, and the closest road and the closest visible road were selected for each bear location. These roads were classified into three groups based on the average daily volume of traffic, with data obtained from the public repositories of the Spanish Ministry of Public Works and Transport (MITMA, https:// www.mitma.es/carreteras), the Castilla y Le on government (Junta de Castilla and Le on, https://carreterasytransportes.jcyl. es/web/es/carreteras-transportes.html) and the Principado de Asturias government (SADEI, https://sadei.es/inicio). The unit of traffic volume was daily traffic intensity (IMD, i.e. vehicles/day, hereafter, v/d), taking as the value for each road the mean value of the average monthly IMD of the years 2012 to 2018, as these were the years with the most complete and modern data. The types of roads based on the IMD were established according to the 'natural breaks' method, which is based on the nature of the data and identifies important jumps in the sequence of values, optimizing the grouping of similar values and maximizing the differences between classes (Jenks, 1967). The road types for the roads considered according to the ranges of IMD values established by the method are the following: Type 1 = 0-337 v/d; Type 2 = 338-1411 v/d and Type 3 = 1412-5129 v/d.
The topographic variables and viewsheds (to assess visibility) were calculated from a Digital Elevation Model and a Digital Surface Model respectively. Both models had a 30 m resolution and were based on a LiDAR (Light Detection and Ranging) data layer obtained from the CNIG (https:// centrodedescargas.cnig.es/CentroDescargas/index.jsp), which has a minimum density of one point every 2 m and an altimetry precision of 15-20 cm. These data allow for the derivation of more precise viewsheds (Aben et al., 2018;Lagner et al., 2018), as they incorporate elements of the physical environment that may interfere with the field of vision of the observers (e.g. vegetation, topographic features or infrastructures).
To obtain visibility of the nearest road from each bear location, a raster of visible cells of the surface from one or more bear locations was calculated (i.e. the viewshed), with a range of 1000 m (visible cells = 1, non-visible cells = 0). The surface model takes into account the real height of vegetation and topographic elements through point clouds with threedimensional coordinates obtained by airborne LiDAR sensors (Lorite et al., 2015). As far as we know, there are no published data on the visual range or acuity of brown bears. However, different thresholds of proximity to roads have been described in which road avoidance can be detected, most of them being at least 500 m (Mace et al., 1996;McLellan, 1989;Torres et al., 2016;Waller & Servheen, 2005). As a result, we decided to conservatively set our road visibility threshold to at least twice the distance proposed for grizzly bears. Next, we merged viewsheds and the road layer (along with road type information), whereby road sections were assigned values ranging from 1 to 3 depending on the traffic load, with a value of 0 for non-visible ones. Finally, we calculated the distance from each bear location to the nearest road, which represented our response variable. All spatial variables (topographic, landcover, visibility and distance from bear location to the nearest road) were calculated using ArcGIS 10.5 software (ESRI, Redlands, USA).

Road network effects on brown bear spatial distribution
We built a candidate set of linear mixed models with the distance of each bear observation to the nearest road as a response variable and the topographic, landcover, traffic level and visibility variables as additive predictors. We logtransformed the response variable to approximate normality in the residuals and re-scaled the numeric predictor variables to mean = 0 and standard deviation = 1, in order to standardize their values. To account for spatial autocorrelation of bear presence data, we included in our model the autocovariate (ac) using the spdep R package (Bardos et al., 2015;Bivand et al., 2019). To exclude collinearity among the predictors, we calculated the variance inflation factor (VIF; Belsley et al., 2005), but no variable exceeded a threshold of 2 (maximum VIF = 1.57). Models were built using the lme4 R package (Bates et al., 2015) and once we generated the sets of competing models, we employed the Akaike information criterion (AIC), using the value of DAIC < 2 as the criterion for selecting the most parsimonious model (Anderson & Burnham, 2004;Richards et al., 2011). All analyses were performed using R v. 3.5.1 statistical software (R Core Team, 2018). Concerning our model, we assessed whether our visual observations, which comprise most of the location data, were biased towards open areas by comparing the same model built with and without visual observations. As both models produced comparable results, we believe observational bias, if any, has a negligible effect on model performance (see Fig. S1). To detect such a bias deriving from our data (consisting mainly of visual observations) we built a model with visual observations and another with only locations based on footprints, which are not biased by detectability since they are not more easily detectable in open terrain. Although the second model had few observations (n = 189), especially in some categorical variables due to the small number of cases, we observed similar trends and responses of the variables, thus discarding biases, and therefore we built the distance to road models with all the observations (Fig. S1). According to the AIC-driven model selection, a random year in the model was not deemed necessary.

Potential impact of road networks on brown bear habitat suitability
To assess habitat suitability for brown bears in relation to the study area road network, we used the software MaxEnt version 3.3.3 k called from the R environment with the packages dismo (Hijmans et al., 2017) and ENMeval (Muscarella et al., 2014) following the methodology applied in Zarzo-Arias et al. (2019). To identify an optimal model structure for each input bear dataset, we evaluated candidate models with all types of feature combinations, each run over a set of regularization multipliers ranging from 0 to 19 (for more detailed information see Zarzo-Arias et al., 2019). We included the same variables as in the previous model (landcover, elevation and slope) and distance to the nearest road at a 500 m scale, but removed ruggedness due to its high correlation with slope to avoid any potential multicollinearity (VIF < 5 for the remaining variables; Zuur et al., 2010). Following the aforementioned methodology, we also included urban area as an additional landcover, not represented in the distance to road models (see Table S1). We made a model for each bear class in each season considering as a presence all 500 m grids with at least one occurrence, except for females with yearlings in hyperphagia as we only had four locations. In these models, we used 5000 iterations, a convergence threshold of 10 À5 , the checkerboard1 method for data partitioning, 50 replicates and centre coordinates from all cells in the study area as background to build the models.
We identified the best combination of feature types and regularization multiplier using AICc. We considered models within 2 AICc units of each other to have equivalent empirical  Table S2). We obtained the output maps of the best model for each group with the complementary log-log (cloglog) format, which allows interpreting values given to each cell from 0 to 1 as a likelihood of bear occurrence. We extracted all cells crossed by a road for each bear type in each season and analysed their occurrence likelihood (Max-Ent cloglog values). For selecting a threshold to identify the most suitable bear range, we extracted the values predicted for occurrence locations of each bear class used to build the models from their corresponding model and calculated the mean of the values predicted by the model for all bear occurrence points, rounded up to a 70% likelihood of bear occurrence (cloglog value = 0.7, Table S3). We selected this value as the threshold and identified those cells above this threshold traversed by the different types of nearest road (see Model covariates). Finally, we performed a niche comparison analysis between models including and excluding the variable distance to roads based on Warren's-I statistic (Warren et al., 2008), which ranges from 0 (no overlap) to 1 (identical).

Results
More than half (52.9%) of bear locations were of solitary bears, 39.7% to females with yearlings and 7.4% to females with cubs. Observations were more frequent during the mating and hyperphagia seasons (41.2% and 38.4% respectively), followed by locations recorded during the denning season (20.4%). Most of the bear locations occurred within forests (52.2%) and shrublands (28.9%) and, to a lesser extent, rocky areas (6.6%), farmland (5.7%) and pastures (5.5%). In relation to bear distance to the nearest road, most locations (72.5%) occurred next to low-traffic roads, followed by moderate-traffic (15.7%) and high-traffic (11.8%) roads. The mean distance (AESD) of bear locations to the nearest road was 968 m AE 804 m (min. = 0 m, max. = 4595 m) and the median (IQR) was 773 m (IQR = 1051; see Tables 1 and 2).

Road network effects on brown bear spatial distribution
The most parsimonious model showed that distance to the nearest road depended on bear class, season, landcover type, traffic level and visibility, as shown by the exclusion of zero in the 95% CI of the estimated parameter for these predictor terms (Table 3). Adjusted R-squared of the most parsimonious model with DAICc < 2 was 0.445. All these variables were present in the five competing models except for bear class, which was absent in one of them (Table 3). Specifically, females with yearlings and females with cubs were found closer to roads than solitary bears, with females accompanied by yearlings the closest to roads (Fig. 2a). Bears also were closer to roads during the mating season and they were found further away during hyperphagia and winter (Fig. 2b). Bears were closer to roads of intermediate traffic intensity (Fig. 2c) and closer to visible roads compared to non-visible roads (Fig. 2d). Also, distance to the nearest road was lower in farmland, pastures and bare ground than in forests (Fig. 2e). These effects were found while controlling for confounding variables known to affect distance to roads, such as elevation (variable with more weight in the models), slope or ruggedness (Fig. 2f-h).

Potential impact of road networks on brown bear habitat suitability
Most models including and excluding the variable distance to the nearest road presented very similar values in terms of variable contribution (Table S4) and habitat suitability, as almost all bear classes had a weak positive relationship with this variable (Fig. S2) and its contribution to the models was generally low (<1.5%), denoting a small impact of road proximity on bear habitat suitability (Fig. S3). In turn, the most important variables were slope, which showed the highest proportional contribution to all models, followed by elevation and forests. In contrast, the model for females with yearlings during the denning season showed a higher contribution of the variable distance to the nearest road (13.4%; Table S4).
The 500 m cells crossed by a road were not classified by the models as highly suitable habitat for bears (less than 16.5% of those cells had more than 70% likelihood of bear occurrence), with little difference between seasons (Table S5, Fig. 3). Despite this, models including and excluding the distance to the nearest road variable presented similar habitat suitability. Actually, the mean niche overlap (based on Warren's-I statistic) of cells crossed by a road between models was >0.99, showing a very high similarity.

Discussion
Here, we used two approaches to explore the potential impacts of road networks on the spatial distribution of brown bears belonging to a population inhabiting a human-modified landscape. Both approaches suggest that roads seem to have a rather small effect on bear spatial distribution and/or habitat suitability.
Several studies have previously described roads as a factor affecting brown bear behaviour (e g. Kite et al., 2016;Morales-Gonz alez et al., 2020;Proctor et al., 2020;Skuban et al., 2017), spatial distribution and habitat use, either by analysing bear distances to roads or by quantifying the density of road networks. However, in our study, we found that a large proportion of the variation in distance to the nearest road was explained by elevation. Taking into account the characteristics of the study area, where roads are generally located along valley bottoms as in other similar areas (Carpathians, Papp et al., 2022), this relationship seemed to be circumstantial since brown bears tend to be more frequently located in higher areas so as to avoid human settlements and infrastructure along valley floors (Fra z ckowiak et al., 2014;Goldstein et al., 2010). This possible explanation is also supported by the fact that elevation and distance to the nearest road presented a substantial Pearson's correlation (0.58) in our habitat suitability models.
Contrary to our expectations, habitat type, visibility of the nearest road and season did not seem to meaningfully affect brown bear distance to the nearest road. Moreover, our habitat suitability models were not affected by the inclusion of roads, except for the model that considered females with yearlings in the denning season. Despite the small sample size of this model (n = 93), we suggest that the impact of roads on the habitat suitability of females with yearlings during this season may be due to their higher activity level compared to other bear cohorts that are usually hibernating. In the southernmost populations and with mild winters, many bears do not hibernate, with the exception of pregnant females, and females with yearlings have been reported to hibernate less than all other bear classes (Gonz alez-Bernardo, Russo, et al., 2020). Additionally, as the habitat suitability models showed, road surroundings are frequently characterized by low suitability, probably because less favourable habitats (i.e. farmland) and human settlements and activities surround roads. Finally, visibility of the nearest road had a very small effect on the models. To our knowledge, this parameter has not been taken into account in other studies to determine avoidance from a bear's point of view, although it has been considered from a human perspective in the analysis of traffic mortality sites (Huber et al., 1998;Psaralexi et al., 2022). Also, visibility of the nearest road is related to the type of habitat in which the road is located. It has been reported that bears near roads preferentially use mature and dense forests during the daytime, which limits their visibility (Roever et al., 2010). Gibeau et al. (2002) and Roever et al. (2008a) observed that Canadian brown bears were closer to roads, or were more likely to cross them, if there were high-quality habitats around them. It has been suggested that closeness to roads may be related to the presence of attractive food (Huber et al., 1998;Morales-Gonz alez et al., 2020;Stewart et al., 2013), especially in spring and early summer (Roever et al., 2008b), as well as during hyperphagia (Graham et al., 2010;Zarzo-Arias et al., 2018). Yet, as roads increase mortality rates (Bourbonnais et al., 2014;McLellan, 2015), vehicle networks represent Table 3 (a) LMs built to investigate variables influencing distance to the nearest road in the Cantabrian brown bear population (n = 2722 locations). The table shows the 5 models with the lowest AICc. Variable codes: elevation = altitude above the sea level, slope = slope (%), class = class of bear, season = period of the bear lifecycle, land = landcover, type_near = type of the nearest road (by traffic intensity, IMD), vis_near = visibility of the nearest road. (b) Effects of explanatory variables on the most parsimonious model. For each explanatory variable, we report the estimate (b), standard error (SE), significance (P), confidence intervals (CI) and explained variance (R 2 ) of the most parsimonious model sink-like areas (Braid & Nielsen, 2015;Falcucci et al., 2009) or ecological/evolutionary traps (Ciarniello et al., 2007;Penteriani et al., 2018). This is particularly relevant since collision mortality is greater in near-road habitats with low visibility, such as those with high vegetation cover, or in areas with greater animal mobility and poor light (Psaralexi et al., 2022). The different bear classes express little dependence in relation to distance from the nearest road, contrary to our expectations. Brown bear sex and age classes have previously been described as showing different patterns of road avoidance and crossing, but they appear to be dependent on local features as described patterns are not consistent between studies (e.g. Chruszcz et al., 2003;Gibeau et al., 2002;Graham et al., 2010;Proctor et al., 2018;Roever et al., 2008b;Steyaert et al., 2016;Waller & Servheen, 2005). Actually, the lack of uniformity in sex and age responses to road networks may be due to local features of the habitats surrounding roads, which would offer different trophic resources and/or shelter throughout the year (Morales-Gonz alez et al., 2020). For example, females with cubs can select areas close to roads as a mechanism of adult male avoidance to prevent infanticide, which means that the relationship between certain bear classes and roads may be dependent on the local context of individual interactions (Graham et al., 2010;Penteriani et al., 2018, but see Waller & Servheen, 2005. However, although we expected this pattern to be repeated in the Cantabrian population, it was not the case. We suggest that females in the Cantabrian mountain range use steep and inaccessible areas as refuge from dominant males, especially in the first few weeks after exiting the winter den with cubs . In these areas, it is to be expected that there are no or very few roads, so their effect would not be very noticeable. In this study, in the solitary bear class, no distinction was made between females and males, so that some sex-related features might have been overlooked. Considering previous studies on the effects of road networks in other brown bear populations in Europe and North America, the observed different response of the species between these two continents stands out. Roads showed a more marked negative effect in North American populations but see G€ uthlin et al. (2011) andFra z ckowiak et al. (2014)). The historical coexistence of humans and brown bears in Europe, contrasts with the shorter and more intense interaction experienced by North American populations, where widespread human presence and activities are much more recent Støen et al., 2020). In Europe, human density and encroachment are also higher than in North America Swenson et al., 2000). Thus, different levels of exposure to human activity and persecution may have driven the observed different behavioural responses Ordiz et al., 2014). Actually, in a previous study on this Cantabrian population, it was recorded that bears do not seem to modify their surveillance behaviour with respect to the distance to roads . In other European populations: (a) bear occurrence was positively related to distance to roads in the Italian Central Apennines (Maiorano et al., 2019); (b) lack of habitat avoidance/attraction was recorded near roads in Slovenia (Kaczensky, 2000;Kaczensky  Brown bear road ecology E. Gonz alez-Bernardo et al. et al., 2003); (c) habitat selection in the reintroduced population of brown bears in the Pyrenees was not affected by road density (Martin et al., 2012) and (d) the proximity of roads in Romania appeared not to affect bear activity (Roellig et al., 2014). Thus, our approach seems to suggest that traffic intensity does not explain distance to the nearest road in the Cantabrian Mountains, despite the expected avoidance based on traffic level that has been recorded elsewhere in Europe (Find'o et al., 2019;Huber et al., 1998;Kaczensky et al., 2003;Skuban et al., 2017). We suggest that in our study area, which is primarily mountainous and lacks large urban areas, levels of vehicle traffic might not reach the avoidance threshold that makes bears sensitive to road presence.
In contrast, in North America (mainly British Columbia and Alberta in Canada and Alaska and Montana in the United States), greater avoidance of the area adjacent to roads has been reported (Ciarniello et al., 2007;McLellan, 1989;Waller & Servheen, 2005). This negative selection of habitats near roads has been described as occurring irrespective of traffic (Jacobson et al., 2016;McLellan & Shackleton, 1988) or as a traffic-dependent phenomenon (Chruszcz et al., 2003;Gibeau et al., 2002;Mace et al., 1996Mace et al., , 1999Northrup et al., 2012;Proctor et al., 2012), with low frequencies of road crossing when traffic level increases (Chruszcz et al., 2003;Waller & Servheen, 2005). Faster movements and lower activity rates have also been described (Proctor et al., 2012;Roever et al., 2010), as well as more nocturnal behaviour (Northrup et al., 2012;Waller & Servheen, 2005). These apparent differences between Europe and North America, therefore, suggest that projecting management or mitigation measures at the same spatial scales for Europe and North America might not be appropriate due to both the different sensitivity that bears on the two continents seem to exhibit to roads and the different degree of human landscape modification and encroachment.
Because discerning the role of road networks in large carnivore ecology is crucial, and particularly relevant for isolated and small populations like the Cantabrian one, we consider it important to stress here a potential limitation of our approach. This work is based on opportunistic observations and indirect evidence, which do not provide as accurate information about movement patterns and rhythms of activity around roads or potential road avoidance strategies as does telemetry. However, (a) the initial filtering of the data for reliability minimizes possible biases due to the opportunistic or non-professional collection of the data (Gantchoff et al., 2022), and (b) the yearround monitoring of bears in most of their distribution range within the Cantabrian Mountains should prevent spatial and temporal biases. As an end consequence, the information presented here has the potential to (a) have crucial implications in the management of this isolated and endangered bear population, for example, mitigation measures of human infrastructure; and (b) represent a first step towards detecting potential habitat loss in bear distributions and locations where road networks may have a greater impact at a spatial level. In animal populations inhabiting territories that have been highly humanized for centuries, such as the one studied here, the scale and entanglement level of human activities in bear habitat should always be taken into account. females with yearlings during the denning season, (e) females with yearlings during the mating season, (f) solitary adults/subadults during the denning season, (g) solitary adults/subadults during the mating season and (h) solitary adults/subadults during the hyperphagia season. The curves show that the predicted likelihood of presence changes as each environmental variable is varied, keeping all other environmental variables at their average sample value. Figure S3. Jackknife tests of variable importance for regularized training gain of the models for: (a) females with cubs during the denning season, (b) females with cubs during the mating season, (c) females with cubs during the hyperphagia season, (d) females with yearlings during the denning season, (e) females with yearlings during the mating season, (f) solitary adults/subadults during the denning season, (g) solitary adults/subadults during the mating season and (h) solitary adults/subadults during the hyperphagia season. Light blue bars represent the regularized training gain for each model without the variable, while dark blue bars represent the regularized training gain for each model with only this variable. The lower red bars represent the regularized training with all the variables. The environmental variable with the highest gain (biggest dark blue bar) when used in isolation contains the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted (smallest light blue bar) possesses the most information that is not present in the other variables. Table S1. Mean values of each environmental variable included in the MaxEnt models in the cells crossed by the different types of roads. Elevation values are indicated in metres, while the values of the different landcovers represent percentages. Table S2. Description of the best model for each bear class in each season. *Indicates that there was more than one model within 2 units of DAIC from the best one and thus the selected model was the simplest one (fewer parameters, and if the same, a smaller number of features and regularization multiplier). fcden, female with cubs during denning; fchyp, females with cubs during hyperphagia; Fcmat, females with cubs during mating; fyden, females with yearlings during denning; fymat, females with yearlings during mating; oden, other bears during denning; ohyp, other bears during hyperphagia; omat, other bears during mating; rm, regularization multiplier. Model features: linear (L), hinge (H), quadratic (Q), threshold (T) and product (P). Table S3. Mean habitat suitability values predicted by the models of each bear class in each season and for all bear occurrence points. fcden, female with cubs during denning; fchyp, females with cubs during hyperphagia; Fcmat, females with cubs during mating; fyden, females with yearlings during denning; fymat, females with yearlings during mating; oden, other bears during denning; ohyp, other bears during hyperphagia; omat, other bears during mating. Table S4. Percent contribution of each of the variables to the MaxEnt models (r: including distance to the nearest road, nr: not including this distance) for each bear class in each season. The variable distance to the nearest road is highlighted, and the model in which it contributed the most (females with yearlings in denning) is in bold type. fcden, female with cubs during denning; fchyp, females with cubs during hyperphagia; Fcmat, females with cubs during mating; fyden, females with yearlings during denning; fymat, females with yearlings during mating; oden, other bears during denning; ohyp, other bears during hyperphagia; omat, other bears during mating. Table S5. Percentage of the number of cells crossed by a road with more than a 70% likelihood of bear occurrence (cloglog value > 0.7), and thus high habitat suitability, for each bear class in each season. fcden, female with cubs during denning; fchyp, females with cubs during hyperphagia; Fcmat, females with cubs during mating; fyden, females with yearlings during denning; fymat, females with yearlings during mating; oden, other bears during denning; ohyp, other bears during hyperphagia; omat, other bears during mating.