The meadow viper Vipera ursinii includes four subspecies with five allopatric areas of distribution in Europe. It is currently considered one of the most threatened reptile species on the continent, mainly because of its patchy distribution and concurrent habitat loss. Taking advantage of a database composed of occurrence data from bibliographical sources and field observations, we present the first European-scale assessment of the historical knowledge and chronogeonemy of this species. In addition, we evaluate the habitat use and coverage of protected areas with regard to both actual occurrences and modelled potentially suitable areas. This was done for Vipera ursinii s.l. as well as for each of the four subspecies. Our results show different patterns of historical knowledge as well as different degrees of legal protection, depending on the country and subspecies considered. Furthermore, most of the occurrences are from habitats which are classified as vulnerable. A gap analysis reveals an inadequate protection status for modelled areas of potential suitability and a heterogeneous coverage of protected areas, again depending on the subspecies considered. Our findings assist towards a more focused conservation management of all V. ursinii subspecies in the next future, which could take place by connecting landscape-scale research with field studies to update management strategies of protected areas. For these latter, Europe-wide coordinated actions are required to promote plans targeting the same conservation goals.
Since the late eighties, reptiles are receiving attention, because many populations are assumed to decline and many species are considered threatened (Christainsen 1981, Honegger 1981, McNeely 1992, Langton and Burton 1997, Gibbons et al. 2000, Reading et al. 2010, Todd et al. 2010, Böhm et al. 2013). During these years, research on the European herpetofauna increased greatly, especially focusing on distribution and conservation status (Sillero et al. 2014a). Based on the work of many local and national associations, which contributed significantly in filling distributional knowledge gaps (Jacob et al. 2007, Loureiro et al. 2008, Creemers et al. 2009, Corti et al. 2010), nowadays the current spatial distribution of European reptiles (and amphibians) is available online (Sillero et al. 2014a, b). According to this database, there are 145 species of reptiles in Europe, and almost 20% of them are considered threatened according to the Biodiversity Information System for Europe (BISE 2018). The majority of these threatened species are endemic to both Europe and the 27 European Union member states (EU 27), with particularly high endemism richness in the Iberian peninsula, in the Balkans and within Mediterranean islands (Cox and Temple 2009, Silva et al. 2009). Reptiles often show very specific habitat requirements and are particularly vulnerable to anthropogenic threats, so that they require urgent conservation measures (Böhm et al. 2013).
The meadow viper, Vipera ursinii (Bonaparte, 1835), is one of the most threatened reptiles in the whole Europe (Újvári et al. 2002, Filippi and Luiselli 2004, Santos et al. 2006, Gvoždík et al. 2012, Péchy et al. 2015). In fact, it is classified as ‘Vulnerable’ in the IUCN Red List under the B2ab(iii) criterium, which refers to range fragmentation and continuing decline in population size and/or habitat quality, both for the entire Europe and for the EU 27 (Cox and Temple 2009, Joger et al. 2009). The taxonomy of V. ursinii was strongly debated in the last decades; consequently, three taxa previously classified as V. ursinii subspecies were recently elevated to the species level (Nilson et al. 1995, Mizsei et al. 2017, 2018): V. eriwanensis (Reuss, 1938), V. graeca (Nilson and Andrén, 1988) and V. renardi (Christoph, 1861), whose distribution is reported in Fig. 1a. Vipera ursinii shows a patchy distribution, ranging from southeastern France to eastern Romania (Mizsei et al. 2017, 2018; see also Fig. 1 and 5). Within this large area, two main habitats are occupied by the Meadow viper: the warm-dry lowland steppe grasslands (0–400 m a.s.l.) and the montane subalpine-alpine meadows (above 1000 m a.s.l.) (Nilson and Andrén 2001, Crnobrnja-Isailović 2002).
Many Vipera ursinii populations experienced severe declines during the last decades due to several reasons. Habitat loss is the most evident one across the whole distribution range (Újvári et al. 2002, Santos et al. 2006): intensification of agriculture alters viper habitats in lowland grasslands (Péchy et al. 2015, Mizsei et al. 2018), while forestation in montane meadows reduces the number of available basking sites (Filippi and Luiselli 2004). Habitat fragmentation is another critical factor. In fact, road construction in grasslands leads to the segregation and subsequent isolation of small sub-populations (Filippi and Luiselli 2004) which is followed by loss of genetic variation, potentially causing local extinctions (Újvári et al. 2002). High densities of European wild boar Sus scrofa, an increase of winter sport activities with the accompanying infrastructures, intentional killings, as well as climate change (which could boost upward shift of tree line, favouring forestation of montane meadows), also continue to threaten meadow viper populations (Filippi and Luiselli 2004, Edgar and Bird 2005, Joger et al. 2009, Zamfirescu et al. 2011).
Considering all these threats, high attention is addressed to the conservation of the meadow viper in many National and European protected areas (PAs), such as National Parks and Sites of Community Importance (SCI) of the Natura 2000 ecological network. Despite these noticeable conservation efforts, population declines are observed (Újvári et al. 2000, Baillie et al. 2004, Edgar and Bird 2005, Zamfirescu et al. 2012). In the following, we evaluate coverage of PAs of the observed and the potential distribution of this species by means of a large dataset of occurrences. We build accumulation curves of knowledge through historical and current field observations as well as museum records. We describe the corresponding chronogeonemy, i.e. the spatial distribution of the recorded observations over time. Finally, we evaluate the current risk status, in terms of EUNIS habitats, for habitats occupied by V. ursinii and we reveal potentially suitable areas with the help of species distribution modeling (SDM) techniques calibrated using a range of climatic variables.
Material and methods
Target species and study area
The meadow viper Vipera ursinii (Bonaparte, 1835) includes four subspecies, namely V. u. macrops (Méhely, 1911), V. u. moldavica Nilson, Andrén and Joger, 1993, V. u. rakosiensis Méhely, 1893 and V. u. ursinii (Bonaparte, 1835), hereafter named macrops, moldavica, rakosiensis and ursinii, respectively. Meadow viper populations from France have been long recognized in the past as a separate subspecies (V. u. wettsteini Knöpfler and Sochurek, 1955), but recent researches dealing with V. ursinii taxonomy based also on molecular data grouped such populations with the Italian ones as V. ursinii ursinii subspecies, so we perform our analyses according to this latter taxonomic framework (Ferchaud et al. 2012, Mizsei et al. 2017, 2018). We exclude from our analyses V. eriwanensis, V. graeca and V. renardi, because all three former subspecies are currently considered valid species. We analyse the four records of V. ursinii located in Croatia in the Velebit area separately, because the taxonomy of these populations still remains unsolved (Ferchaud et al. 2012, Mizsei et al. 2017). Vipera ursinii exhibits a relatively large but highly fragmented distribution (Fig. 1). Within the study area large environmental datasets are available from the European Union, such as the Natura 2000 and the European Nature Information System (EUNIS). Moreover, the distribution of Vipera ursinii falls within these boundaries, with the exception of macrops, whose range spans also the non-EU countries Bosnia and Herzegovina, Montenegro, Serbia, Albania and Macedonia.
Data description
A dataset comprising 560 occurrences records was generated, ranging over a time frame from 1867 to 2018. It was created by integrating data from bibliographical sources, museum records and own field work (Fig. 1). Bibliographical sources included data from grey literature such as conference papers and technical reports. Data with high geographical uncertainties, e.g. from old museum records reporting the occurrence within a wide area, were excluded from the dataset building process. In particular, the data selected for the analyses were defined either by GPS coordinates or by a fine-scale geographic precision (approximately 1 km2), corresponding to small and well-demarcated areas (e.g. specific toponyms).
A validation of these records was carried out through the Mizsei and colleagues' revised distribution atlas of meadow and steppe vipers in Europe (Mizsei et al. 2018). Within this atlas, grid cells (resolution: 50 × 50 km) are assigned the status of ‘verified' (a category reporting confirmed observations over time), ‘new' (reporting new records from areas where the species was not previously recorded), ‘unverified' (no evidence of presence), ‘historical' (records before 1992) and ‘error’ (erroneous identification of species and/or location). We take advantage of these information to verify whether our records fall within these cells, so as to validate if the occurrence localities we gathered and used for our analyses are comparable with the current knowledge about Vipera ursinii in the whole Europe. Unfortunately there is a risk of threatening a population of a rare and/or vulnerable species even further through the publication of detailed occurrence records (Luiselli 2004, Tulloch et al. 2018). Therefore, we do not publish the detailed dataset which anyway remains available upon request.
Gap analysis
To assess the current protection status of V. ursinii, a gap analysis was performed in ArcMap 10.0 (ESRI 2010), evaluating the overlap between the existing network of protected areas (PAs) and species' occurrence records, considering Nationally Designed Areas, also named Common Database on Designated Areas (CDDA), and European Natura 2000 sites. This spatial information was downloaded from the websites of the European Environmental Agency (< www.eea.europa.eu/data-and-maps/data/nationally-designated-areasnational-cdda-12>) for CDDA sites and of the European Commission (< http://ec.europa.eu/environment/index_en.htm>) for Natura 2000 sites. In the latter, both special protection areas (SPA) and sites of community importance (SCI) were taken into account. SPAs indicate territories recognised as important for the conservation of avian species, designed under the ‘Birds Directive' (79/409/EEC, amended in 2009 and then included in the Natura 2000 programme). SCIs describe areas protected at the European level to ensure preservation, or restoration to favourable conditions, of the habitats and species listed in Annexes I and II of the ‘Habitats Directive' (92/42/EEC). The analyses were carried out with records collected between 1992 and 2018 to ensure the evaluation of the actual (and most recent) status of conservation. Duplicated occurrences, i.e. localities monitored over multiple years, were excluded. Occurrence records outside of PAs were analysed separately, calculating the Euclidean distances from the boundaries of PAs (‘Euclidean distance' raster tool, setting 20 km as maximum distance, in ArcMap 10.0) and then extracting these distance values through the ‘Extract values to points’ tool in ArcMap 10.0.
A further gap analysis was performed to evaluate the degree of protection of the potentially suitable areas of distribution identified for each subspecies through an ensemble modelling process with the aforementioned information on PAs.
Accumulation curves and chronogeonemy
We built accumulation curves to assess the increase of knowledge on meadow vipers' occurrences in Europe, excluding records with no information on the year of observation. Moreover, a chronogeonemy analysis which returns a geographical representation of the increment of knowledge about species in their distribution area was performed in ArcMap 10.0 (ESRI 2010), for each of the four subspecies.
Current habitat protection status and species distribution modelling
The EUNIS habitat classification was downloaded from the geo-portal of the European Environmental Agency (< www.eea.europa.eu/data-and-maps/data/ecosystem-types-ofeurope>) in raster format. In particular, the EUNIS habitat classification level 2 was used in a resolution of 100 m, to assess viper's habitat preference per risk category; based on the forecasts reported in Janssen et al. (2016). The evaluation was performed on the ‘recent' dataset (1992–2018), using the ‘Extract values to points’ tool in ArcMap 10.0 (ESRI 2010) to extract records per risk category.
Model building and calibration
Species distribution models (SDMs) use information about the occurrence/abundance of the target species at surveyed sites and values of a set of variables (usually referred to as ‘predictors') on such sites to estimate the relationships between environmental conditions and species' occurrence/ abundance (Elith and Leathwick 2009). Once these relationships have been estimated, the models can be projected in geographic space to depict the species' potential distribution. Here, SDMs were built with the aforementioned dataset of occurrence records, considering the 19 bioclimatic variables available from the Worldclim.org repository (Hijmans et al. 2005) as potential predictors at a spatial resolution of 30 arc-seconds. To avoid possible multicollinearity of variables, a correlation matrix among all 19 predictors was built in ArcMap 10.0 (ESRI 2010) to extract variables' pairs whose values show Pearson's |r| > 0.85 across the study area extent (Elith et al. 2006, D'Alessandro et al. 2018, Iannella et al. 2018a, Cerasoli et al. 2019). Within each pair, the variable having less ecological importance based on available bibliographical information (Tomović et al. 2004, Filippi et al. 2011, Strugariu et al. 2011, Lyet et al. 2013, Mizsei et al. 2016) was excluded from model building (Supplementary material Appendix 1a). To avoid spatial autocorrelation of occurrence data, the occurrence records were processed through the ‘spThin' package in R (Aiello-Lammens et al. 2015) for each subspecies. Afterwards, a Moran's I test was performed in ArcMap 10.0 (ESRI 2010), to further check for possible residual spatial dependence using the sum of occurrences in each cell (dimension: 30 arc-seconds) for each subpecies’ range as the dependent variable.
Species distribution models can be built through various model classes: envelope algorithms, which take into account only environmental conditions at occurrence localities to derive regularly-shaped niches in environmental space; regression-based techniques, which implement linear, quadratic and/or polynomials functions to fit the predictors–response relationships; machine learning techniques, which do not assume any ex ante data structure and aim at finding dominant predictors–response patterns directly from the data they are trained on (Elith et al. 2008, Qiao et al. 2019). Since the different model classes were shown to provide more or less accurate predictions based on several factors (e.g. complexity of the real-world relationships between environmental variables and species response, quality and size of occurrence datasets) (Elith and Graham 2009), combining SDMs built through various algorithms into so-called ensemble models has been suggested as a promising strategy to investigate the potential distribution of species (Araújo and New 2007).
We built ensemble models (EMs) through the ‘biomod2' package (Thuiller et al. 2016) in R environment (< www.rproject.org>), combining the single SDMs obtained using: generalized linear models (GLM) and multivariate additive regression splines (MARS), as representatives of regression-based techniques; gradient boosting models (GBM, commonly known as boosted regression trees, BRT) and Maxent as representatives of machine learning algorithms. Parametrization of each algorithm is reported in Supplementary material Appendix 1b. Ten sets of 500 pseudo absences (Barbet-Massin et al. 2012) were generated by means of the ‘surface range envelope' (SRE) algorithm, setting the quantile at 0.95, so that pseudo-absences were selected outside the 95% of the multivariate climate envelope built considering only occurrence records (Thuiller et al. 2016). Model calibration was then performed through the ‘BIOMOD_ Modelling’ algorithm.
Model evaluation and post-modelling analyses
For each algorithm and each set of pseudo absences, five evaluation runs were performed, each time using a randomly chosen subset (80%) of the available occurrences to calibrate the model and the remaining 20% for the model validation, yielding a total of 200 models (4 algorithms × 10 sets of pseudo-absences × 5 evaluation runs) for each of the four subspecies. Discrimination performances of the SDMs were assessed through the true skill statistics (TSS) (Allouche et al. 2006) and the area under the curve (AUC) of the receiver operator characteristics curve (Phillips et al. 2006). Only the models exceeding a threshold of TSS > 0.85 and AUC > 0.7 were selected to enter the ensemble modelling process. There is a tradeoff between TSS, which is less prone to overfitting, and AUC, which represents a powerful threshold-independent discrimination metric (Iannella et al. 2018b). The EMs were built using the ‘wmean' (weighted mean of probabilities) and the ‘cv' (coefficient of variation) functions. ‘Wmean' leads to an EM in which the contributions of the single models are weighted based on their respective AUC and TSS, while the ‘cv’ is used to highlight areas with high or low reliability of the EM's predictions based on the agreement among individual SDMs (Thuiller et al. 2016).
All spatial information resulting from the described modelling framework was further analysed in GIS; rasters of potentially suitable habitats were processed with CDDA and Natura 2000 shapefiles using the ‘Intersect' and ‘Extract by mask’ tools in ArcMap 10.0 (ESRI 2010).
Results
The validation of the database of Vipera ursinii s.l. occurrence records (excluding duplicate occurrences; n = 515) resulted in 94.6% of records falling within the cells reported in the revised distribution of Mizsei et al. (2018): 452 points fall within ‘verified' atlas' cells, 28 in the ‘new' cells, 0 in ‘unverified' cells, 6 in ‘historical’ cells and 0 in ‘errors'. The remaining 6.4% of our dataset was outside from Mizsei et al. (2018) range; indeed, these occurrences, taken from verified bibliographical sources, correspond to extinct populations (Kovács et al. 2002, Edgar and Bird 2005). After correcting for spatial autocorrelation, the remaining number of occurrence records were: macrops n = 50, moldavica n = 62, rakosiensis n = 44 and ursinii n = 109. Moran's I test resulted in: I = –0.034 (expected = –0.022), z-score = –1.100 and p = 0.271 for macrops; I = –0.014 (expected = –0.016), z-score = 0.391 and p = 0.695 for moldavica; I = –0.009 (expected = –0.009), z-score = 0.014 and p = 0.988 for rakosiensis; I = –0.020 (expected = –0.023), z-score = 0.274 and p = 0.783 for ursinii, showing no residual spatial autocorrelation among corrected occurrence records (random distribution) for all subspecies.
The gap analysis performed on Natura 2000 sites shows that 383 points of our dataset fall inside special protected areas (SPA), 395 in sites of community importance (SCI) and 375 are covered by both levels of protection (i.e. occurrences located in territories where SPAs and SCIs overlap), whereas six points fall within special areas of conservation (SAC) (Fig. 2a). Further, following the classification of the IUCN's categories for CDDA sites, the gap analysis resulted in 332 points falling within category ‘II', 34 in category ‘IV', eight in category ‘V', two in category ‘NA', one in category ‘Ia' and zero in the ‘Ib’, ‘III' and ‘VI' categories (Fig. 2b). Figure 2c–d, illustrate for each subspecies and each category of protection, the contribution in percentage per country. It is worth to notice that macrops' occurrences falling within CDDA sites span several countries and IUCN categories (Fig. 2d), with an overall low level of protection (Fig. 2b). Notwithstanding rakosiensis is distributed in Hungary and Romania, its occurrences have low protection status (in terms of both Natura 2000 and CDDA sites), with no CDDA sites protecting this subspecies in Romania (Fig. 2d). On the other hand, Romania hosts all the Natura 2000 and CDDA sites which protect moldavica (Fig. 2c–d). The Natura 2000 framework protects, on average, the 50% of all the subspecies' occurrence localities (Fig. 2a). Moreover, ursinii shows a disjunct distribution between Italy and France, with the two countries differently contributing to its protection, depending on the considered PAs categorization: Natura 2000 sites protect the subspecies mainly in Italy; France hosts the total of IUCN's V category, while the remaining categories (i.e. II and IV) are mainly found in Italy. Finally, of the four records of the Velebit area (Croatia) belonging to V. ursinii with unclear taxonomy, one falls within a SPA, one within a SCI, one is covered by both SPA and SCI and one is outside the Natura 2000 network. Concurrently, two of these localities report a ‘doubled' degree of protection, as they are covered by the Natura 2000 network as well as by local protected areas with a ‘Not Reported’ IUCN category status, while the remaining two occurrences are not covered by any CDDA category.
Considering the distribution of the occurrences outside the current network of PAs, a high number of localities can be found within 1–2km from Natura 2000 sites' boundaries, especially for moldavica and rakosiensis (Fig. 3) with a sharp decrease of occurrences with increasing distance from the boundary. Nevertheless, many macrops' records, as well as some occurrence records of moldavica, are found very far from the boundaries of PAs (Fig. 3).
The accumulation of occurrence records over time (Fig. 4) for ursinii starts around 1920 (apart from the species' description in 1835), slowly increases during the twentieth century until the beginning of the new one, where it shows a rapid increase of ‘published’ records, followed by a lower increase that seems to level out until today. The pattern obtained for rakosiensis displays a lack of information from the 1860s until the 1890s. Records begin to accumulate afterwards with a steep increase until the early 1900s, followed by a constant accumulation over time until today. The accumulation over time for moldavica shows a gradual increase during the 1960s–1990s and a sharp increase around 2000, levelling out after 2010. In contrast, occurrence records for macrops present an initial phase of slow increase from the end of the nineteenth century to the end of the twentieth century, followed by a far steeper increase during 1990s which continues up until today.
The chronogeonemy analysis performed on Vipera ursinii s.l. resulted in four classes of observations' dates between 1835 and 2018 (Fig. 5a–d): 1835–1910, 1911–1950, 1951–1991 and 1992–2018. Most of the older data stems from eastern Europe, with another small set of old data from the Balkan peninsula. Records from intermediate time classes are mainly reported from France and, to a lesser extent, from Italy. Recent occurrences are homogeneously distributed between the Balkan peninsula and eastern Europe.
Considering the chronogeonemies of the different subspecies, macrops (Fig. 5a) displays a non-homogeneous increase in published knowledge of occurrence records, with only few historic data and the majority of information provided recently. The chronogeonemy of moldavica presents a complete absence of old information (the oldest occurrence is dated to 1957); on the contrary, this subspecies has been intensively investigated in recent years (Fig. 5b). The chronogeonemy of rakosiensis shows a quite uniform trend (Fig. 5c), with groups of both old and recent data available. A lack of information on intermediate records is also observed for this subspecies. Finally, a lack of occurrence records between old and recent data in Italy is observed for ursinii, with a small number of old and intermediate date records (between 1870 and 1910, and between 1911 and 1950) and many recent data (after the early 1990s) (Fig. 5d, lower). On the contrary, most of the gathered occurrences in France refer to old or intermediate date records, with few recent data available (Fig. 5d, upper).
The analysis of EUNIS habitats harbouring Vipera ursinii s.l. indicated that 278 occurrences correspond to habitats belonging to the category ‘Vulnerable’, which comprises the species' most used habitats such as mesic grassland (E2) and littoral zone of inland surface waterbodies (C3); 129 occurrences fall in the ‘Least Concern' category, primarily including the inland cliffs, rock pavements and outcrops (H3), arctic, alpine and subalpine scrubs (F2), alpine and subalpine grasslands (E4); 41 occurrences are assigned to the ‘Least Concern-Near Threatened' category, represented by the miscellaneous inland habitats with very sparse or no vegetation (H5); 33 occurrences are classified as ‘Endangered' including arable land, market gardens (I1) and seasonally wet and wet grasslands (E3); lastly 34 occurrences correspond to buildings of cities towns and villages (J1), and low density buildings (J2) categories, which are classified as ‘Not Applicable' by IUCN (Fig. 6a). The evaluations of EUNIS habitats for the single subspecies are represented in Fig. 6b–e; a major percentage of occurrences belonging to ‘Vulnerable' habitats is found for all the subspecies, along with a considerable proportion of records belonging to ‘Endangered’ habitats.
The ‘wmean' ensemble models (EMs) (details are reported in Supplementary material Appendix 1b for models' performance, Supplementary material Appendix 1c–j for marginal response curves and Supplementary material Appendix 1k for variables’ contribution) resulted in high values of TSS and AUC for all the four subspecies analysed (mean TSS across subspecies = 0.94 ± 0.02, mean AUC across subspecies = 0.99 ± 0.01), coupled to a clear separation between areas with high and low predicted suitability (Fig. 7).
The gap analysis applied to these potentially suitable habitats resulted in most of the PAs covering scarcely suitable territories (Fig. 7). In particular, macrops and ursinii present a high level of protection, in terms of area covered by Natura 2000 sites, but generally associated to low habitat suitability. moldavica and rakosiensis show a similar trend, but the protection in this case refers to CDDA areas. On the other side, for each subspecies, some PAs also cover territories with high potential habitat suitability, even though they have smaller extent.
Discussion
The protection status of the meadow viper Vipera ursinii was analysed at a continental scale, targeting the overall V. ursinii complex as well as each of the four subspecies inhabiting Europe. Both the Natura 2000 sites and the national protected areas were considered to assess the degree of protection the meadow viper currently receives within the different EU countries. Moreover, current conservation status and future perspectives of EUNIS habitats hosting V. ursinii occurrences were investigated.
The gap analysis performed on the Natura 2000 sites highlighted that a considerable part of the currently known occurrence records is protected, also considering the ‘double degree’ of protection provided by SCI and SPA sites. Similarly, CDDA classification depicts a good level of protection for V. ursinii s.l.; in fact, most of the occurrences fall within National Parks (IUCN category II) followed by the habitat/species management areas (category IV). Nonetheless, all V. ursinii s.l. occurrences, except those located within strict nature reserves (category Ia), are affected by human presence. This means that the authorities responsible for the management of these areas should take into consideration appropriate conservation actions, to avoid habitat loss and intentional killings.
Analysing the legal protection of the four subspecies reveals some interesting trends. It is evident that the results of both Natura 2000 and CDDA gap analyses are affected especially by ursinii (Fig. 2a–b), which seems to currently receive stronger protection than the other subspecies. This pattern can be linked to the high number of available ursinii records as well as to the high number of inhabited habitats falling into PAs. Conversely, the other three subspecies contribute less to the overall degree of protection of V. ursinii s.l., sometimes exhibiting a lack of protection coverage which is dependent on the characteristics of their respective distribution. In particular, results of the gap analysis from macrops occurrences shows a complex protection status. Indeed, this species is partially preserved by Natura 2000 sites (both SCI and SPA) in Croatia as well as Bosnia and Herzegovina, but it does not benefit from considerable protection in these countries by the national designed areas (CDDA), which are absent in Croatia while are comprised in the category ‘Not Applicable’ in Bosnia and Herzegovina. On the other hand, the CDDA sites protect this subspecies in Albania, Serbia and Kosovo while in Montenegro the protection status within CDDA sites results in the category ‘Not Applicable'. This situation is due to the differences between the respective state laws, within and outside of the European Union. A different conservation status emerged for ursinii, whose occurrence localities mostly fall inside both Natura 2000 and CDDA sites, and many of them are covered by multiple levels of protection; in particular, most of these protected occurrences refer to central Apennines populations, which are protected as a whole. The remaining protection of ursinii populations refers to French PAs, which accord special attention to these populations because of the recognised threats such as habitat modification, fragmentation and the corresponding genetic depression (Ferchaud et al. 2011).
Low levels of CDDA protection with regard to the occurrences falling within Natura 2000 sites emerged for both moldavica and rakosiensis. However, the very short distances between the boundaries of PAs' and the occurrence sites of moldavica and rakosiensis suggest that a better protection may be obtained for these subspecies with little efforts by including them in spatial planning of both CDDA and Natura 2000 sites. A partially similar situation is observed for ursinii, for which some of the French localities are located just few kilometres outside one or more PAs.
On the contrary, a completely different protection strategy is needed for macrops. Only few of its occurrence sites fall within PAs or in their respective neighbourhoods; many localities are far away from PAs, at least 20km away from their borders.
Analysing the accumulation of occurrence records over time for each subspecies, it is evident that discovery of new records of rakosiensis has been more gradual than the one observed for other subspecies, resulting in a linear but stepped trend for about 100 years between two main increases of knowledge about new occurrence localities. The steep increase of new records in the 2000s is probably due to the extensive field campaigns which took place in Hungary and Romania under the LIFE+ projects NAT/HU/000116 (European Union 2004) and NAT/RO/006404 (European Union 1999), respectively, to estimate densities and viabilities of local populations.
On the other hand, the 2000s resulted in a pivotal point for the research performed and published on macrops, moldavica and ursinii. The trend in accumulation of records over time for macrops reflects the social and political issues which affected the countries hosting this subspecies during the second half of the 20th century. Indeed, the war events leading to the raise of the Socialist Federal Republic of Yugoslavia after the end of World War II, as well as those of the Yugoslav Wars (1991–2001), may have slowed down research on biodiversity within the countries involved in these conflicts, resulting in few novel macrops occurrences recorded until the late 1990s. Furthermore, the occurrence of macrops in six countries, despite some of them being part of the European Union where a harmonised approach could be expected, leads to different criteria for nature protection, which in turn affect the individual and possibly non-coordinated studies (Vernes et al. 2017). Differently, the noticeable increase in the accumulation of moldavica and ursinii's records since the late 1960s and especially in the 2000s reflect the great interest for these two subspecies in Romania, France and Italy, where PAs represent important centres for studies and projects dedicated to the meadow viper.
The chronogeonemy of Vipera ursinii s.l. is the distribution of occurrence records according to sampling dates, describing the ‘geographical knowledge’ over time. The non-homogeneous distribution of dates of observation is evident for the species complex as well as for each of the four subspecies.
The chronogeonemy of macrops highlights an important number of recent data (from 2007 to 2019) from the Balkan peninsula, which also confirm the oldest records, dating back to 1897. As mentioned before, the lack of information between 1913 and 1990 is probably due to the harsh sociopolitical conditions affecting Balkan countries during most of the 20th century. A considerable number of recent observations, resulting from an evidently increasing research interest, emerge from the chronogeonemy of moldavica; on the other hand, this pattern reveals scarce information of past knowledge on the distribution of this subspecies. No data are reported before 1980s from the Danube delta and no data before 1957 are available from north-eastern Romania.
It is worth to notice the temporal gap within rakosiensis chronogeonemy at the border of Austria and Hungary, with Austria hosting only old observations due to the extinction of the local populations in the period 1890–1910 (Kovács et al. 2002, Péchy et al. 2015). In the remaining part of its distribution, rakosiensis has mostly recent observations. On the other hand, the typonominal subspecies ursinii presents a good level of knowledge distributed evenly over time in the Apennines populations, with a noticeable increase of records in the last thirty years due to specific field studies promoted and published by local PAs.
The meadow vipers resulted to occur in habitats with a high risk of negative change in the future (Janssen et al. 2016). 61% of occurrence records fall in two threatened categories, with 55% classified as ‘Vulnerable' and 6% as ‘Endangered'. This means a serious risk of extinction for populations in habitats E2 (mesic grasslands) and I1 (arable land with unmixed crops grown by low-intensity agricultural methods), due to the forecasted alterations (Janssen et al. 2016). In fact, these two habitats are nowadays subject to rapid modifications and predicted to undergo dramatic changes in the future, due to agricultural intensification and loss of traditional farming which lead to alteration of vegetation structure and chemical pollution of soils because of the widespread use of fertilizers (Janssen et al. 2016). Focusing on the EUNIS habitats, the results are also negative for each subspecies, especially for rakosiensis and moldavica. The former has 68% of occurrences in habitats categorised as threatened, ‘Vulnerable' (E2) and ‘Endangered' (I1), while the latter reports 63% of occurrences in the same two categories. Although ursinii and macrops currently appear to be less vulnerable than rakosiensis and moldavica to habitat modifications, more than 50% of occurrence records fall into critical categories (‘Vulnerable’ and ‘Endangered'). The occurrences of moldavica, rakosiensis and ursinii falling within urban habitats (coded as J1 and J2, where the habitat threat categorization is ‘Not Applicable') are probably linked to occasional findings of individuals moving from (and towards) more suitable habitat patches, or to the presence of natural habitats next to human facilities. Despite the few occurrences falling into urban patches do not imply the common use of these habitats, these findings highlight the possibility of intentional killing as an additional menace, which vipers (as well as other snake species) commonly suffer.
Finally, the results of the gap analysis performed on potential distributions derived from SDMs show that possible migration areas in territories with equally or more suitable climatic conditions for the meadow viper would not be sufficiently covered by the PAs belonging to Natura 2000 and CDDA networks. This is evident for all the subspecies, suggesting the need of a future rearrangement of PAs' boundaries to cautiously manage the legal protection of the meadow vipers. The necessity of carefully planning appropriate conservation strategies for the meadow viper is reinforced by previous researches indicating that in part of its distribution some historic populations went extinct and the remaining ones are mostly geographically isolated and small-sized (Kovács et al. 2002, Luiselli 2004, Krecsák and Zamfirescu 2008, Ferchaud et al. 2011, Lyet et al. 2013). Moreover, the research carried out by Ferchaud and colleagues (2011) on the French ursinii populations demonstrated noticeable inter-population genetic divergence coupled with geographic isolation, making it necessary to treat small groups of close populations as separate management units. Vipera ursinii s.l. is considered one of the most threatened reptiles in Europe, given the high fragmentation of its distribution and the various pressures threatening its preferred the habitats. Indeed, the montane open meadows typically inhabited by ursinii and macrops are predicted to become more and more constrained by forestation linked to the abandonment of traditional shepherding practices and to the upward migration of treeline due to climate warming (Lyet et al. 2013, Ferrarini et al. 2017). On the other side, eastern European lowlands hosting moldavica and rakosiensis are menaced by intensive agriculture and urbanization (Zamfirescu et al. 2011, Péchy et al. 2015). Our analyses highlight that, notwithstanding the effort in conforming PAs categories and management at an international level, the protection status of the meadow viper heavily depends on the country it lives in. As an example, a thorough assessment of the protection status and size of ursinii populations was recently carried out in France within Natura 2000 sites, in the context of the European LIFE 06 NAT/F/000143 programme (Lisse et al. 2012), while fewer monitoring campaigns on other subspecies such as macrops and moldavica have been performed up until now. This supports our statement that additional detailed analyses are desirable also in Balkan countries, due to a current lack of PAs covering meadow viper distributions.
In conclusion, we suggest that future research should focus on both local- and landscape-scale studies on all V. ursinii subspecies. Indeed, we encourage research combining information on size, viability, genetic variability and connectivity of populations, with evidences from landscape-scale research, thus defining how land use and other human activities affect V. ursinii habitats. This would permit to define optimal conservation strategies for the currently-established PAs and for those planned to be established in the next future. Finally, we encourage stakeholders (parks, environmental organizations, local managers) to undertake awareness campaigns throughout all the countries hosting the meadow viper, to sensitize people to the detrimental consequences of habitat degradation upon this and other rare species as well as to prevent intentional killings, which unfortunately still represent a common practice in some portions of V. ursinii distribution.
Acknowledgements
Authors are particularly thankful to all colleagues who helped during field campaigns and to the Editor whose comments improved the clearness of the manuscript. Conflict of interest – The authors declare no conflict of interest.
References
Appendices
Supplementary material (available online as Appendix wlb-00604 at < www.wildlifebiology.org/appendix/wlb-00604>). Appendix 1.