Introduction

Roads impact biodiversity because these linear perturbations negatively affect the ecosystems they traverse both during their construction and later use (Forman & Alexander, 1998; Van der Ree, Smith, & Grilo, 2015). Some of the main ecological effects that roads have on wildlife are (a) an increase in habitat loss and fragmentation (Andrews & Gibbons, 2005; Jacobson, Bliss-Ketchum, de Rivera, & Smith, 2016; Shepard, Dreslick, Jellen, & Phillips, 2008), (b) degradation of habitat quality adjacent to the road (Clevenger, Chruszcz, & Gunson, 2002; Roger & Ramp, 2009), (c) changes in the composition and structure of communities due to colonization by invasive or generalist species that are able to increase their distributions and abundances through the use of borders or by disturbances caused by human colonization of adjacent areas (Coffin, 2007; Laurance & Balmford, 2013), and (d) direct effects due to mortality resulting from collisions of vehicles with wildlife (Maschio, Santos-Costa, & Prudente, 2016; Pinowski, 2005).

In the case of reptiles, the consequences of the effects of roads range from increased mortality from road kills, disruption of gene flow between populations, alteration of dispersion patterns, and regional changes in species composition (Andrews, Gibbons, Jochimsen, & Mitchell, 2008; Shine & Fitzgerald, 1996). Snakes are one of the reptile groups most vulnerable to the direct and indirect impacts of roads (Andrews, Gibbons, & Jochimsen, 2006; Fahrig & Grez, 1996). They exhibit high vulnerability due to their morphology, life history strategies, behavior, and daily and seasonal activity patterns that determine their distributions and habitat use (Andrews & Gibbons, 2005; Bonnet, Naulleau, & Shine, 1999; Shine & Lambeck, 1985). Also, as ectotherms, snakes are influenced by environmental conditions that impact, at different scales, their foraging patterns, fecundity, and survivorships, limiting when and where a species occurs in different components of the landscape (Welsh, Hodgeson, & Lind, 2005).

It has been proposed that functional traits of reptiles help reinterpret the structural patterns of the communities they inhabit in fragmented landscapes (Carvajal-Cogollo & Urbina-Cardona, 2015; Carvajal-Cogollo et al., 2019). Functional traits affect reproduction, survival, growth, or dispersion (Blaum, Mosner, Schwagger, & Jeltsch, 2011; Luck et al., 2010) and may help elucidate the functional role of the species in the ecosystem (Díaz & Cabido, 2001; Violle et al., 2007). In this sense, life history traits also may help explain general patterns of vulnerability of snake species to road kill, such as higher vulnerabilities in actively foraging, diurnal species (Maschio et al., 2016; Ramos y Meza-Joya, 2018), or in terrestrial species (López-Herrera, León-Yusti, Guevarra-Molina, & Vargas-Salinas, 2016; Quintero-Ángel, Osorio-Domínguez, Vargas-Salinas, & Saavedra-Rodríguez, 2012). Similarly, body size may be a reliable indicator of the response of a species to the transformation of its ecosystem (Pfeifer et al., 2017). For example, highways with high vehicular traffic tend to reduce the abundances of snake species with large body sizes, because these species are more vagile and more prone to venture onto the road surface (Jochimsen, Peterson, Andrews, Gibbons, & Drawer, 2004). Previous studies have examined the importance of ecological traits in determining snake road kill patterns in some areas of the Neotropics (Ferreira & Silva-Soares, 2012; Maschio et al., 2016; Quintero-Ángel et al., 2012; Sosa & Schalk, 2016). In this study, we consider species traits too as well as the effects of spatial and temporal variation.

There are approximately 319 snake species registered for Colombia (Uetz, Freed, & Hošek, 2018), of which 49 are known to occur in the Orinoco region of the department of Meta (Ramírez-Villalba, Gómez, Velásquez, & Mendoza, 2015). At the national level, 10 snake species are classified as having some level of risk of extinction (Morales-Betancourt, Lasso, Páez, & Bock, 2015) and one of the main threats to these species is mortality associated with roads, as well as habitat degradation and loss (Lynch, Angarita, & Ruiz, 2014). In the Llanos ecosystem of the Colombian Orinoco, a high incidence of snake road kill was observed in comparison to other vertebrate groups (Astwood-R et al., 2018). However, little is known about the main landscape and ecological factors that relate to the road kill of snakes in this region. Our general objective was to determine how certain intrinsic factors of the species and extrinsic factors such as seasonality and varying traffic levels are related to snake road kill patterns in a heterogeneous Orinoco landscape. The specific objectives were to (a) compare snake diversity found on the highway and in adjacent vegetation cover within landscapes with differing levels of heterogeneity, (b) establish the relationship between the diversity of road-killed snake species and extrinsic ecological factors such as rainfall seasonality, traffic levels, and landscape heterogeneity, and (c) determine to what extent snake road kill patterns could be explained by intrinsic ecological traits such as body length, temporal foraging dynamics, foraging strategy, and habitat preferences.

Methods

The study was conducted in the eastern Llanos of the Colombian Orinoco, in the Department of Meta, on a highway (and in the adjacent vegetation cover) that connects the municipalities of Villavicencio (4°06′51.12″ N y 73°16′17.92″ O) and Puerto López (4°05′25.29″ N y 72°57′26.33″ O) (Figure 1). This is a two-lane highway with an average pavement width of 10 m and is 79.7 km long, crossing an elevation gradient from 192 to 400 masl. In the areas adjacent to the highway, level plains predominate, with savannas dominated by Moriche (Mauritia flexuosa) palms, gallery forests, and agricultural systems such as pastures for raising cattle (dominated by Urochloa decumbens) and fields of commercial crops such as corn (Zea mays), rice (Oryza sativa), soybeans (Glycine max), palms (Elaeis guineensis), and citric fruits (Citrus sinensis, Citrus hibridus, Citrus latifolia, Citrus reticulate, and Citrus limon; Fierro-Patiño et al., 2004). The rainfall regime in the area is divided into one dry and one rainy season each year, with annual mean precipitation of 2,414 mm and a mean temperature of 25.7°C (Instituto de Hidrología, Meteorología y Estudios Ambientales [IDEAM], 2005).

Figure 1.

Location of the study site in the Department of Meta (green), Colombia. The highway connecting the municipalities of Villavicencio (yellow) and Puerto López (orange) appears in red.

We selected 13 landscapes along the highway based upon inspection of satellite images (Landsat 8, scale = 500 m, January 2016) and initial visits to the area (Supplemental Material, Appendix 1). The criteria used to select and delimit each landscape were that (a) they were located at least 1 km apart from each other in order to have independence among sites (taking into account reptile movement; Mendenhall, Karp, Meyer, Hadley, & Daily, 2014), (b) there were differences among them in traffic levels, and (c) they were comprised of different proportions of vegetation covers (savanna, forest, human settlements, and pasture-fields). Each landscape window was characterized in terms of traffic levels, the nearest distance to a population center, landscape structure and composition, and rainfall during each field trip. Traffic level was estimated by conducting visual counts of the number of vehicles (including motorcycles) passing through the center of each landscape for 15 m periods 3 times a day on 13 different days during each of four seasons (52 total days); with these data, we calculated the median traffic level for each season for each landscape window. A classified image with a 1:25,000 spatial resolution from the GEF-Paisaje Palmero Biodiverso project was used to evaluate landscape structure and composition. ARCMAP 10.2.2 software (ESRI 2014) was used to classify areas into four types of landscape cover: savanna, forest, human settlements, and pasture-fields. For each landscape, we used the FRAGSTATS 3.4. software (McGarigal & Marks, 1995) to calculate 15 metrics at the landscape and class levels for the area defined by a 1 km radius from the centroid of each highway transect. The distance to the nearest population center was measured, only considering settlements that had an area of more than 700 m², measured from polygons in Google Earth Pro (Landsat 8/April 1, 2018, spectral resolution of three bands: 4, 3, 2, at a resolution of 30 m, zoom 500 m, NASA Landsat Program-USGS).

In each landscape, we surveyed snakes in a 1 km permanent transect along the highway and along 10 permanent transects located randomly within the vegetation cover types adjacent to the highway. In this study, each landscape exhibited a different landscape structure and traffic level, comprising a gradient when considered together (Supplemental Material, Appendix 1).

Six field trips of 13 days each were conducted from March 2017 to January 2018. Historical rainfall data were obtained from the nearest meteorological stations (n = 5) of IDEAM (the Instituto de Hidrología, Meteorología y Estudios Ambientales de Colombia;  https://www.datos.gov.co/Ambiente-y-Desarrollo-Sostenible/Catalogo-Estaciones-IDEAM/n6vw-vkfe). To determine whether a field trip had been conducted in the rainy or dry season, the median historical precipitation for each month was calculated and compared to the monthly total for each field trip. We considered any monthly value above the historical annual median to correspond to a rainy season inventory for this study, and monthly values below this median to correspond to a dry season inventory (Supplemental Material, Appendix 2).

During each field trip, two people sampled each landscape. The order of sampling for the 13 landscapes was determined randomly for each field trip. The transects were sampled using the visual encounter survey method (Crump & Scott, 1994) from ground level up to 2 m into the understory (Urbina-Cardona, Londoño-Murcia, & García-Ávila, 2008). Live specimens encountered in the adjacent vegetation cover transects were captured. Each day, snake surveys were conducted from 06:30 a.m. to 12:00 p.m. and from 04:30 p.m. to 07:30 p.m. The total sampling effort for the entire study was 1,404 person-hr, of which 156 hr corresponded to sampling the highway transects and 1,248 hr corresponded to sampling the adjacent vegetation cover transects. On the highway, the corpses of road-killed snakes were collected, measured for total length and then fixed in 10% formal, and deposited in the Museo de Herpetología, Universidad de Antioquia. In the adjacent vegetation cover, each captured individual was measured for body length and photographed before its release. When it was not possible to identify a snake taxonomically in the field, the individual was sacrificed with an injection of anesthetic, according to the protocol approved by the Committee for Ethical Experimentation with Animals of the Universidad de Antioquia in Act No. 106 of October 13, 2015. The specimen was then fixed in 10% formal and deposited in the same museum collection. All specimens were identified to the lowest taxonomic level possible using herpetological keys (Campbell & Lamar, 2004; Dixon, 1989; Dixon, Wiest, & Cei, 1993; Köhler, 2003; Pérez-Santos & Moreno, 1988; Savage, 2002).

Intrinsic factors for each of the species encountered were defined based on four traits defined by consulting the scientific literature for the species (or closest taxonomic level attained with the keys) (Supplemental Material, Appendix 3): (a) Habitat preference, categorized as Fosorial (F), Semifosorial (SF), Terrestrial (T), Arboreal (Ar), Semiarboreal (SAr), Aquatic (Ac), or Semiaquatic (SAc); (b) Temporal foraging dynamic, which refers to the typical hours for active foraging (diurnal, nocturnal, cathemeral); (c) Foraging strategy, in terms of the foraging technique used by the species (active searching or sit and wait); and (d) Total length (TL). This last variable was obtained from the specimens collected in the field, separating the values obtained from corpses on the highway versus individuals captured in adjacent vegetation cover. When it was not possible to obtain a value for TL from the field data (because corpses were too damaged, were tailless, or when live individuals escaped), we used the mean value for TL for that species reported in the literature (n = 5 individuals).

Statistical Analyses

Landscapes were classified into categories based on rainfall values, vehicular traffic levels, distance to the nearest human center, and landscape metrics. With those variables, a matrix of Euclidian distances of the normalized values was calculated and evaluated to inspect for significant multivariate structure diversity among the landscape windows, using a similarity profile routine (SIMPROF; Clarke, Somerfield, & Gorley, 2008) with 9,999 Monte Carlo simulations (Supplemental Material, Appendix 4). This analysis was conducted using the PRIMER software with the permutational multivariate analysis of variance (PERMANOVA) add on (Clarke, Gorley, Somerfield, & Warwick, 2014) allowing us to identify homogeneous landscape groups. To estimate the number of snake species expected in the landscapes (both on the highway and in the adjacent vegetation cover), we constructed species accumulation curves using the EstimateS 9.0 software (Colwell, 2013). We predicted the number of species for the highway road kill data as a function of the species accumulation curves, using the nonparametric richness estimators Chao 1, Chao 2, and Bootstrap. We used the estimators Chao 1 and Chao 2 because they assume environmental homogeneity in the sample (Magurran, 2004) and therefore are appropriate for estimating snake species richness on the highway, and Bootstrap because of its appropriateness for estimating community richness when there is a large number of rare species (Magurran, 2004). We conducted the same analysis for the expected species richness in the adjacent vegetation cover, and in the total area of study, using the estimators Jack 1 and Jack 2, because they assume heterogeneity in the sample (Magurran, 2004). Based on the richness values predicted by these estimators (Chao 1 and 2, Jack 1 and 2, and Bootstrap), we estimated the percent of species that were represented in the actual inventory in this study (completeness, Soberón & Llorente, 1993) for the highway and the adjacent vegetation cover. To compare abundance patterns, composition, and species uniformity on the highway among the landscapes, we calculated range-abundance curves (Feinsinger, 2001) using the STATISTICA 13.3 software (TIBCO Software Inc., Palo Alto, CA), transforming the relative abundance for each species into a logarithmic scale against the hierarchical order of the species, from greatest to least.

We evaluated whether there were spatial or temporal differences in snake assemblages encountered in each landscape by using a PERMANOVA based on the Bray-Curtis similarity matrix and type III partial sums of squares. First, we evaluated whether there were differences in the structure of the snake faunas killed on the highway versus that registered in the adjacent vegetation cover (fixed factor, n = 2, highway and adjacent cover); for this analysis, we first standardized snake abundances for sampling effort in each transect and landscape window. Upon determining that the snake assemblages on the highway and adjacent vegetation cover were different, we decided to conduct separate analyses thereafter; mainly because of intrinsic natural rarity and low detectability of snakes in adjacent vegetation cover types. The second analysis of variance had a three-factor design: landscape diversity (fixed factor, n = 8, different homogeneous landscape groups based on the SIMPROF routine, see results, Figure 2), landscape (random factor, n = 13) nested within landscape category classification through SIMPROF routine (fixed factor, n = 8, different landscape types), and climatic season (fixed factor, n = 2, rainy and dry seasons). The original matrix used to calculate the Bray-Curtis matrix was based on the absolute abundances for each transect in each landscape window. The PERMANOVAs were run for three response variables: assemblage structure, snake species richness, and total snake abundance. We conducted χ² analyses with the Poisson distribution, followed by an estimation of the dispersion index (variance/median) over all data on snake abundances for the highway, to explore the spatial distribution patterns in the data that were not detected statistically using the PERMANOVA. In addition, we conducted a Pearson correlation analysis of the values for vehicle traffic levels and snake road kill abundances on the highway to inspect for a relationship between these two variables.

Figure 2.

Analysis of principal coordinates of snake assemblage structure found (a) killed on the highway and (b) alive in the adjacent vegetation cover. The colors illustrate the landscape windows (capital letters) grouped by similarity profiles (SIMPROF). The first two components for the highway assemblage (a) explained 38.9% and for the adjacent cover assemblages (b) 49% of the variation in snake assemblage structure. The blue letters (lower case) are the codes for each snake species. The length and direction of the arrow (blue) indicate a greater association of the species with the grouped landscape windows.

Using the matrix of ecological traits for the snake species registered as killed on the highway and as alive in the adjacent vegetation cover, we calculated a Gower distance matrix and constructed a regression and classification tree (LINKTREE) using the PRIMER v7 software and the PERMANOVA add on (Clarke et al., 2014). The tree allowed us to determine whether the values of the intrinsic ecological traits allowed for the detection of a classification pattern of the ecological groups of snakes encountered on the highway and registered in the adjacent vegetation cover. The ecological groups were statistically validated through a similarity profile routine (SIMPROF; Clarke, Somerfield, & Gorley, 2008) with 9,999 Monte Carlo simulations. Once the species were classified into ecological groups, we graphed the number of species per ecological group on each landscape window. Finally, to determine the level of association of the landscape windows based on snake species structure and their ecological groups, we generated a heat map (Somerfield & Clarke, 2013) that represented the association index (Whittaker, 1952) of the species by landscape window and compared the landscape windows based on a Bray Curtis similarity matrix.

Results

Diversity of Snake Species on the Highway and Adjacent Vegetation Cover

All of the snakes encountered on the road were road-killed corpses, and all of the snakes encountered in the adjacent vegetation cover were live captures. The structure of the snake assemblages on the highway and the adjacent vegetation cover was different (df = 1, pseudo-F = 23.4, p-perm = .0044), so analyses were conducted separately for each assemblage (Figure 2). With 1,404 person-h of sampling, we registered a total of 119 individual snakes belonging to 33 species. Of these, 92 individuals belonged to 28 species found on the highway (Table 1, Figure 3) and 27 individuals belonged to 11 species found in the adjacent vegetation cover (Table 2, Figure 3), of which five species were not encountered on the highway. In addition, another five species were recorded for the general area while we were in the field, but not during the actual sampling of the transects, both as road kills on the highway and alive in adjacent cover (Tables 1 and 2).

Table 1.

Abundance (N = 92) and Richness (S = 28) of Snake Species Registered in 78 Samplings of 1 km Transects on the Highway in 13 Landscape Windows Between Villavicencio and Puerto López, Department of Meta, Colombia.

Figure 3.

Some of the snake species encountered as road kills (left) and in adjacent vegetation cover (right) of the Villavicencio—Puerto López highway, Department of Meta, Colombia. Left: (a) Amerotyphlops reticulatus; (b) Pseudoboa neuwiedii; (c) Epicrates cenchria; (d) Atractus sp; Right: (a) Helicops angulatus; (b) Atractus univittatus; (c) Corallus ruschenbergerii; (d) Ninia atrata.

Table 2.

Abundance (N = 27) and Richness (S = 11) of Snake Species Registered in 780 Samplings of 100 m Transects in the Adjacent Vegetation Covers in 13 Landscape Windows Along the Villavicencio—Puerto López highway, Department of Meta, Colombia.

The most abundant taxa found on the highway was comprised of Small Burrowing Colubrid Snakes of the genus Atractus, making up 15.2% of all records, followed by the Brown-banded Water Snake (Helicops angulatus), with 13% of the records. In the adjacent vegetation cover, the most abundant taxa was Helicops angulatus, comprising 29.6% of the records, followed by the Cat-eyed Night Snake (Leptodeira annulata) and the Golden Tree Snake (Chironius carinatus) with 14.8% of records each). The taxa that were only encountered in the adjacent vegetation cover were the Yellow Rat Snake (Spilotes pullatus), Racer Snakes of the genus Mastigodrias, Blind Snakes of the genus Epictia, Chironius carinatus, and “Colubridae morpho7.” This last record was a skeleton of an individual belonging to the Family Colubridae found in a rice field that could not be identified to genus.

Based on the richness estimators Chao1, Chao2, and Bootstrap, in the highway inventory we achieved representativeness of 88.8%, 87.0%, and 86.7% of the expected species, respectively (Figure 4, Supplemental Material, Appendix 5). In the adjacent vegetation cover, the richness estimators Jack1 and Jack2 indicated representativeness of 72.5 and 92.6% of the expected species, respectively (Figure 4, Supplemental Material, Appendix 6). The snake species richness was higher on the highway than in the adjacent vegetation cover, and for both areas, the species accumulation curves almost attained an asymptote, with representativeness of 92.6% of the species in the adjacent vegetation cover and 88.8% of the species for the highway.

Figure 4.

Species accumulation curves in adjacent vegetation cover (a) and on highway (b) between the municipalities of Villavicencio and Puerto López, Department of Meta, Colombia. Number of species observed (real = blue line and green line) and the number of species estimated based on the richness estimators for the adjacent cover (▪ = Jack1, ▲= Jack2) and the highway (♦ = Chao1, ● = Chao2, X = Bootstrap). The curves almost attain an asymptote, achieving maximum representativeness of the sampling of 92.6% of the species for the adjacent cover and 88.8% of the species of the highway.

The curves of the range of abundances of snakes on the highway exhibited slopes that differed for the 13 windows nested within landscape diversity, showing differing patterns in the abundance distributions of dominant and rare species. The taxa that most dominated and were most shared among landscape windows were Atractus sp. and Helicops angulatus, dominating in four windows and sharing distributions in five. One dominant species that was present in only one window was the Rainbow Boa (Epicrates cenchria) (Window 11) (Figure 5).

Figure 5.

Range-abundance curves for snake species on the Villavicencio—Puerto López highway, Department of Meta, Colombia. The relative abundance of each species (n/N) was converted to a logarithmic scale against the range of species ordered from most to least abundant. The number (upper line) indicates the number of the landscape window. The letter (lower line) and the color of the curves indicate the group the landscape window was classified into based on the metrics, rainfall, and local conditions of the highway.

Relationships With Extrinsic Factors

The assemblage structure of snakes found road-killed on the highway and encountered in the adjacent vegetation cover was different (Pseudo-F = 23.4; p(perm) = .004). We found no effect of landscape heterogeneity or precipitation seasonality on species richness, snake abundance, or the general structure of the snake assemblages on the highway or adjacent vegetation cover (Supplemental Material, Appendix 7). However, the PERMANOVA for the richness and abundance values for the highway comparing windows demonstrated variation among the different landscape window groups, as follows: Highway richness: Groups 2 and 11 (t = 2.43, p(Mc) = .044), Groups 4 and 6 (t = 3.8, p(Mc) = .033), Groups 4 and 11 (t = 4.14, p(Mc) = .014), and Groups 10 and 11 (t = 2.43, p(Mc) = .046); Highway abundance: Groups 7 and 11 (t = 3.45, p(Mc) = .037), Groups 4 and 6 (t = 4.67, p(Mc) = .016), Groups 4 and 11 (t = 5.15, p(Mc) = .0058), and Groups 10 and 11 (t = 2.41, p(Mc) = .044) (Supplemental Material, Appendix 8). There was a positive relationship (p-Pearson = .92) between the patterns of snake species richness and abundance on the highway (Supplemental Material, Appendix 9). The χ² test (α = .05, df = 3, p < .05) lead us to reject the null hypothesis of random spatial distribution of the abundance data, with the dispersion index yielding a value greater than 1 (ID = 2.52), indicating an aggregated abundance distribution for snakes on the highway, mainly among the Windows 1, 2, 5, and 9 (Figure 6). The Pearson correlation test indicated a positive relationship (r = .76, p < .05) between the total abundance of snakes that were road-killed on the highway and traffic levels.

Figure 6.

Spatial distribution of snake abundances in 13 landscape windows for the highway data, according to the Poisson χ² and Dispersion index. The colors indicate from below (red) to above (blue) the most to least abundance aggregations of snakes on the highway. The bars indicate the traffic level quantified for each of the 13 landscape windows, which ranged from 500 to 4000 total vehicles.

Relationships With Intrinsic Factors

The regression tree (LINKTREE) showed that there were seven ecological groups, based on the ecological traits of the snakes registered on the highway and adjacent vegetation cover, with the first division in the classification being 95% based on the differences in foraging strategy (active searching vs. sit and wait). The second division in the classification of the species was based on temporal foraging dynamics, total length, and habitat preferences, together accounting for 90% of the differences. The third division was based on temporal foraging dynamics, accounting for 58% of the differences. The fourth division was based on habitat preferences, accounting for 55% of the differences. The fifth division was based on temporal foraging dynamics, accounting for 46% of the differences. The sixth division was based on habitat preferences, accounting for 36% of the differences (Figure 7). With this classification, ecological Groups 6, 4, and 5 contained the largest number of species (with 10, 9, and 6 species, respectively), while Groups 1, 3, and 7 only contained a single species (Figure 7). There were no taxonomic effects apparent in the ecological groups in terms of road kill patterns, with each ecological group containing species from the same taxonomic group that were only registered on the highway or only in the adjacent vegetation cover.

Figure 7.

Regression and classification tree (LINKTREE) of the intrinsic ecological traits of the species killed on the highway (<▴>) and documented alive in the adjacent vegetation cover (<▪>). The classification produced seven ecological groups of species. (a) Sit and wait foraging strategy, terrestrial, semiarboreal, cathemeral, median total length less than 1.06 m. (b) Sit and wait foraging strategy, terrestrial, semiaquatic, cathemeral, median total length less than 1.06 m. (c) Sit and Wait strategy, arboreal, nocturnal, median total length greater than 1.06 m. (d) Active foraging strategy and diurnal. (e) Active foraging strategy and cathemeral. (f) Active foraging strategy, nocturnal, terrestrial. (g) Active foraging strategy, nocturnal, arboreal. R = nonparametric median of the multivariate difference (or degree of separation) among species; B% = the mean of the differences in ranges among groups, using the original ranges of the similarity matrix (modified Gower), scaled to give a value of 100% if the first division is a perfect division (that is, R = 1). Below each significant binary division, R and B% are reported and below those values, we reported the trait attributes responsible for this division.

None of the sections of highway sampled by the 13 landscape windows contained all of the ecological groups of road-killed snakes (Figure 8). The sections that had the highest number of species of road-killed snakes (landscape Windows 9, 11, 1, and 2) were not the ones that had the greatest number of ecological groups. For example, Window 11 had seven species of road-killed snakes, but these species belonged to only three of the seven ecological groups, and Window 9 had nine species of road-killed snakes, but only representing four ecological groups (Figure 8). In contrast, Windows 3, 4, 5, 7, 8, and 13, despite having low species richness values, were comprised of species that each represented a different ecological group (Figure 8).

Figure 8.

Number of road-killed snake species by ecological group in each landscape window. The ecological groups were defined in the caption of Figure 7 according to the routine LINKTREE (see Figure 7). E.G. = ecological group.

As a summary, a heat map shows the 13 landscape windows, that despite containing eight homogeneous landscape groups (validated by the SIMPROF similarity profiles: Supplemental Material, Appendix 4), did not have a coherent structure in terms of snake species road kills; this is also shown by the fact that the landscape windows with the greatest Bray Curtis similarity in snake species (Windows 9 and 10, 1 and 5, and 2 and 11) belonged to different groups in terms of landscape heterogeneity (Figure 9). Also, the grouping of snake species killed on the highway that were associated with a particular landscape window was not consistent with the ecological group of the species (Figure 9).

Figure 9.

Heat map showing the association between species and landscape windows. Each geometric symbol allocated to the species name represents its ecologic group (see Figure 7). Capital letters in parentheses represent the similarity group (Supplementary Material, Appendix 1).

Discussion

In this study, we detected 79.5% of the snake species registered for the Llanos region in the Department of Meta (Ramírez-Villalba et al., 2015). This percentage indicates that an important proportion of these species occur in foothill savanna habitat and that the heterogeneous landscapes such as those adjacent to the highway studied, despite having few remnant gallery forest patches and being dominated by crops and pastures, are occupied by a diverse snake fauna (Jochimsen et al., 2004; Jochimsen, 2005; Welsh et al., 2005). This suggests that the snakes apparently use all elements that make up their landscapes (including the matrix) and do not depend exclusively on remnant forest patches, similar to other vertebrate groups like birds and mammals (Van Dorp & Opdam, 1987; Vieira et al., 2009), permitting them to persist in these areas despite the disturbances.

In total, we encountered the remains of 92 road-killed snakes on the highway, with the taxa Atractus sp. having the most records (15.2%). The dominance of this genus in the road-killed assemblage is consistent with reports by Quintero-Ángel et al. (2012) and López-Herrera et al. (2016) for the Andes region. The similarity of our results with previous studies may be due to the methods employed, as all three studies searched highways by walking transects, permitting the detection of species with small body masses and medium-sized mean body lengths, such as species in the genus Atractus (Köhler, 2003; Roze, 1996). Similarly, the fact that Atractus is the most species-rich genus of snakes in the Family Colubridae, with approximately 65 species reported for Colombia alone, may explain its high detection rate, regardless of the sampling method (Uetz et al., 2018). Helicops angulatus was the second most common taxa in our records from the highway (n = 12) and also was the taxa with the highest number of encounters in adjacent areas as well (n = 8). This could be explained because this species is a generalist in terms of habitat use and diet (Carvalho Teixeira, de Assis Montag, & dos Santos-Costa, 2017). This trait may permit the establishment of this species in disturbed areas such as near urban centers or regions with scarce vegetation cover like flooded pastures, where it still might find some of its preferred prey species such as fish and anurans (Cortes-Ávila & Toledo, 2013).

This is the first study of snake road kills in this region that also reports diversity patterns for areas adjacent to the highway. The species accumulation curves indicated high representativeness of our samples, both for the highway and the adjacent vegetation cover, with the estimated richness being greater on the former. Although road kill surveys are known to underestimate diversity because they miss remains that are consumed by scavengers or are completely destroyed by vehicles (Mumme, Schoech, Woolfenden, & Fitzpatrick, 2000), sampling roads does permit better detection of rare species. For example, in our study, we registered the remains of the Reticulate Worm Snake (Amerotyphlops reticulatus), and in a similar fashion Quintero-Ángel et al. (2012) recorded other species difficult to detect in natural habitats on a highway located in the Central Andes of Colombia, such as Blind Snakes of the genus Leptotyphlops. This pattern also could be related to the low detectability of snakes in natural habitats in general, and especially during the day, and to their irregular distributions, reserved and cryptic natures, and often nocturnal tendencies (Lind, Welsh, & Tallmon, 2005), traits that together make them difficult to encounter in areas adjacent to the highway where the vegetation cover is heterogeneous. Given our results, we recommend the use of a sampling design that combines searching highways for snake remains with inventories of live individuals in areas adjacent to the highway, to achieve the most complete documentation of species richness of a region, including the detection of rare species on the highway or species difficult to detect alive in other habitats.

Most species were registered both on the highway and in the adjacent vegetation cover (N = 6, out of a total of 11 species registered in the adjacent cover). However, one of the species not shared by both location types was Chironius carinatus, documented by one individual encountered dead on the highway, but not during the transect surveys. This contrasts with the results of Astwood-R et al. (2018) that were obtained from the same highway in 2015, where this species was the most common road kill (N = 17). This difference between studies may be explained because Astwood-R et al. (2018) sampled more intensively during the rainy season, when species of this genus are more active (Angarita-Sierra, 2014), as this time corresponds to the mating season (Marques, Almeida-Santos, Rodrigues, & Camargo, 2009) and coincides with a time of greater food abundances (frogs; Marques, Eterovic, & Endo, 2001).

We observed a pattern in the structure of the snake assemblage with the range-abundance curves for the highway data, with differences in dominance and in the detection of rare species among the landscape windows, even among windows belonging to the same type of landscape group (SIMPROF). This was also apparent in the PERMANOVA analyses of snake abundance and species richness values among landscape windows. This was because some landscape windows had greater values of abundance and richness than others and had different dominant species. For example, landscape Window 11 was characterized as having predominantly an adjacent vegetation cover of pasture-fields, almost exclusively rice fields, and only 9.7% forest. In contrast, other landscape windows, such as Window 7, belonged to the same SIMPROF group, yet this was the window that had the lowest abundance and richness values on the highway, but with an adjacent vegetation cover dominated by more perennial crops (palm and citric fruits) and a greater proportion of forest (12.9%). This suggests that in more disturbed habitats (pasture-fields), the specific crop present may influence the number and diversity of snake road kills. Also, landscapes with a low diversity of adjacent vegetation cover dominated by crops, but with some significant percentage of forest cover, may constitute dispersal corridors that are associated with movements of some species related to feeding, reproduction, shedding, thermoregulation, and occasionally foraging. For example, Epicrates cenchria was the dominant species in Window 11, probably because it was drawn to this site due to the abundance of mice associated with rice crops (M. T. R. A., personal observation). Freitas, Sousa, and Bueno (2013) also found higher reptile mortality on roads where the adjacent vegetation was herbaceous, perhaps because roads in such areas are more used for thermoregulation. They also suggested that for roads that cross level landscapes, adjacent forest cover could represent preferred dispersal habitat for this group.

Our results indicated an aggregated pattern of road kills, with Windows 1, 2, 5, and 9 having the greatest numbers of road kills. We believe that this aggregation is related to factors like traffic levels and distance to the nearest population center, which were positively correlated with the total abundance of road-killed snakes (Windows 1 and 2 were near Villavicencio, Window 5 near Pompeya, and Window 9 near Pachaquiaro, all important population centers within the study area). This is consistent with previous studies of road kills (Fahrig, Pedlar, Pope, Taylor, & Wegner, 1995; Rudolph, Burgdorf, Conne, & Schaefer, 1999; Shepard et al., 2008) that also found a higher incidence of road kills where traffic levels were higher. However, De Souza, Da Cunha, and Markwith (2014) found a negative relationship between vehicular traffic levels and road kill rates of vertebrates in the Pantanal in Brazil, where elevation played a greater role in explaining these rates. It seems intuitive that road kill rates should be related to traffic levels, but apparently other factors such as elevation, season, habitat characteristics, and even sampling methods may obscure this effect (Monge-Nájera, 2018).

In the case of snakes, the effect of vehicular traffic has been related to their relatively slow rates of locomotion, their tendencies to thermo-regulate on road surfaces, and their reactions to passing vehicles (Andrews & Gibbons, 2005), as well as to cases of intentional road kills (Andrews & Gibbons, 2005; Ashley, Kosloski, & Petrie, 2007; Secco, Ratton, Castro, da Lucas, & Bager, 2014). Our results of the first PERMANOVA showed that the highway and adjacent vegetation cover exhibited different structure patterns of the snake faunas. This could be because several snake species have ecological characteristics that increase their vulnerability on the roads (e.g., tendency to thermo-regulate on the road, peaks of activity corresponding to peaks in traffic levels, habitat requirements that vary according to the season and the reproductive mode; Jochimsen et al., 2004; McCardle & Fontenot, 2015; Rudolph et al., 1999; Seigel, 1986), likely causing biases in road samples towards species exhibiting these types of traits. Additionally, the difficulty for us of sample at night (due to poor security conditions in the area) added to the low detectability of snakes (Lind et al., 2005), perhaps impeding our detection of more species in the adjacent vegetation cover. This result differs from that of Vargas-Salinas, Delgado-Ospina, and López-Aranda (2011), who found similar species compositions in the reptile communities within the forest and on a road in a reserve in the Cauca Valley of Colombia. They attributed this result to the large home ranges snakes require for carrying out activities such as locating mates, food, or shelter. However, the landscape structures sampled in their study were primary forest and did not vary substantially, as did the heterogeneous landscapes in the areas adjacent to the highway sampled in our study.

We did not find a relationship between landscape diversity and snake diversity on the highway based on the second PERMANOVA, despite the variation among landscape windows. This may be explained because the structure of snake communities varies at local scales depending on the distribution of resources such as shelter, mates, and prey, rather than being directly related to landscape diversity. For example, the species most commonly found in our study (the Brown-banded Water Snake Helicops angulatus, the Black-headed Snake Tantilla melanocephala, the Parrot Snake Leptophis ahaetulla, Small Burrowing Colubrid Snakes of the genus Atractus, Racer Snakes of the genus Mastigodrias, and Cat-eyed Night Snake Leptodeira annulata) are habitat generalists with high vagilities (Cortes-Ávila & Toledo, 2013; Silva, 2004) capable of satisfying their needs across a variety of landscapes (Jochimsen et al., 2004), a finding similar to Freitas et al. (2013) and Astwood-R et al. (2018). In contrast, we presume that the movement patterns of other taxa such as Epictia sp. are limited, restricting them only to the highest quality habitats.

We expected to find an effect of rainfall seasonality on the diversity of road-killed snakes, as suggested by other studies with Neotropical snakes. For example, Sosa and Schalk (2016) found road mortality to be higher during the rainy season. However, our analysis failed to show such a relationship. This result may be related to the anomalous number of days with lower than average rainfall that occurred during the 2017 rainy season (between 3 and 9 days with no rainfall even in the wettest months; IDEAM, 2017), which would reduce any seasonal effects that might otherwise exist.

Based on the regression tree analyses, we did not find patterns associated with ecological groups that would permit us to separate the species found on the highway from those found in the adjacent vegetation cover (Figure 7) or detect apparent patterns in terms of landscape window groups. These results are in contrast to reports from other studies of road-killed snakes in the Neotropics, where a greater incidence of road kills was reported for species that were diurnal, semi-arboreal, terrestrial, and with large body sizes (Hartmann, Hartmann, & Martins, 2011; López-Herrera et al., 2016; Maschio et al., 2016; Ramos & Meza-Joya, 2018; Sosa & Schalk, 2016). In Santa Cruz, Bolivia, Sosa and Schalk (2016) reported that terrestrial and semi-arboreal snake species experienced more road kills (71.6 and 8.2%, respectively), with fossorial and semiaquatic species being less affected (3.7%, 2.5%, and 0.4%, respectively). However, Sosa and Schalk (2016) suggested that the road may act as a barrier to fossorial species, which contrasts with our study where the greatest abundance of road-killed snakes in the area of study belonged to the fossorial taxa Atractus sp.

Finally, in our study species with an active foraging mode were detected more often both as road kills and in the adjacent vegetation cover. This is similar to results found by Hartmann et al. (2011) and Sosa and Schalk (2016) for other Neotropical snake faunas, where species with active foraging characteristics also were the most vulnerable to road mortality. Similarly, for temperate zones, most studies show the species with active foraging strategy are more prone to road kills (Bonnet et al., 1999; Forman et al., 2003). This seems to be related to the fact that species with this type of foraging strategy explore a wider range of habitats, are more mobile and, therefore, have a greater risk of mortality when venturing onto roads (Hartmann et al., 2011). However, our findings also could be because the majority of the species known to occur in the study area share an active foraging strategy. Also, species with the sit and wait strategy can be principally nocturnal (as with vipers and some boids), and they may have morphologies or cryptic coloration that makes detection more difficult (Henderson, Pauers, & Colston, 2013; Vitt & Caldwell, 2013; Wasko & Sasa, 2009), so we may have undersampled this component of the snake fauna in our study, where transects were searched only during the day.

Implications for Conservation

Mortality of snakes on the Villavicencio—Puerto López highway is common. In this study, only 78 days of sampling 13 different 1 km sections of a highway documented the remains of 92 individual snakes, belonging to 28 species that died due to collisions with vehicles. Also, this result was not adjusted based on estimates of the detection probabilities of the species (depending on the methods employed), estimates of the rate of removal of cadavers from the road (by scavengers or due to complete destruction by subsequent vehicular impacts), or estimates of the proportion of mortality due to intentional road kills (Ashley et al., 2007; Teixeira, Coelho, Esperandio, & Kindel, 2013). Although we did not find an influence of landscape heterogeneity on snake mortality (apparently because of the nondiscriminating way most species use the adjacent vegetative cover), in areas with intense crop cultivation, such as rice fields, that were associated with forest patches, the abundance/dominance of some species may increase. Thus, a study of the effects of landscape heterogeneity on a larger scale may yield a clearer picture of the importance of certain types of vegetation cover and landscape connectivity on patterns of snake mortality on highways.

Our study also documented that road kills were more common in sections of the highway with higher levels of traffic, regardless of the time of year or the morphological or ecological traits of the snake species. Identifying these points of the highway where mitigation measures may be implemented is crucial, since the development and implementation of strategies and measures to reduce the effects of roads on fauna are costly, both in terms of time and money, impeding implementation of adequate management measures. Thus, our main recommendations for attempting to reduce the threat of highways to snakes in this area (and surely also for other vertebrate groups) are, in order of priority, to install speed radars in areas with the highest traffic levels, to post advisories to warn drivers of the risks of colliding with wildlife species on the highway (as also suggested by Machado Fontes, Moura, Mendes, & Romao, 2015), and to conduct environmental education programs to raise awareness of local drivers about the impacts of road kills on local biodiversity. This is particularly relevant for snakes, as intentional road kills for this group has been demonstrated in previous studies (Ashley et al., 2007; Crawford & Andrews, 2015; Secco et al., 2014). We consider that education programs should place special emphasis not only on current drivers but also on children, emphasizing the role and value of snakes for the ecosystems, following Crawford and Andrews (2015).

All the recommendations mentioned earlier must be accompanied by road management, mainly in areas where nonstructural measures such as control radars and fences may be insufficient, such as, for example, in the places identified in this study as landscapes with low diversity of adjacent vegetation cover dominated by crops, but with a significant percentage of forest cover that may constitute dispersal corridors that are associated with movements of some species, for example, Epicrates cenchria. Proposing to take structural measures is relevant for these areas due to their special characteristics, such as, for example, the lack of illumination at night. These measures could be implemented by the land owners in these areas, thus involving the community to raise awareness about the importance of this biological group.

Supplemental Material