Road mortality of water snakes in light of landscape structure and traffic intensity in north-eastern Hungary

Study area

The studied road sections are located in the Tokaj-Hegyalja region, in the south of the Zemplén mountains and north of Tokaj town (E 21.4124, N 48.1215) in north-eastern Hungary. In the south they are bordered by two floodplain areas, the Taktaköz in the southwest and the Bodrogzug in the southeast that share most of the habitat types of the region but differ in hydrological regime. Both areas are dominated by floodplain meadows, forests and arable lands, whereas the hillside areas are dominated by meadows, forests, and vineyards (Fig. 1). The construction of a system of levees and drainage canals led to the cessation of floods in the Taktaköz floodplain (referred to as “flood-protected area” heretofore), whereas the Bodrogzug area has remained an active floodplain (Boldogh et al., 2016).

Study species

Primarily, we were interested in understanding which species were most readily susceptible to road mortality. To accomplish this, we surveyed every snake carcass encountered within the parameters of the study site. Results revealed that the majority of carcasses belonged to two semi-aquatic snake species, the grass snake Natrix natrix and the dice snake Natrix tessellata. We found only few carcasses of other species occurring in the area which were eventually excluded from this study due to low sample sizes. Both Natrix species are usually confined to aquatic habitats, such as rivers, lakes and nearby habitats. Natrix natrix is a habitat generalist often found far from water in forests and meadows where its amphibian prey occur, while the more habitat specialist piscivorous N. tessellata spends more time in and around water (Speybroeck et al., 2016).

Field methods

We surveyed a total of 58 km of roads which transecting areas where the land use types characteristic to the region are well represented (Fig. 1). All roads surveyed were similar in width with two lanes and two-way traffic.

We estimated snake road mortality by counting the number of carcasses as an index of roadkill incidents, due to sampling errors not considered in the estimated rates. Surveys were conducted biweekly throughout the entire activity period (late March to early November) of snakes in 2012, 2013 and 2020. One of the authors (MS) surveyed all road sections in one direction in a one-day duration riding a bicycle with a steady speed (c. 20 km/h). Upon finding a snake carcass on the road, the surveyor recorded the specimen’s location with a GPS device (Garmin eTrex 10, accuracy ± 3 m), identified the species, measured its size (snout-vent length, SVL) and determined its sex when possible. Snake carcasses were also removed from the road to avoid multiple counting. Using literature sources (Madsen, 1983; Ciesiołkiewicz, Orłowski & Elzanowski, 2006; Borczyk, 2007; Ajtić et al., 2013) and SVL measurements, we applied a proxy age classification based on body size, where snakes were classified into three categories as juveniles (SVL <20 cm), subadults (20 cm <SVL <30 cm) or adults (SVL >30 cm), both for females and males. Although females are usually larger than males and mature at a larger size, this classification allowed us to estimate the age of snakes when sex determination was not possible due to the bad condition of carcasses (n = 87 N. natrix, n = 19 N. tessellata). We determined the sex of adult specimens by visually examining the cloacal region and inferred male sex for individuals that had a thicker tail base. We used the same age classification and sex determination for each studied species.

Analyses of temporal and spatial patterns

We divided the roads into 58 1-km-long sections and summed the number of observed snake carcasses per section per survey as primary datapoints. We used this method instead of kernel density estimation or another method with data transformation because we assumed that exact numbers in road sections might be more informative for local stakeholders when the results are presented for conservation interventions. We used the R 4.0.2. statistical environment (R Core Team, 2021) for statistical analyses and QGIS 3.18 (QGIS Development Team, 2021) for spatial analyses and visualisations.

We assessed temporal and spatial variation in road mortality by comparing the number of snake carcasses found per month and per road section. We built generalized linear mixed-effects models (GLMMs) using the ‘lme4′package of R (Bates et al., 2015), to analyse the effects of possible interactions between months and age/sex categories in the first GLMM and between road sections and age/sex categories in the second GLMM on the number of road-killed N. natrix and N. tessellata, with year as random factor. We also compared the number of snake carcasses per month to a random equal distribution across months using χ² tests. A similar analysis for road sections was not performed due to many zeros, which violated the assumption of tests based on the χ² distribution. χ² tests were also used to compare whether the sex ratio of roadkilled snakes differed from unity.

To make our results comparable with other road mortality studies we calculated observed road mortality rates both temporarily and spatially by dividing the number of carcasses with the number of days in the study period (from the first day of sampling to the last day summed for the three years) and with the number of road sections for both species (Medrano-Vizcaíno et al., 2023a; Medrano-Vizcaíno et al., 2023b).

Analyses with environmental drivers of land use, flood regime and traffic intensity

To evaluate the effect of land use on snake road mortality, we used the Ecosystem Map of Hungary (Ministry of Agriculture of Hungary, 2019; Tanács et al., 2019). This map contains recent (2010-2018) data on land use types in raster format, with a cell size of 20 × 20 m. For analyses, we reduced the number of land use types from 39 to 22 by merging similar types, and types with a low, 1–2% spatial share into broader categories (e.g., Low buildings and High buildings into Buildings). We then calculated the proportion of land use types in 2-km buffer zones around the 1-km road sections. While the home range of these species might be smaller (Wisler, Hofer & Arlettaz, 2008), we chose this distance to include most of the known and suspected hibernating sites in the foothills and the foraging sites such as oxbows and canals in lower-lying areas and to include individuals that may move beyond their home ranges, e.g., during seasonal migration to/from winter hibernating and summer feeding sites. Because floods in road surrounding areas could also influence the number of road-killed snakes, we divided the road sections into two groups based on the flood regime of the neighbouring areas. Snake movement might be affected as floods maintain more wetlands which are important for Natrix species (Speybroeck et al., 2016). Thirty-six road sections were adjacent to regularly flooded areas (Fig. 1, sections 18–39 and 45–58) and 22 sections (sections 1–17 and 40–44) were surrounded on both sides by areas protected from floods by levees. To characterise traffic intensity, we obtained daily traffic volume from the website of Hungarian Public Roads (2020). Traffic data were available in 5 to 10-km sections, and one such section could cover several of our road sections. When road section borders in the traffic data fell within one road section in this study, we averaged the two values for that road section.

The Hungarian Public Roads website and the Ecosystem Map of Hungary only shows recent data, thus we assumed that land use and traffic volume did not vary significantly over the study period. Land use change in the region mostly occurred in the 1990s (Prokopová et al., 2018) and land use was rather stable in the study period (M.S., pers. obs.). To check whether traffic volume was the same between the years, we used traffic data from a 17-km road section (sections 19–35) for which data were available from 2013 and 2021. We found that traffic volume did not differ between 2013 and 2021 in these sections (t-test: t₁₁ = 0.249, p = 0.808).

To model the effect of landscape structure on snake road mortality, we further reduced the 22 land use types into four latent variables (LV1–LV4) in a Bayesian ordination approach using the ‘boral’ package of R (Hui, 2016). Latent variables are used in a similar way as principal components from a PCA but have an advantage over them because they are not linear combinations of the original variables yet they still reduce variance inflation arising from multicollinearity among the original independent variables (Mizsei et al., 2020). For the interpretation of results, we calculated Pearson correlation coefficients between LV1-4 and the original data on the proportion of land use types. These correlations showed that LV1 represented a gradient from rural to urban habitats, with significantly negative correlations with habitats such as sessile oak-hornbeam forests, fruit plantations and arable land and significantly positive with habitats like buildings, paved roads and urban forests. Within the same logic, LV2 represented a gradient from dry to wet habitats from dirt roads, arable lands and railways (negative correlations) to rivers, riverine willow-poplar forests and mixed woodlands (positive correlations). LV3 represented a gradient from artificial to natural habitats with complex cultivation patterns, vineyards and other artificial areas (negative) to tall marsh vegetation standing in water, wet meadows with periodic flooding and lakes (positive), and LV4 was a gradient from closed (woody) to open (treeless) habitats from urban forests, urban green areas and mixed woodlands to arable land, closed mountain grasslands and lakes (Fig. 2).

We assessed the effect of environmental drivers (landscape structure, flood regime and traffic volume) on snake road mortality by constructing a GLMM for each species using the ‘lme4′package of R. We built GLMMs with the number of individuals per section per survey as discrete dependent variable, and flood regime (flooded vs. flood-protected as binary variable), traffic volume, survey week and the four latent habitat variables (LV1-LV4) as fixed explanatory variables and by specifying Poisson error distribution with logit link function. Road section ID and survey year were used as random effects to control for spatial and temporal non-independence of road sections and survey dates, respectively. Then a model selection approach was used in an information-theoretic framework (Burnham & Anderson, 2002) to identify models with substantial empirical support based on Akaike differences (Δi = AIC_i –AIC _min <2.0) and to find the best set of explanatory variables for the number of carcasses found in each species. After this we performed model-averaging of the best models using the ‘MuMIn’ package of R for each species (Bartoń, 2019). When either of the LVs was significant in the averaged model, simple linear regressions were conducted to evaluate which land use type had direct positive or negative effect on the number of road-killed snakes. In these analyses, we only tested variables that showed the three strongest correlations with the number of road-killed snakes either in the positive or the negative direction (Fig. 2), and we adjusted significance levels for multiple comparisons by the False Discovery Rate method (Benjamini & Hochberg, 1995).

Numbers of road-killed snakes

We found a total of 1,715 carcasses of snakes: 1,441 belonged to Natrix natrix, 214 to Natrix tessellata, 45 to Zamenis longissimus, 14 to Coronella austriaca and one was Vipera berus. Because sample sizes were low for the latter three species, we restricted statistical analyses to N. natrix and N. tessellata.

Temporal patterns

By visual examination we identified two peaks in the number of carcasses found (Fig. 3A), which differed significantly from a random equal distribution across the months (N. natrix: χ² = 790.45, df = 8, p < 0.0001; N. tessellata: χ² = 137.93, df = 8, p < 0.0001). For N. natrix, the first peak was in May and June, and the second was in September and October (Figs. 3B and 3D), whereas for N. tessellata, the first peak was between March and May, and the second from July to September (Figs. 3C and 3E).

It was possible to measure the length and to assess the age of 461 roadkilled individuals of N. natrix (247 adults, 66 subadults, 148 juveniles) and 111 individuals of N. tessellata (57 adults, 16 subadults, 38 juveniles). Road mortality did not influence either the age categories (GLMM, N. natrix, F_3,20 = 2.58, p = 0.082; N. tessellata, F_3,17 = 1.153, p = 0.356) or sex (N. natrix, F_2,21 = 2.902, p = 0.077; N. tessellata, F_2,18 = 0.611, p = 0.554) of the carcasses (Figs. 3B–3E). Within adults, sex ratio was 1:1.42 females to males in N. natrix and 1.3:1 females to males in N. tessellata, which did not differ from unity (N. natrix: χ² = 2.469, df = 1, p = 0.116; N. tessellata: χ² = 0.689, df = 1, p = 0.540).

The total number of days in the study period was 605, observed mortality rate per day (roadkilled individuals/days in the study period) was 2.382 for N. natrix and 0.354 for N. tessellata.

Spatial patterns

The number of snake carcasses in the 1-km sections ranged from 0 to 202 in N. natrix and 0 to 53 in N. tessellata. The number of carcasses found suggested large differences among the sections and was significantly different from a random, i.e., equal, distribution of carcasses in road sections in both species (N. natrix: χ² = 956.21, p < 0.0001; N. tessellata: χ² = 69.5, p < 0.0001) (Fig. 4A).

We recorded particularly high numbers of carcasses of both species between road sections 18 and 24, and in sections 27 (highest number for N. tessellata), and 35 (Fig. 4A). Road mortality did not influence either the age categories (Figs. 4B–4C; GLMM, N. natrix, F_2,226 = 1.508, p = 0.224; N. tessellata, F_2,102 = 0.884, p = 0.416) or sex (Figs. 4D–4E; N. natrix, F_2,112 = 0.757, p = 0.386; N. tessellata, F_1,68 = 0.094, p = 0.761) of the individuals indicating that risk of mortality was similar regardless of age category and sex on all sections (Figs. 4B–4E).

Observed road mortality rate per sections (divided by the 58 1-km sections) was 24.845 for N. natrix and 3.69 for N. tessellata. Observed mortality rate per sections per day was 0.041 for N. natrix and 0.006 for N. tessellata.

Environmental drivers of road mortality

The GLMM models suggested that the number of roadkilled N. natrix was significantly positively influenced by the closed to open LV4 and week of the year, and significantly negatively affected by the dry to wet LV2 and the artificial to natural LV3, while flood regime, traffic volume and the rural to urban LV1 had no effect (Table 1). The number of roadkilled N. tessellata was significantly positively influenced by flood regime and week of the year, and significantly negatively affected by traffic volume (Table 1). We did not find significant relationships with any of the latent variables.

Finally, the simple linear regression analysis suggested that the number of N. natrix carcasses found was positively influenced by the area of mixed woodlands, riverine willow-poplar forests, tall marsh vegetation standing in water, and complex cultivation patterns, and negatively affected by the area of railways, dirt roads, urban green areas, arable lands, and vineyards (Fig. 5, Table 2).