In this research, carried out in Tuscany, Central Italy, we considered four species of snakes that are commonly widespread throughout the Italian territory. The Atlas of amphibians and reptiles of Tuscany (Vanni and Nistri 2006) reports eight species of snakes in the whole region and four are expected in our sampling site according to the Atlas and excluding any human introduction: three colubrids [Hierophis viridiflavus, the European Whip snake, Natrix natrix, the grass snake (now N. helvetica as in Kindler et al. 2017) and Zamenis longissimus, the Esculapian snake] and one viperid (Vipera aspis, the asp viper). Only the European whip snake (H. viridiflavus) and the Aesculapian snake (Z. longissimus) are included in the Habitats Directive, but we have also applied the monitoring protocol to the other two species.

The studied area is an ecotone characterised by the transition from the Tombolo mixed forest (Pisa, central Italy) and the Tyrrhenian coast (43°38'52.63"N, 10°18'18.28"E) (DMS), next to the Arnino Laboratory (Department of Biology, University of Pisa), formerly a cultivated field and now abandoned. The ecotone is north-south orientated (facing westwards), with a small pathway dividing the area in two contiguous segments of 500 m length each, for a total of 1000 m, 2 m width and an overall estimated surface area of 0.2 ha. It falls inside the range of the SIC (“Sito di Importanza Comunitaria”), now ZSC (“Zona Speciale di Conservazione”) “Selva Pisana” (code IT5170002), an area of the Regional Park “Migliarino San Rossore Massaciuccoli”.

Sampling and data recording were carried out from 1995 to 2017 with different sampling efforts. Continuous and comparable research efforts, during the whole of the snake’s active season, were carried out in 1997, 1999, 2002, 2004–2005, 2016–2017: sampling has been carried out from the beginning of the snake’s outdoor activity to the end of active season (March-October), on average twice a week by two or three operators (one being the same during the whole study period, MALZ) for 5 hrs per sampling session. We performed the sampling by repeated surveys along the ecotonal transect, alternating North to South and vice versa, patrolling the direction and ensuring that the sampled area was covered during all day hours. Snakes were visually detected, captured by hand and carefully examined, sexed, measured and marked for further recognition, with ventral clipping and a serial enumeration system, suitable for at least four years (Brown and Parker 1976). They were released at the capture point after data recording.

The long term analysis considered all different snake age classes, we recorded (newborns, juveniles, adult males, adult females) and we aimed at verifying general life-history dynamics of each species, that is population stability, decrease or increase.

The short term analysis was carried out with the aim of controlling the area of previous research, after ten years of no research, in order to verify any change in occurrence patterns and species abundance. In parallel, the monitoring protocol at national level (Mezzasalma et al. 2016) provides that transects should be four per area, repeated at least six times per site within a six year span. Such a protocol provides confirmation of the presence of snake species, when observations gave positive results during consecutive surveys. Therefore, we have tested the monitoring protocol procedure and suitability on our target ecotone, even if it represents only one site. Therefore, we have also estimated abundance indices upon repeated counts along the selected transect, during the period of maximum activity (from mid April to the first days of June). We repeated the transect six times/year (Mezzasalma et al. 2016), for a total of 12 monitoring days in 2016–2017. Given the very low number of recaptures during the 2016–2017 period (only one H. viridiflavus recaptured, only one N. helvetica with two recaptures and two N. helvetica with one recapture each), it was not possible to analyse the recapture matrix with specific software as Capture, Mark and others. To estimate snake density, we have, therefore, considered the number of captured and recaptured snakes (divided per species) relative to the overall considered surface area, as in Rugiero et al. (2002).

In our analyses, we used the large dataset of seven years of investigation distributed over a 21 years period. The Aesculapian snake (Zamenis longissimus) has been excluded from the analysis due to the very small sample size (only a few records in the last ten years). We applied a χ² on all the data that meet the test assumption, i.e. applicable when not more than 20% of all the cells have an expected frequency lower than five (Siegel 1980). Additionally, in short-term research, we performed a χ² to test possible fluctuations on the total number of captures per each species. We have then analysed average size variations between sampling years per each species using one-way two factors analysis of variance. Data were analysed using the software SPSS IBM (2011). All tests were considered significant if P < 0.05.

In the long-term analysis, we obtained a different pattern depending on the species considered (Fig. 1). Hierophis viridiflavus (8.7 ± 4.7 captures per year) was more common in the first years while it constantly decreased in the last period; on the contrary, N. helvetica (3.7 ± 4.2 captures per year) increased gradually in abundance; V. aspis (12.1 ± 7.4 captures per year) appeared relatively constant throughout the years. The total number of captures per species and year is reported in Table 1. Despite the exclusion of the years 1997 and 1999 (due to the absence of Grass snakes, that prevented the application of a χ² test), we obtained a significantly different distribution amongst years and species (χ² = 24.5; df = 8; P = 0.002). From total snakes observed, per species, we obtained an average occurrence of 6.1 H. viridiflavus/100 m, 2.6 N. helvetica/100 m and 8.5 V. aspis/100 m.

Figure 1. 

Distribution (as numbers) of captures per year of sampling from 1997 to 2017.

Considering the pattern for each year, we recorded 0.3 to 1.7 H. viridiflavus/100 m, 0.1 to 1.0 N. helvetica/100 m, from 0.3 to 2.3 V. aspis/100 m. If we considered the available suitable surface of our transect, approximately of 0.2 ha, we could infer a yearly density variability of 15–85 H. viridiflavus/ha, 5–50 N. helvetica/ha and 15–115 V. aspis/ha. If we consider a theoretic hectare of suitable habitat (i.e. five times our study area), we have to consider the ecotonal area bordering an hypothetical figure, as for simplicity of calculation, on a pentagon of 1000 m side or on a quadrate of 1250 m side: in this hypothetical scenario, our target habitat would represent the ecotones of a much larger area of ca. 172 or 156 ha, respectively. The estimated minimum-maximum densities should therefore be 0.09–0.49 or 0.1–0.54 H. viridiflavus/ha, 0.03–0.29 or 0.03–0.32 N. helvetica/ha and 0.09–0.67 or 0.1–0.74 V. aspis/ha, respectively.

In 2016–2017, during 29 surveys, we captured, marked and recaptured 164 snakes, 60 of them individually marked. Only one H. viridiflavus, three N. helvetica and seven V. aspis were recaptured at least one time, only three V. aspis were recaptured at least three times, unfortunately preventing any further population analyses (e.g. recapture rate, survival etc.). The χ² test confirmed that no significant variation in population density occurred (χ² = 1.611; df = 2; P = 0.447) (Table 2).

In addition, population average body size varied significantly amongst years for each species (H. viridiflavus: F=2.244, df=6, P=0.037; N. helvetica: F=5.892, df=4, P=0.003; V. aspis: F=7.299, df=6, P < 0.0001) (Table 3). Hierophis viridiflavus differed especially in 1997, 1999, 2004 and 2017 with larger size in 1999 and smaller in 2017. Natrix helvetica differed amongst years, with larger size in 2017. Vipera aspis differed amongst all years, with larger size in 1999 and 2016 and smaller in 2002. Notably, body size variations recorded in 2002–2005 are referred to several newborns and juveniles with respect to the other years.

For each year, we considered six randomly selected surveys to test the validity of monitoring techniques. During these 12 monitoring surveys, we captured 31 snakes in total, six H. viridiflavus, six N. helvetica and 19 V. aspis, that is on average 2.58 snakes/day (0.5 H. viridiflavus, 0.5 N. helvetica and 1.58 V. aspis per day).

Our study highlights that monitoring and analysing dynamics and ecology of snake species-community needs a long-lasting sampling effort to assess patterns of variation in species dynamics, in changes of population structure or even in ecology and phenology phenomena (Fitch 1999; Jäggi and Baur 1999; Gentilli et al. 2006; Rugiero et al. 2013). We observed even slight, but significant, variations in population size. In general, H. viridiflavus, N. helvetica and V. aspis represent common species of ecotonal areas of Mediterranean habitats, widely known to be generalist snakes (Filippi and Luiselli 2000). In our research, the average snake occurrence per year was 24.6 snake/yr (172 snakes in seven years), as in other European studies (Luiselli et al. 2011, Natrix tessellata, 21.5 snake/yr; Moreno-Rueda and Pleguezuelos 2007, Malpolon monspessulanus, 16.3 snakes/yr; Capula et al. 2014, Hierophis viridiflavus, 22 snakes/yr), but less than in others (Rugiero et al. 2002, Hierophis viridiflavus, 66 snakes/yr, Zamenis longissimus, 41 snakes/yr; Leliévre et al. 2013, Hierophis viridiflavus, 59 snakes/yr, Zamenis longissimus, 65 snakes/yr). Our data represent, however, low estimates when focusing on a single species, ranging from about four N. helvetica to 12 V. aspis, likely reflecting the carrying capacity of our small area. Our estimation from the hypothetical area with bordering ecotones, shows maximum densities of 0.49–0.54 H. viridiflavus, 0.29–0.32 N. helvetica and 0.67–0.74 V. aspis per ha. Rugiero et al. (2002) stated: “If we pool data for four habitat types, the mean density per year of C. viridiflavus adults was 2.32 ± 1.05 specimens × ha^-1, and the mean density of E. longissima adults was 1.34 ± 0.89 specimens × ha^-1.” Analysing their Table 1 or also the number of marked H. viridiflavus (respectively 730 sighted or 511 marked snakes in 11 years, pooled from four different habitats of 20 ha each; see Rugiero et al. 2002, materials and methods and results), there emerges an almost certain error, with a much lower density estimation, that is from 0.83 to 0.58 H. viridiflavus/ha, respectively and on average ≈ 66 (66.4) snakes captured each year (error confirmed by the author, L. Luiselli, pers. com.). In the paper of Leliévre et al. (2013), of 1608 captured snakes in 13 years, 764 were H. viridiflavus (calculated from Table 1, p. 587). Those snakes were captured and recaptured in four different areas, whose overall surface (calculated according to the reference scale [as in Fig. 1, p. 587], as area 1 = 11.27 ha; area 2 = 7.69 ha; area 3 = 6.93 ha; area 4 = 69.43 ha) is of 95.32 ha. From this overall surface, the estimated density for H. viridiflavus in north-western France (Chizé, Deux Sevres) in 13 years was ≈ 0.6 (0.62) snake/ha with, on average, ≈ 59 (58.8) snakes captured each year. Derived data from Rugiero et al. (2002) and Leliévre et al. (2013) provide results that are, however, much more comparable and in line with our findings, as is logical to be expected in the same species and in a very similar habitat structure, despite the fact that they were recorded in different climatic and geographical areas (Mediterranean climate in central Italy vs. Continental climate, in North-Western France). However, we should underline that, in our research, when comparing the observation of snakes throughout different years (Table 1), the patterns amongst species did differ significantly. Hierophis viridiflavus becomes less frequent as regards the number of individual snakes captured/recaptured. Natrix helvetica, on the contrary, was found for the first time in 2002 and seems to be increasing in number in the following years. Vipera aspis shows a relatively stable occurrence throughout the whole considered period, despite a relatively larger number of newborns captured and recaptured during 2017. In a long-lasting research (Rugiero et al. 2013), it was shown that climatic changes, namely global warming, induced an increase in delay in entering into hibernation and a decrease in aboveground activity of asp vipers (annual onset, sensu Rugiero et al. 2013), but no effect of such phenomenon on snake density has been recorded. We found some similarities in the occurrence patterns only in the impressive and longest snake community research published by Fitch (1999), despite studied species and their life-history traits not being directly comparable. Fitch (1999), in fact, reported strong density variations of Coluber constrictor in relation to the different fields and crops management, with increasing or decreasing herbaceous structure and small mammals’ densities and, accordingly, increased or decreased snake densities.

During 2016–2017, in our study area, the average body size of H. viridiflavus and N. helvetica did not change differently from the recorded significant variations in population abundance, underlining significant effects only on the snake occurrence, as found in Fitch (1999). More importantly, V. aspis, albeit stable throughout the long period, has significantly decreased its density and its body size in the last two years (we mostly recorded sub-adults and juveniles; data not shown). As many snake species are quite long-lived organisms (Madsen 1987; Flatt et al. 1997; Fornasiero et al. 2016) and some of our studied individuals have been found after 9–10 years since first capture, it may be logically hypothesised that the causes that have induced such population variability may have occurred several years before our analysis (Gentilli et al. 2006; Rugiero et al. 2013), most probably involving both H. viridflavus and V. aspis, whose adults’ occurrence in the last two years was actually lower than reasonably expected. Only Jäggi and Baur (1999) were able to support robust hypotheses on V. aspis extinction at a local scale (as forest overgrowing increased humidity and decreased sun exposition) and even our data do not provide information about any possible ongoing extinction process. Prey availability and/or predators presence was not recorded for such a long period (see Saviozzi and Zuffi 1997), thus, as far as we are concerned, it is not easy to formulate any valid hypotheses for what we recorded. Alternatively, we can hypothesise that a fluctuation in snake composition, like the one we found, may represent a natural pattern, whose slight variations are, on average, difficult to detect over a few years, but are quite common and easy to detect in long term research (Fitch 1999; Jäggi and Baur 1999; Sullivan 2000; Rugiero et al. 2002; Moreno-Rueda and Pleguezuelos 2007; Pringle et al. 2009; Leliévre et al. 2013). Actually, natural habitats of and around our study area (i.e. wooded area surface, ecotonal structure) have not shown any evident or measurable change during the period, as emerged from historical vs. recent aerial photographs (Google Earth images), supporting our hypothesis of no habitat effect on snakes dynamics, different from what was observed in northern Italy (Gentilli et al. 2006) and in Switzerland (Jäggi and Baur 1999). The appearance of N. helvetica may represent a new colonisation of the area, given several humid and marshed areas and canals very close to our study area, the transformation of which should be studied in future research.

Finally, we are therefore confident that the monitoring system protocol of the Italian Environment Ministry (Stoch and Genovesi 2016) is a practical, highly informative and easy to apply tool to be carried out by professional herpetologists, even if it would deserve quite a large effort in the field before assessing true population dynamics.