Influence of Landscape Factors on Amphibian Roadkills at the National Level

Roads exert multiple effects on wildlife, from animal mortality, habitat and population fragmentation, to modification of animal reproductive behaviour. Amphibians are the most frequently road-killed animal group. Many studies have attempted to analyse the factors driving amphibian casualties on roads, but these factors are limited to the roads themselves (e.g., traffic, type of roads, roads crossing water bodies) or to structures along them (e.g., ditches, walls). Sometimes, roadkills are related to land use along the roads. We analysed the influence of landscape factors on roadkill hotspots at the national level (Slovenia). Specifically, we aimed at: (1) identifying hotspots of roadkills, (2) analysing whether records of amphibian presence on roads are related to the distribution of water bodies and (3) analysing which factors (proximity to water bodies or human factors) explain the distribution of hotspots. Hotspots were identified by Malo’s method. Roadkills were modelled with Maxent for the first time in Slovenia. The relationships between amphibian presence and hotspots with factors were analysed with GLM. A total of 237 road sections were identified as hotspots, corresponding to 8% of road sections and containing 90% of road-killed individuals. Proximity to forests, meadows and agricultural land were the most important variables in Maxent models. The number of roadkills depended on the proximity to agricultural land, forests, water bodies and wetland areas, while the number of hotspots additionally depended on the proximity to urban settlements.


Introduction
Roads have strong negative effects on biodiversity [1] and animal roadkills are one of the most crucial mortality factors caused directly by humans [2,3]. Amphibians are the most frequently road-killed group, with casualty rates of 60-90% in vertebrate road ecology studies [4,5]. Roadkills are determined by many factors. Some factors are driven by road characteristics (type of road, number of cars per time unit) and also by driver behaviour [6][7][8][9][10][11]. Indeed, studies have shown that roadkills of snakes are intentional [12][13][14]. Even the level of conservation protection affects the risk of being road-killed: protected areas can have a higher proportion of roadkills [15]. Furthermore, roadless areas have the highest conservation capacity [16][17][18][19]. Factors related to landscape structure and composition determine the presence of animals on roads and, therefore, the probability of being killed by a vehicle [10,[20][21][22]. In the particular case of amphibians, roadkills are associated mainly with traffic intensity, quality of adjacent reproduction places, dispersion patterns and capacity, and activity patterns [7,10,[23][24][25]. Amphibians depend on a complex landscape structure in order to complete their life circle, moving between reproductive and feeding habitats [11]. Roads act as barriers [21], increasing the impermeability among habitats and mortality [20,21]. Road fragmentation of habitats consequently produces fragmentation of populations and modification of animal behaviour [26][27][28].
Many conservation measures have been used to mitigate roadkills [29][30][31]. As resources are limited, mitigation measures must be applied at least in places with higher roadkill rates, the so-called hotspots [11,21]. Hotspots can be identified with spatial statistics and modelling procedures, such as probabilities based on Poisson distribution [32], binary logistic regression, ecological niche factor analysis, kernel density estimation, nearest-neighbour hierarchical clustering [33], Ripley's K function [8,34,35] or Getis-Ord Gi* [36]. Hotspots depend on: (1) points of animal crossings; (2) barriers for crossing roads safely, where the animal is blocked in the road; and (3) high roadkill rates [11].
In Slovenia, studies on amphibian road mortality have been carried out at the local, regional and national levels in the last 20 years [37][38][39][40]. So far, there are more than 1500 locations of roadkills registered in the country [41]. In addition, some long-term volunteer conservation actions have been conducted for several years [40,42,43], with new actions starting every year. Despite all conservation efforts, only few amphibian mitigation measures have been implemented on state roads [40]. In light of this, identification of hotspots is a valuable tool for national planning and prioritising locations for road mitigation measures on state roads in the future.
We analysed the influence of human or landscape factors on roadkill hotspots at the national level (Slovenia). Specifically, we aimed at: (1) identifying hotspots of roadkills; (2) analysing whether records of amphibian presences on roads are related to the distribution of water bodies; and (3) analysing which factors (proximity to water bodies or human factors) explain the distribution of hotspots.

Study Area
Slovenia (20,272 km 2 ) lies at the junction of four main biogeographical regions-Alpine, Dinaric, Mediterranean and Pannonian ( Figure 1). Thus, Slovenia's nature is influenced by all these regions and the transitions between them, which are also reflected in its herpetofauna [44]. With 20 reported species of native amphibians, Slovenia ranks among the richest Palearctic countries in terms of its amphibian diversity [45], especially considering its small size.
Owing to the mountainous character of the Alpine and Dinaric regions, Slovenia's average elevation is slightly above 556 m [46]. The Alpine region in the northern and western parts of the country (42.1% of its total area) includes the highest mountains and the highest-lying valleys of the major rivers, while high plateaus, deep and narrow valleys and flat karst poljes (karst plains) characterize the physical relief of the Dinaric region (28.1%) in the south [47]. With their intermittent rivers and lakes, these poljes are the only surface wetlands in the major part of this region. Towards the east of the country, the average elevation drops as the mountainous terrain turns into the Pannonian lowland. This region (21.2%) comprises the lower reaches of major rivers that dissect wide plains and smaller hilly landscapes. As this region also contains numerous fishponds, gravel pits, reservoirs and other artificial bodies of water, it is the richest in terms of surface wetlands. The average elevation also drops in the southwestern part of the country towards the Adriatic Sea, where the Mediterranean region starts, characterized by plains and low hills. Mediterranean climate prevails along the coast, but its impact reaches further inland. This region has higher temperatures and is under a strong Mediterranean influence. The rest of Slovenia has mild continental climate, while in the Alps and in the higher parts of the Dinaric Mts, mountain climate prevails. Local climate is significantly influenced by altitude and ragged relief [44,48].
Throughout history, human activities have played a very important role in terms of amphibian distribution and richness. From 1970 onwards, a significant decline in the diversity of landscapes and habitats has been witnessed. This is mostly the result of urbanization and intensive agriculture. Consequently, the fragmentation of suitable habitats for amphibians has increased, raising the risk of extinction, especially for less common amphibian species [44]. More than half of Slovenia's land area is covered by forests (56% or 58%, including shrubland), while other mainly natural areas, natural grasslands, wetlands, water bodies and open spaces with little or no vegetation, cover 4% of the land. Thirty-five percent of the surface area is intended for farming, while slightly less than 3% is built-up area [49]. Many aquatic ecosystems are degraded due to an uncontrolled and dispersed urbanization pattern. Slovenia has c. 2 million inhabitants and more than half of the country's population live in settlements with less than 2000 inhabitants. This dispersed and uncontrolled pattern of settlement is also the reason for urbanization of floodplains and marshlands [44], which has led to the disappearance of many amphibian habitats in the Slovenian lowlands. Slovenia has almost 39,000 kilometres of public roads and tracks, including 6707 km of state roads, managed by the Slovenian Infrastructure Agency [50], that are included in this study. Especially after joining the European Union, road transport in Slovenia has increased threefold. The growth can be ascribed to the increase in domestic transport as well as to the transit arising from the geographical position of Slovenia. As a result, roads have been recognised as the factor with a great negative effect on amphibian populations in Slovenia [41].

Amphibian Data
The Centre for Cartography of Fauna and Flora (CKFF) is maintaining the national amphibian database for Slovenia. Collected data are the results of more than 20 years of the Centre's amphibian conservation projects in cooperation with Societas herpetologica slovenica.
In 2018 [40], CKFF collected all available data on amphibian migration across the state-managed roads in Slovenia, which are presented as follows. More than 90% of presence-only data of amphibians on roads (dead or alive) are the results of many projects, where extensive road surveys have taken place in the last 15 years [40]. The survey method in all projects was the standard night driving [51] in suitable weather conditions, while the amphibian sightings have been recorded on state roads. There are at least three projects where surveys have been conducted on a large national or regional scale. The first national overview was prepared in 2000 [37]; at that time, the amphibian road mortality was recorded on 7.4% (478 km) of the state-managed roads. In 2003 followed a cross-border project, "Amphibians and Traffic in the Alps-Adriatic Region" [38], and in the 2005-2007 period the Interreg IIIA Slovenia-Austria project "Amphibian and Bat Conservation in the Alps-Adriatic Region" [39], when more than 2500 km of roads were intensely surveyed in the northern part of Slovenia. Additionally, we included data from 13 road sections where volunteer conservation actions took place and temporary fences were used [40,42,43].
We also collected publicly available data from citizen science projects (published and web resources, ca. 5%) as well as random observations by direct contacts (ca. 5%) with several individuals, societies and nature conservation organisations, all dealing with amphibians and road problem in Slovenia. This is the result of growing amphibian conservation awareness in society, which has also significantly increased the number of volunteer conservation actions on Slovenian roads in the last 10 years. All species data were verified by amphibian specialist (by photo or trusted source), and only spatially accurate data (GPS-obtained data and/or exact plotted data on the maps and/or data with exact description, sometimes with official state road markings) were used in the analysis.
In the analysis we included two datasets: (a) roadkills: amphibian presence-only data on roads (alive or dead; 2991 species records, in total 166,906 individuals from 1691 locations) and (b) amphibian sites: amphibian presence-only data in neighbouring habitats within 1 km buffer area from the road (16,291 species records from 7520 locations) ( Figure 2). Species record is observation of the number of individuals of a single species on one location on a given date. In case of several records for the same location, we selected records with maximal abundance for single species. We extracted data for amphibian sites from the national database. Amphibians were surveyed in accordance with a licence of the Centre for Cartography of Fauna and Flora (CKFF) issued by The Ministry of Environment and Spatial Planning (35601-121/ 2008-4) Figure 2. Slovenian state roads network (6707 km) [52] and amphibian locations (7520) in 1 km buffer zone.

State Roads and Land Use
The analysis included state roads managed by the Slovenian Infrastructure Agency (SIA), that is, motorways/expressways, main roads and regional roads, in total length of 6707 km. Vector road data were provided by SIA [52], which also included traffic information-average daily number of vehicles per road section. Traffic was used in GLM analysis.
Land use data for the whole of Slovenia are freely available from the Ministry of Agriculture, Forestry and Food [53] and were used to relate roadkills and hotspots with landscape factors. Agricultural and Forestry Land Use Database is an official national database of agricultural and forestry land use in Slovenia. The basic element in the database is a polygon that represents the unique part of land with the same land use, which is captured from aerial imagery (with a scale 1:5000) and additionally verified by field-gathered data and other possible provided information. It is regularly updated and we used the last available version from 2017. Data are aggregated in the following classes: Agriculture, Forest, Meadows, Urban, Water and Wetlands. For the analysis, we used land use data within 1 km buffer area on the state roads. In order to avoid categorical variables, the distance of each land use class to the closest road section was calculated.

Determination of Roadkill Hotspots
We divided roads in 1 km sections, using the GRASS function "v.split.length" in QGIS. As the road vector data was composed of different features (line elements), some sections overlapped at road intersections. Thus, the number of road sections was higher than the number of total road kilometres. We counted the number of roadkills in each section (but only once when several sections overlapped), independently of time. We identified roadkill hotspots (sections with a higher number of roadkills than expected by chance) following the Malo's method described in [32], which detects spatial clusters of roadkill records. Malo's method [32] compares the frequency distribution of roadkill records with a null Poisson distribution. The number of GPS points belonging to a particular road section was counted with the command "Count points in polygons" in QGIS. The spatial pattern of roadkills was compared with a random model where the likelihood of roadkills for each road section follows the Poisson distribution. Using this hypothesis and given the mean number of roadkills per road section (λ), the probability of any road section having x number of roadkills is: In a random distribution of road-killed specimens, the hotspots of roadkills should be distributed at random across the roads, and their aggregation would be extremely unlikely. Hotspots of roadkills were identified using a graph with the probability values (p(x)) versus the number of roadkills (x). These sections were located in the graph zone were the curve reaches the asymptote. Thus, all roadkill records above the point where the curve reaches the asymptote are considered hotspots, as their probability of occurring is very low. Hotspot records constitute the Hotspots database. In order to analyse only the most important hotspots, we selected the top 5% and top 10% of all hotspots with respect to the total number of roadkills.

Ecological Niche Models
The ecological niche model was used [54,55] to analyse the relationship between Roadkills, Hotspots and Amphibian sites datasets and land use data. As absences of roadkills cannot be guaranteed, we used a correlative realised niche algorithm [54] for the presence-only records: Maximum Entropy (Maxent), a general purpose machine-learning method that uses presence-only occurrence data and background data [56][57][58]. Maxent is particularly well suited to noisy or sparse information and is capable of simultaneously dealing with continuous and categorical variables [56][57][58]. It looks for the statistical model with the most uniform distribution but still infers as accurately as possible the observed data, selecting at random uniformly distributed data from the background pixels.
Here, background sample does not mean species absence at the selected sites, but rather provides a spectrum of the available conditions [59]. In fact, [59] obtained better results extracting data from those pixels of the background close to the species' presences. Maxent starts with a uniform probability distribution (gain = 0) and iteratively alters one weight at a time to maximise the likelihood of the occurrence data set. Since the algorithm converges to the optimum probable distribution, the gain can be interpreted as a representation of how much better the distribution fits the sample points than the uniform distribution does [56][57][58]. Maxent output represents the habitat suitability, ranging from 0.0 to 1.0.
As Maxent is a probabilistic modelling method, we calculated the arithmetic mean and the standard deviation of a set of 100 models per species and dataset through an iterative process [60]. We chose 100 models as a compromise among statistical analysis power, computation time and physical storage. We ran Maxent in clog-log format with default parameters using all presence records as training data [58]. Duplicated records (i.e., records inside the same pixel) were removed. The realised niche models were performed with Maxent 3. 4.1 software [61].
The Maxent models were tested with receiver operating characteristic (ROC) plots. The area under the curve (AUC) of the ROC plot was taken as a measure of the overall fit of the models [62]. Random models have an AUC equal to 0.5; the closer an AUC is to 1, the better the model is. AUC was selected because it is independent of prevalence (i.e., the proportion of presence in relation to the total dataset size [63]); but see [64]. In addition, we calculated a set of 100 null models, following the methodology outlined by [65]. We used Maxent models to identify the importance of each independent variable in explaining species distribution through factor analysis: (1) a jackknife analysis of the average AUC using training and test data and (2) a calculation of the average percentage contribution of each variable to the models. For this purpose, the variables were excluded in turn and a model was created with the remaining variables; then, a model was created using each individual variable.

GLM for Number of Roadkills
We analysed the relationship among the number of roadkills and hotspots per road section with traffic and land use data using Generalised Linear Models (GLM) with a Negative Binomial distribution as the response variable is discrete. We selected randomly the same number sections without roadkill events than roadkills and hotspots in order to avoid zero-inflated problems. The process was repeated 10 times, obtaining therefore 10 different datasets, each with different sections without roadkills and the same sections with roadkills. GLM is a flexible generalisation of ordinary linear regression that allows for response variables that have error distribution models other than a normal distribution. GLM allows the linear model to be related to the response variable via a link function (negative binomial distribution) and the magnitude of the variance of each measurement to be a function of its predicted value. Number of roadkills and hotspots were correlated with each independent variable previous to the GLM analyses. All statistical analyses were performed using R software.

Results
We collected 2991 roadkill observations, corresponding to 167,417 roadkill individuals on the roads analysed in 9642 1 km road sections in Slovenia: 1782 with fatalities and 7860 without any roadkill records.

Determination of Roadkill Hotspots
The graph obtained by Malo's method identified as hotspots the sections with ≥34 roadkills (Figure 3), resulting in 237 sections (7.92% of total road sections). This corresponds to 149,989 roadkills (893 roadkills points, 89.59 % of total road-killed individuals). We also identified the top 5% and 10% of hotspots in terms of total roadkills (Figure 4), corresponding respectively to 24 sections (≥1300 road-killed individuals) and 12 sections (≥2800 roadkill individuals). Hotspots are widespread in Slovenia but are more abundant in the eastern part of the country (Figure 4).

Ecological Niche Models
All independent variables had a pairwise correlation lower than 0.7, and were included in the models. All Training and Test models obtained an AUC higher than 0.7 (Table 1). AUC Null models were significantly lower than AUC from Roadkills, Hotspots and Amphibian Sites models (Kruskal-Wallis chi-squared = 149.26, df = 1, p-value < 0.0001; Null model: 0.66 ± 0.03; Roadkills: 0.84 ± 0.01; Hotspots: 0.90 ± 0.01; Localities: 0.79 ± 0.01). Figure 5 shows the spatial output of Maxent models for the three datasets. Proximity to meadows was the most important variable in the Hotspots and Roadkills models, while proximity to forest was the most important variable in the Amphibian Sites model (Table 1). For the Hotspots and Amphibian Sites models, the second most important variable was proximity to agriculture; for the Roadkills model, this variable was proximity to forest (Table 1).

Generalised Linear Models
All Roadkills and Hotspots GLMs showed similar patterns of low correlation between the number of roadkills and hotspots per section and the independent variables. All correlations were not significant ( Table 2). Proximity to agriculture, meadows and urban areas were not significant variables in the GLM for Roadkills dataset (Table 3; we show only one of the models). Traffic and proximity to wetlands presented negative trends, while proximity to forest and water presented positive trends. In the case of Hotspots models, traffic was the only significant variable with a negative trend. Proximity to water presented a marginal significance (Table 3). Table 3. Generalized Linear Model (GLM) results for Roadkills and Hotspots datasets. Only one of the 10 models performed using 10 random absence datasets are shown here. Significance levels are indicated with an *: p < 0.05; ***: p < 0.001.

Determination of Roadkill Hotspots
We successfully used Malo's method [32] to identify hotspots of amphibian roadkills in Slovenia. Malo's method has been suggested to be the best one for hotspot identification [33], in comparison with other statistical methods such as binary logistic regression (BLR), ecological niche factor analysis (ENFA), kernel density estimation and nearest-neighbor hierarchical clustering (NNHC). K-function [32,34,35] and Getis-Ord Gi* [36] have also been used as spatial clustering methods for roadkills.
We identified 237 sections as hotspots, corresponding to~8% of road sections and 90% of road-killed individuals included in the analyses. All of the most species-diverse road sections corresponded to the top 5% of all hotspots, but there are top 5% and 10% hotspots that did not correspond to the most diverse road sections. Some top hotspots corresponded to sections with only one species. Although hotspots are widespread, they are more common in eastern Slovenia. However, this can be an artefact due to bias in the sampling effort. Unfortunately, logistic limitations hamper a homogeneous effort across the entire country. The result is higher than the number provided by the national overview of amphibian road mortality on state roads [40], where hotspots were registered on 17.7% of roads using empirical methods. The number of hotspots was considerably reduced when selecting the top 5% and 10% with a higher number of roadkills. However, they were all identified also in the top 10% road sections in the empirical study [40]. Hotspots were concentrated in central, eastern and western Slovenia. Other studies provided higher values of number of hotspots. For example, higher proportions of hotspots were found in Hoces del Alto Ebro y Rudrón Natural Park (Spain), using Kernel density functions [66]: 35.12% of road sections (69.1% of all roadkills of amphibians and reptiles). However, other studies also presented low numbers of hotspots. In north Portugal, using Malo's method as well, 8.72% of road sections (35.76% of amphibian roadkills) were considered as hotspots [11]. In Salamanca (Spain), hotspots corresponded to 3.77% of road sections and 56.09% of amphibian roadkills [21]. Regarding amphibians and birds in Thousand Islands Parkway, Ontario, Canada, 79% of roadkills were found in two out of six road sections [67]. Only three road sections in Romania were considered hotspots for tortoises in almost 4000 km of roads [68]. In California and Maine, using Getis-Ord Gi*, 10% of roads were identified as hotspots for mammals [36].
Only 8% of road sections contain 90% of road-killed individuals found in Slovenia; this value, however, is unsustainable for amphibian populations [69]. Slovenia is considered a well-conserved country, with large extensions of forests. This may be the reason for having such high numbers of roadkill hotspots, as animals are prone to be road-killed with higher intensity in preserved regions [15].

Ecological Niche Modelling of Roadkills and Hotspots
This is the first time that Maxent was used to model roadkills in Slovenia. Maxent was used to model roadkills of mammals and birds in California [22], of the Iberian Linx in Spain [70] and leopards in Iran [71], and as a permeability layer for road connectivity studies [20]. Other methods of ecological niche models [54], such as logistic regressions, have been used to model roadkills of the Moose in Sweden [72], of the Eurasian Lynx in Switzerland [73], of the Grizzly Bear in Canada [74] or of several marsupials in Australia [34]. For example, binary logistic regression and ENFA modelling have been used for recognizing roadkill deterministic factors [33].
The proximities to forests, meadows and agricultural land were the most important variables in Maxent models of Roadkills, Hotspots and Amphibian sites datasets. Effectively, amphibians depend on water bodies that are very common in forests and agricultural areas. Other studies with amphibians have obtained similar results: anuran roadkills are determined by proximity to water bodies [10,25] and to unaltered habitats [11,21]. Roadkills are located mainly along migratory paths [21] and the factors driving Amphibian sites and Roadkills Maxent models were similar: the environment of amphibian populations located nearby roads is similar to the conditions where amphibians are road-killed.
Studies of other fauna groups exposed several factors that drive the occurrence of roadkills: traffic volume, vehicle speed, occurrence of fences [72], distance to vegetative cover and distance to wildlife passages or culverts [8] determine the cause of roadkills. In the particular case of amphibians, road ditches are positively associated with the amphibians' road mortality as they frequently store water during rainy nights [11]. Amphibians specialized in reproducing in temporary waters actively look for road ditches, where they are road-killed [11]. Indeed, the phenology of species activity is the main factor determining roadkill temporal patterns [75]. Furthermore, roadkills tend to occur far from wildlife passages or culverts [8]. Independent variables like distance to forest cover, density of intersections between forest edges and animal abundance are the factors explaining the animals' presence on roads [72]. For example, birds were road-killed on roads surrounded by dense vegetation in Wales [69] and on roads with roadside vegetation with low height in Australia [76]. Mammals were road-killed with higher intensity in areas with protective cover and abundant forage on roadside verges in Australia [76]. Distance to town and distance to water also contributed to explain most of the variation when modelling roadkill occurrence [34].

Modelling of Number Of Roadkills and Hotspots
All Roadkills and Hotspots GLMs showed similarly low values of correlations, although significant probably because the sample size was very high. Number of roadkills depended on traffic and proximity to forest, water and wetland areas, while the number of hotspots depended exclusively on traffic. The negative trend presented by traffic (roadkills increase when traffic decreases) is a result of the skewed distribution of the Roadkills dataset. Although an increment of roadkills with traffic is expected [10,11], our model failed to predict the very high values of roadkills in Slovenia, deeply weighted by a large number of low values of roadkills. Thus, the model predicts more easily the sections without hotspots. The proximities to water and urban and agriculture areas were also important in Maxent models. Both GLM and Maxent analyses for the Roadkills dataset provided similar results. Again, our results are related to habitat quality [10,11,21,25].
Amphibian roadkill hotspots were associated with broadleaved forests [11] with natural vegetation [10], which are main habitats for amphibians. Roadkills in most vertebrate groups are associated with dense tree cover and lower number of buildings [32,34,69,74] as well as with protected areas [15]. However, it is not only habitat conditions that drive roadkills. For example, roadkills of owls are associated with good habitat conditions for species occurrence and also with specific conditions that promote hunting behavior near roads [33].

Conclusions
Slovenia has a huge number of amphibian roadkill hotspots, which likely affects amphibian populations. Unfortunately, we still do not know whether amphibians living nearby roads suffer population decline. Studies analyzing amphibian population trends around roads are needed. Our results confirmed that habitat conditions and presence of water bodies are the main factors driving roadkills. Determination of hotspots is considered crucial in order to introduce mitigation measures on roads [40,77]. The results of this analysis and empirical study [40] are very similar in determining the most critical road sections (5% hotspots), which we believe should be the first priority in conservation efforts on Slovenian state roads. However, sites with good habitats but low current roadkill rates are also important locations to introduce restoring measures [67].