Prioritizing road‐kill mitigation areas: A spatially explicit national‐scale model for an elusive carnivore

Roads impact wildlife in different ways, among which road mortality has been the most studied. Budgets in conservation biology are usually small, and macroecological approaches have been employed in recent years as the first steps towards guiding management. Carnivores are particularly vulnerable to mortality on roads due to their elevated ecological needs (low population density, often low fecundity and relatively large home ranges). Our aim was to develop a ranking methodology to prioritize specific areas for road‐kill mitigation.

Budgets for conservation are usually constrained. In recent years, macroecological approaches that have been frequently based on open-source data, such as citizen science projects, have been increasingly employed to guide managers before implementing expensive mitigation actions. For instance, Brehme, Hathaway, and Fisher (2018) used macroecological approaches to rank the reptiles and amphibians most susceptible to road mortality and fragmentation in California based on their life histories and space-use characteristics after a literature review. In a spatially explicit approach, González-Suárez, Zanchetta Ferreira, and Grilo (2018) explored how life history traits could explain road-kill risk for birds and mammals in Brazil. On the other hand, Beston, Diffendorfer, Loss, and Johnson (2016) developed a prioritization system to identify the avian species most likely to experience population declines from wind facilities based on their current conservation statuses or population sizes and their expected risk from fatalities due to collision with turbines.
Other authors like D' Amico et al. (2019) used atlas data and species traits from the literature to rank species based on their risk of population extinction from collision with power lines and to identify the areas where high richness of these species overlapped with current electricity grids. In all cases, these approaches were the first step before implementing field-level mitigation actions, such as fencing, underpass construction or wire marking.
Here, we present a method aimed at identifying areas with greater concentrations of road sections with high risk for road-kills at a large scale, with the elusive polecat (Mustela putorius) as the case study and continental Italy as the study area. For this declining species (Croose et al., 2018), road mortality represents one of its most important threats (Barrientos & Bolonio, 2009;Barrientos & Miranda, 2012;Croose et al., 2018). Specifically, we aimed to (a) identify the variables, both habitat-and road-related, determining the polecat road mortality in Italy; (b) generate a spatially explicit road-kill risk map for the polecat at the national level; and (c) identify the top-ranking cells concentrating road sections where mitigation to avoid road-kills should be prioritized.

| Summary of the modelling framework
The map of road risk for polecats was built by combining several raster maps, each describing a specific factor that could affect polecat collision risk for this species. These included road characteristics and ecologically related factors like habitat suitability and prey availability and five environmental factors related to the matrix permeability to polecat movements ( Figure 1). Habitat

K E Y W O R D S
citizen science, Italy, Mustela putorius, prioritization, Road Ecology, road-kill, species distribution model suitability was derived from a spatial distribution model that combined the polecat occurrence data with a set of environmental variables related to polecat biology. Prey availability was derived from the combination of available expert-based suitability maps of 26 preferred prey species. The relationships among the road risk locations and the eight variables were analysed with generalized linear models (GLM). The variables included in the best performing model were then combined to produce the road risk map and to identify the areas with the highest densities of risky road sections (see Figure 1 for flow chart).

| Road-kill and presence data
The locations of road-kill incidents (n = 212) and live individuals (n = 271) ( Figure 2) were obtained from regional authorities, literature (Bizzarri, Lacrimini, & Ragni, 2010;Bon, 2017;Fusillo & Marcelli, 2014;Ragni et al., 2014;Rondinini, Ercoli, & Boitani, 2006), researchers, NGOs and from those citizen science projects (Appendix S1) providing data validated by experts. Moreover, we used only data that referred to direct observations and discarded records of indirect signs, such as scat or footprints. Previous studies have shown F I G U R E 1 This flow chart summarizes the various procedures used to create the road risk model for the polecat in Italy. CLC, Corine Land Cover; GLM, generalized linear model; NEN, National Ecological Network; SDM, species distribution model; DEM, digital elevation model F I G U R E 2 Locations of the polecat occurrences (black triangles; n = 271) used to generate the SDM, and the polecat road-kill events (white circles; n = 212) used to produce the road risk model. Data are mapped on the Italian range of the species based on IUCN (Skumatov et al., 2016) the contribution of citizen science to a number of research topics, including the study of road fatalities (Barrientos & Miranda, 2012;Heigl, Horvath, Laaha, & Zaller, 2017;Périquet, Roxburgh, le Roux, & Collinson, 2018;Tiago, Pereira, & Capinha, 2017). The road-kill dataset included 60 records in which the specimen was found on the road, although it was not explicitly reported as a road casualty. The dataset covered the time spans of 1990 to 2017 for road-kills and 1968 to 2017 for direct observations, considering that the road network has undergone few changes in these years (Pinto & Franchin, 2010). However, the oldest data represented only a small proportion of the dataset (road-kills: 1990-1999 = 7%; direct observations: 1968-1999 = 8%).

| Habitat suitability
We used the 271 live polecat records to build a species distribution model ("SDM", hereafter) for the polecat in Italy. As there is a strong correlation between the road-kill numbers and the species ). Specifically, we included the following seven algorithms into the ensemble modelling approach: general linear models, generalized additive models, generalized boosted models, classification t analysis, artificial neural network, multiple adaptive regression splines and random forest models (see Appendix S1 for the model calibration).
We randomly splitted each the occurrence dataset into an 80% sample for model calibration and 20% for model validation; the procedure was repeated ten times, and the results were averaged (Russo et al., 2014). The predictive performance of the models was assessed by measuring the area under the receiver operating characteristic curve (AUC; Hanley & McNeil, 1982)

| Prey availability
Because "prey availability" may represent a variable associated with road-kills (Barrientos & Bolonio, 2009;Barrientos & Miranda, 2012), we considered the habitat suitability models provided by the National Ecological Network (NEN) project (Boitani, Falcucci, Maiorano, & Montemaggiori, 2003) for 26 potential polecat prey species (Baghli, Walzberg, & Verhagen, 2005;Lodé, 1997;Prigioni & De Marinis, 1995;Sidorovich & Pikulik, 1997; see Appendix S1). The NEN rasters are provided at a resolution of 100 × 100 m grid cells, where each cell has a suitability score from 0 (unsuitable) to 3 (highly suitable). The 26 suitability maps were combined into a single raster by summing the suitability score of each cell. This map was then used as a proxy for prey availability in a 50 m buffer radius generated around the random and road-kill points.

| Road risk model
To produce the road risk model (RRM), we related a set of environmental characteristics with 212 road-kill locations and an equal number of randomly chosen points along the road network. One random point per road-kill was generated within a section of road included within a circular buffer, with a radius of 1.5 km, centred on the corresponding road-kill. We selected this radius because it is the average home range size for polecats (Baghli et al., 2005;Rondinini et al., 2006) and therefore represents a conservative approach for selecting control points within the potential polecat habitat. Within each buffer, we recorded the SDM value and the prey availability value. Also, we used "type of road" as a categorical variable because it is related to road-kill rates (Barrientos & Miranda, 2012). It had three levels, that is, "state," "provincial" and "local" roads. State roads are the main traffic routes, have speed limits of 90-100 km/hr and are often bounded by barriers (e.g. Jersey barriers). Provincial roads connect various cities in the same region and have speed limits of approximately 70 km/hr. Finally, local roads are urban roads with speed limits of 50 km/hr.
All types of roads are usually two-lane roads. Motorways were not considered due to their lack of representation in the data (only one random point and three road-killed polecats were found along motorways). Since urban habitats are positively related to carnivore road casualties (Barrientos & Bolonio, 2009;Grilo et al., 2009) and can be positively selected by polecats (Baghli et al., 2005;Virgós, 2003;Zabala, Zuberogoitia, & Martínez-Climent, 2005), we calculated "human influence" as the distance in metres to the nearest town and divided it by the population density (ISTAT, 2011); accordingly, surroundings of more populated towns have greater weight than less populated ones. As polecats are known to move across specific habitat types (Baghli et al., 2005), we used the most relevant land cover category within a buffer of 50 m radius around each point to determine the coverage of these preferred habitats, that is, "grasslands" (CLC 231), "pastures" (CLC 321) and "broadleaved forests" (CLC 311), as a proxy of the landscape structure (Barrientos & Bolonio, 2009). In case the buffer did not fall into one of these categories, the value zero was assigned. We used the 1990 CLC (scale 1:100,000) inventory for road-kills from 1990 to 1999 and the 2000 CLC inventory for road-kills from 2000 to 2017. Finally, as polecats are also known to move preferentially along shores (Mestre, Ferreira, & Mira, 2007;Rondinini et al., 2006;Zabala et al., 2005), we measured the distance in metres to the nearest lake or river (hereafter, "water bodies"). All of these spatial analyses were performed using the program QGIS (QGIS Development Team, 2017).

| Statistical analyses
We conducted a regression analysis on the set of eight environmental variables using a binomial distribution (road-kill points vs. random points) and logit link function. We used the best subset procedure and the Akaike information criterion (AIC) to identify the set of models (i.e. combination of variables) that best explained the occurrence of road-kills. We used this technique because it yields consistent results regardless of the order in which variables are included in the model and allows models with different numbers of parameters to be directly compared with each other (Burnham & Anderson, 2002). We used AIC values without correction for small sample size because the number of observations was high (i.e. the ratio between the number of observations and explanatory variables was over 40) for the number of explanatory variables (Burnham & Anderson, 2002). The models with the lowest AIC values represent the best compromise between the maximal fit and the minimal number of explanatory variables (i.e. statistical parsimony). To evaluate the relative explanatory power of the competing best models, the Akaike weights (ω i ) were calculated. The evidence ratio was calculated to compare the Akaike weights of the best model and the competing ones (Burnham & Anderson, 2002). To estimate the relative importance of every variable included in any of the best models (those with ΔAIC ≤ 2), we calculated the sum of the Akaike weights of the models that included these variables (Burnham & Anderson, 2002). The quantitative differences in the explanatory variables between stretches with and without road-kills were evaluated with paired t tests. We analysed the differences in the proportion of road-kills per type of road by using 2 × 2 contingency tables with the Yates correction. All analyses were performed with the Statistica v.10 software (StatSoft).
We first obtained a subset of the models that best separated the road-kills from the random points (Table 1). We assessed the predictive performance of the two models with best AIC values with the same procedure used for SDM (see above). We then compared the AUC values of these two models using the Wilcoxon-Mann-Whitney and selected the model with highest values of the AUC.
We produced a spatially explicit prediction of the road risk at a 100 × 100 m resolution. Finally, we overlapped the Italian Universal Transverse Mercator 10 × 10 km grid with the RRM and ranked the cells based on their summed road risk values.
The two top-ranked models in Table 1

| D ISCUSS I ON
The approach used in the present study enabled for the first time testing of how environmental variables affected road-kill risk at a large scale, which led to a ranking of the areas with high numbers of risky road sections based on relatively easily accessible information. This method was empirically supported as it was based on the traits that differentiated the road sections with recorded road-kills from those without. Since the mitigation actions to prevent roadkills, such as wildlife passage construction or fencing, are expensive (Rytwinski et al., 2016;Smith, van der Ree, & Rosell, 2015;Van der Ree et al., 2015), our macroecological approach represents a promising tool for selecting road sections where these mitigation actions should be implemented, thus reducing costs (e.g. generalized mitigation). The 10 × 10 km cells were ranked by their road-kill risk, which was determined by the total length of risky road stretches within the cells, that is by the abundance of broad-leaved forests bisected by national or provincial roads within suitable polecat distribution range.
Once focused on those cells with higher concentrations of potentially high-risk road sections, the mitigation actions should be prioritized in the road sections with the highest potential risk (those marked in red colours in the example from Figure 4). We would like to add a cautionary note to our model as some landscape traits can be captured with a low resolution by a 50 m buffer. Our approach will become more accurate in the future if fine-scale environmental data become available. However, we chose to use this buffer size as it is the most accurate for achieving truly effective mitigation actions at local scales (Barrientos & Bolonio, 2009;Barrientos & Miranda, 2012).
The importance of broad-leaved forests seems to be related to landscape structure and animal movements. This kind of forest TA B L E 1 List of the twelve models that best separated the road-kills from the random points. The ΔAIC was the difference in the AIC values compared to the estimated best model (lowest AIC) that allowed the ranking of models from an estimated best (top of the table) to worst. Akaike weight is the estimated probability that a model is the best model in the set. The evidence ratio indicates to what extent the best model is better than another. The K is the number of predictors contained in the model. See the full set of models in Appendix S1

Model no.
Variables contained in the model K ΔAIC Akaike weight (ω i ) Evidence ratio includes riparian and gallery woodlands. These ecosystems were selected in other areas by polecats, such as in Luxemburg, Spain, Italy or Portugal, probably because the dense stream margin vegetation provides protection during their movements (Baghli et al., 2005;Mestre et al., 2007;Rondinini et al., 2006;Zabala et al., 2005). This is a general pattern for all carnivores, especially in agriculture-dominated matrixes (Virgós, 2001). Riparian habitats also provide a variety of food resources because they act as ecotones, which increases their value for polecats as hunting areas (Mestre et al., 2007;Zabala et al., 2005). Beech and oak forests were also used at above-average rates in the study of Baghli et al. (2005), especially in summer. Although we do not know the sex of the casualties in our study, road-kill datasets are usually male-skewed (Barrientos, 2015). Interestingly, the preference to move across deciduous forests seems stronger in males (Baghli et al., 2005), likely because males usually move for longer distances in search of receptive females and are more prone to dispersal (Rondinini et al., 2006). A significant relationship between the type of road (highly related to traffic flow) and road-kill rates has been frequently found, including in studies on carnivores (Clarke, White, & Harris, 1998;Philcox, Grogan, & Macdonald, 1999;Seiler, 2005). Namely, roads with low traffic levels (i.e. local roads) do not represent a threat compared to other types of roads, likely because they usually have little traffic that moves at low speed; on the contrary, roads with both higher traffic and speed limits (i.e. provincial and national roads) present higher road-kill rates. Finally, the high traffic intensity on motorways discourages wildlife from attempting to cross them, although we do not have enough data to validate this last conclusion. However, it is also possible that the population F I G U R E 3 Road-kill risk map for the polecat in Italy. Cells (10 × 10 km) were ranked based on their summed road risk F I G U R E 4 Enlargement of one of the cells with the highest concentrations of high-risk road sections (in red) for road mortality of polecat in Italy declined in these areas due to previous road mortality (Ascensão et al., 2019;Teixeira, Kindel, Hartz, Mitchell, & Fahrig, 2017).
The importance of potential prey did not determine the occurrence of road-kills in Italy, as happens, for instance, in Mediterranean Spain (Barrientos & Bolonio, 2009;Barrientos & Miranda, 2012).

This could be because a single prey species (the European rabbit
Oryctolagus cuniculus) is by far the main prey in Spain (Barrientos & Bolonio, 2009), which makes it easier to model the influence of prey abundance on the probability of road-kill. Alternatively, the field sampling carried out by Barrientos and Bolonio (2009) could be more accurate than the map of habitat suitability for prey species (Boitani et al., 2003) that we used in the present study.
Finally, our study emphasized the potential utility of citizen sci-  (Périquet et al., 2018). Additionally, small carcasses such as those of mustelids persist on the road for a shorter period and are detected at lower rates compared to those of large mammals . Those carcasses that are repeatedly run over are flattened, becoming unattractive for citizen scientists, who lose motivation or put less effort into reporting unidentified carcasses or non-flagship species (Lukyanenko, Parsons, & Wiersma, 2016). These biases might have affected the unequal distribution of records of live polecats in the study area, that is the small number of records collected in southern Italy, while the number of road-kill incidents remained high in this region. This in turn might explain the low predictive performance of the SDM in our study. However, according to Schwartz, Shilling, and Perkins (2020) and to our experience with other mustelids (Fabrizio, Di Febbraro, D'Amico, et al., 2019;, this claims for increasing efforts in collecting road-kill data on small-medium mustelids at country scale. We stress the relevance of considering both the RRM and the potential distribution of the target species as the most appropriate method for prioritizing road sections in need of mitigation measures. Despite the utility of citizen science to identify priority areas for mitigation, we should not lose sight of the fact that the final step of any road-kill study should be to apply the corresponding mitigation measures in the field. The implementation, and success in reducing road-kills, of these measures (fencing, construction of wildlife passages) will be the definitive proof that our method has worked.
Finally, medium-to long-term monitoring is necessary for confirming that mitigation has been effective in reducing road-kills, as management actions need time for animals to adapt to them (e.g. Corlatti, Hackländer, & Frey-Roos, 2009; Soanes et al., 2013). This monitoring should be focused on population-level impacts and include fieldwork and genetic analyses (Corlatti et al., 2009).

DATA AVA I L A B I L I T Y S TAT E M E N T
Data can be downloaded in most of the cases from those databases indicated in Appendix S1. In some cases, data cannot be shared because they have been obtained upon agreement from data providers. However, in these latter cases, researchers can get in contact with the persons listed in the Appendix S1 to ask permission to access these data.