Species occurrence points registered in Italy between 2000 and 2021 were obtained from the Global Biodiversity Information Facility (GBIF.org 2021). GBIF occurrence datasets DOIs for each species are available in Suppl. material 3: Appendix 3. Only iNaturalist Research-grade observations were retained. Research-grade observations represent the highest quality occurrence data on the iNaturalist database and therefore, the most reliable. In order to be accepted as Research-grade, an observation must include the date of collection, a picture or sound recording and GPS coordinates (https://www.inaturalist.org). We did not filter data for georeferencing accuracy. All occurrence points not falling within the species IUCN range were removed, as well as outliers and duplicate records. The processed dataset was then implemented in the validation process described in the last paragraph of this section.

Area of habitat (AOH) maps were produced, based on a set of species distribution ranges on a global scale. Global distribution ranges were obtained from the IUCN Red List database (IUCN 2020). We retrieved a total of 25 amphibian and 34 reptile distribution ranges of species occurring in Italy, which represent approximately 60% of all Italian amphibian and reptile species (Sindaco and Razzetti 2021). All species live in terrestrial/freshwater habitats, excluding cave habitats. An additional range was added for Vipera walser, a recently described endemic taxon of viper which only occurs in a narrow area in northern Italy (Ghielmi et al. 2016; see also Doniol-Vacroze et al. 2021 for additional taxonomic discussion) The range of V. walser was obtained from Ficetola et al. (2020). As a taxonomic reference, we relied on two databases: AmphibiaWeb (amphibiaweb.com) for amphibian species and The reptile Database (Uetz et al. 2021, http://www.reptile-database.org) for reptile species. We selected species for which we could retrieve the distribution range from the IUCN Red List database. Additionally, the species’ taxonomy had to be accepted by the referenced database.

Species were linked to their natural habitat using Copernicus Global Land Service Land Cover (CGLS-LC100) 2019 classes as a habitat surrogate. Species-Habitat associations and altitude limits were based on monographs on the biology of target species (Sindaco et al. 2006; Lanza et al. 2007; Corti et al. 2011), while for V. walser, we used data in Ficetola et al. (2020). CGLS-LC100 is a global map of land cover produced by the European Union’s Earth observation programme, freely accessible at the official programme website (land.copernicus.eu). A link to a detailed description of all land cover classes is available in Suppl. material 3: Appendix 3. We selected 13 land-cover classes available on the geographic surface of Italy that could be classified as habitat or lack of habitat for the species. The selected categories are available in Suppl. material 4: Appendix 4. The map is available in a standard latitude/longitude grid (EPSG:4326), ellipsoid WGS 1984. The grid resolution is 1°/1008, approximately 100 m at the Equator. We compiled a binary classification to link species habitat to land-cover classes. A score of either 1 or 0 was assigned to the species, according to whether a specific category represents a viable habitat for the species or not, respectively. Furthermore, we compiled altitude limits (minimum and maximum) for each species. As altitude limits of species distributions can vary according to the local environmental and geological conditions, we implemented altitude-limit categories. Maximum altitude limits were assigned with a minimum of 0 m, progressing in steps of 200 m up to 1800 m. For species which can be found above 1800 m, we set the maximum altitude limit to 2500 m, except for Rana temporaria for which we assigned a limit of 2800 m, based on recently reported altitude limits for the species (Tiberti and von Hardenberg 2012).

In order to generate the AOH maps, we followed the procedure presented in Rondinini et al. (2011) and Lumbierres et al. (2021). First, the CGLS-LC100 map was combined with altitude information to produce a base map. Each raster cell in the base map holds information on the land-cover class and altitude data. The first three digits of cell values represent land cover and the three last digits, the altitude in m/10. Altitude data were obtained from the Shuttle Radar Topography Mission (USGS EROS Archive 2019) map, which has a resolution of approximately 30 m. We resampled the altitude map at the resolution of the CGLS-LC100 using aggregation. The base map can be reclassified in order to retain only areas which meet species’ habitat requirements. Species-habitat association scores and altitude limits were combined to produce reclassification files for each species in R version 4.0.3 (R core team 2021). The R code to produce reclassification files is available in Suppl. material 1: Appendix 1. First, in GRASS GIS version 7.8.5 (GRASS Development Team 2020), we used IUCN species ranges to mask the area of the base map which is found out of the species range. Then, the base map was filtered according to reclassification files. Only land cover areas that met the species’ habitat requirements were retained to produce the area of habitat map. Each pixel in the base map belongs to a single class of the land cover. Therefore, the pixel must necessarily belong to one of the classes which are linked to species’ habitat in order to be retained in the AOH map. We assigned a value of 1 to each raster cell classified as available habitat for the species. For each cell falling within the range of the species, but not matching the species habitat requirements, the assigned value was 0. All cells not falling within the species range were classified as NA. An explanatory scheme of the methodology used to develop the AOH maps is shown in Fig. 1. The complete GRASS code to produce the AOH maps is available in Suppl. material 1: Appendix 1.

Figure 1.

An explanatory scheme of the workflow used to produce the area of habitat maps.

AOH maps were validated for 44 species within the Italian distribution range, as the species-habitat association scores are based on knowledge of species’ habitat requirements in Italy. We used a hypergeometric distribution approach, as suggested by Jiménez and Soberón (2020) and implemented in Dahal et al. (2021). Based on a defined number of occurrence points, the hypergeometric distribution describes the probability of a point falling within the area of habitat of a species. This method allows testing for accuracy of species distribution models through the use of presence-only data. The method is based on point prevalence. After defining a confidence interval, the hypergeometric test calculates the upper and lower limits of the function, based on point prevalence and model prevalence. If the point prevalence exceeds the upper limit of the hypergeometric function, the model can be considered better than random. The model prevalence is the proportion of raster cells within the distribution range which were classified as available habitat for the species (Dahal et al. 2021). The point prevalence represents the proportion of occurrence points out of the total number of occurrences within the range of the species falling inside the available habitat (Dahal et al. 2021). For maps’ validation, we used the dataset of occurrence points obtained from GBIF. We only validated AOH maps of the species for which we could retrieve at least five occurrence points falling within their IUCN range in Italy. This threshold was implemented in previous works to evaluate AOH maps (Ficetola et al. 2015). Before testing for model performance, we applied a buffer of 150 m to the occurrence points. The buffer accounts for miscalculations due to possible inaccuracy of the occurrence records (e.g. limits of the GPS device to record coordinates), as well as for errors due to the resolution of the AOH maps (100 m). Subsequently, model prevalence and point prevalence were calculated to test the model performance. The confidence interval was set to 95%. The full R code used for the analysis with comments on settings and parameters is available in Suppl. material 1: Appendix 1. A summary of the model performance test, including the number of points used for each species, is available in Suppl. material 2: Appendix 2. The complete dataset of occurrence records which were used to validate the AOH maps is included in Suppl. material 3: Appendix 3.

A total of 60 Area of Habitat maps were produced. The mean percentage of area of habitat cover within the range of the species is 48% for amphibians and 44% for reptiles, showing that AOH maps represent an important refinement of the original distribution range of the species. However, for some species, the habitat extent is still overestimated due to the resolution and limitations of land-cover data. For instance, this can be true for semi-aquatic species, as we focused on the terrestrial habitat. The habitats where amphibians breed (e.g. small ponds, vernal pools...) are too small to be detected by remote sensing. Therefore, in this case, the habitat extent should be considered an upper estimate of the available terrestrial habitat; the actual estimate is probably lower because only terrestrial habitats that are spatially connected to suitable breeding sites can be occupied by adults (Becker et al. 2007). A complete list of the species for which AOH maps were produced can be found in Table 1. The 60 available maps include 25 amphibians and 35 reptile species occurring in Italy, corresponding to 60% of all Italian species of reptiles and amphibians (Sindaco and Razzetti 2021). An example of an AOH map produced for the four-lined snake Elaphe quotuorlineata is shown in Fig. 2.

Figure 2.

A map showing the area of habitat (AOH) surface of Elaphe quatuorlineata. Habitat for the species (blue) is identified within the IUCN global range.

A total of 44 AOH maps were validated, while for 16 species, we could not retrieve a sufficient number of occurrence points (≥ 5) to test the AOH model performance. All hypergeometric test results are listed in Suppl. material 2: Appendix 2. Out of the 44 validated AOH maps, the performance of 38 (86%) was significantly better than expected under randomness as the point prevalence exceeded the upper limit of the confidence interval (Jiménez and Soberón 2020). The remaining six maps showed a point prevalence percentage remaining below the upper limit of the confidence interval, with five of them having a point prevalence which exceeded the model prevalence. Four of the six AOH models which did not perform better than expected by chance were assessed with a number of occurrence points equal to 11 or lower.

In Europe, habitat loss has been identified as one of the major threats to many amphibian and reptile species (Cox and Temple 2009; Falaschi et al. 2019). AOH maps can be combined with information on land-use change to evaluate how habitat loss is affecting species distribution (Tracewski et al. 2016), to develop monitoring plans in areas suffering habitat loss and to support conservation actions for threatened species. The produced dataset includes AOH maps for four amphibians and one reptile species which are listed under the threatened taxa in the IUCN Red List (IUCN 2021), as reported in Table 1. Three of the threatened amphibian species are facing population decline due to habitat loss and fragmentation. The dataset also includes AOH maps for five reptile species which are globally assessed as “Near threatened” (NT) in the IUCN Red List (IUCN 2021). Habitat degradation is one of the major threats to four of these species (IUCN 2021). Knowing where the habitat of a species occurs is also crucial for conservation planning, for instance, in the process of designing new protected areas. AOH maps can be incorporated into the identification of Key Biodiversity Areas, as they can be used to infer population size of a species (IUCN 2016). The AOH maps are based on species’ habitat requirements in Italy and were validated within the Italian geographical boundaries, but they also provide an estimate of the distribution of the habitat throughout the whole range of the species. Thus, they can be implemented to apply KBA criteria which rely on this information to detect potential KBAs for amphibians and reptiles in Italy (IUCN 2020). Finally, the AOH maps are a useful tool to investigate general patterns and distributions of herpetological diversity in Italy. They can be used to plan research orientated sampling and fieldwork for target species and to evaluate changes in the habitat available for species if new land-use time series are produced. There are some limitations in the use of AOH maps for conservation, the most obvious one being the resolution of the data. Information on available habitat quantity within the land-cover pixel is missing. Moreover, a binary classification is not optimal, as the habitat quality may vary across the spatial surface. In the future, habitat maps may be improved by including local predictors and landscape predictors, thus providing additional information on habitat availability (Dubos et al. 2022).