Distribution of the grass snake ( Natrix natrix ) and dice snake ( N. tessellata ) in Bulgaria

: We summarise the distribution of the two species of the genus Natrix occurring in Bulgaria, based on records from 147 peer-reviewed publications, grey literature, and data repositories, combined with unpublished data. This is the first extensive mapping for N. natrix; records fall in 560 cells of the 10-km MGRS/UTM grid, of which 102 cells (18.2%) were with published information we could not confirm with new data, 175 (31.2%) were with published and confirmed, and 283 (50.5%) were with new localities. For N. tessellata we increased the number of cells with records by 64% compared to the 2011 mapping, by identifying 445 cells with localities: 162 cells (36.4%) were previously published and unconfirmed, 152 (34.1%) were published and confirmed, and 131 (29.4%) were with new data. Gross climatic conditions for records with exact locations were assigned following the Köppen-Geiger classification; the distribution for both species does not seem to be highly correlated to climate as they were found in 9 of the 12 Köppen-Geiger classes present, only missing from the 3 classes that are limited to high elevations in Bulgaria and account for less than 1% of the area. The vertical distribution of the observations supports our knowledge that the species are most numerous at lower elevations (92.4% of records were <1000 m above sea level for N. natrix and 92.6% were <500 m for N. tessellata ). Higher elevations and some lowlands remain relatively understudied and future sampling will likely reveal new localities for both species

The two study species are among the most common, abundant and widespread ophidian species in Bulgaria (eastern Balkans) . Some notable publications on their distribution include one of the first national mapping efforts by Buresch & Zonkow (1934), and biodiversity assessments of more limited territories, e.g., Eastern and Western Rhodopes (Petrov et al., 2001(Petrov et al., , 2006, Natura 2000 protected areas "Ponor" and "Besaparski Ridove" (Popgeorgiev et al., 2010(Popgeorgiev et al., , 2014a(Popgeorgiev et al., , 2014b and "Oranovski Prolom-Leshko" (Malakova et al., 2018), Vitosha Mountains (Tzankov et al., 2014), Vratchanska Planina Mountains (Naumov et al., 2016), and the Bulgarian part of the Lower Danube River . Multiple other publications provide limited reports. Naumov et al. (2011) published a current state of the knowledge on the distribution of N. tessellata in Bulgaria, but this was prior to an increase in the number of active field researchers and a much improved and centralised data collection through the means of the SmartBirds.org system. No recent and exhaustive national mapping of N. natrix has been made.
However, the perceived abundance and commonness of the study species often lead to field observations about them being unrecorded, even by herpetologists. Thus, the genus Natrix remains understudied and under-reported in Bulgaria, including aspects of the distribution and ecological requirements. Therefore, the main aim of this study was to combine published and unpublished observations to produce a comprehensive and up-to-date database on the distribution of the two Natrix species in Bulgaria. Our goals were to 1) update the distribution of N. natrix and N. tessellata on a 10×10 km UTM grid, 2) evaluate their vertical distribution, and 3) evaluate their occurrence within the Köppen-Geiger climatic classes.
In Bulgaria the study species reach 114 cm for N. tessellata and 163 cm for N. natrix (Naumov et al., 2020). Natrix natrix has a ubiquitous distribution, predominantly at lower elevations, but reaching up to ~2000 m (Buresch & Zonkow, 1934;Naumov & Tomović, 2005); no quantified vertical distribution has been published. Natrix tessellata is also widespread, with ~85% of the known localities in Bulgaria occurring below 500 m, and only two observations recorded above 1100 m -the highest being at 1420 m from Rila Mountains Tzankov et al., 2011). Both species inhabit diverse fresh waters (e.g., streams and river courses, temporary and permanent ponds, spills, natural and artificial lakes, reservoirs, marshes, canals, etc.), as well as brackish ones, such as river mouths at the Black Sea .
For analyses, text data were stored and manipulated with a spreadsheet software (LibreOffice Calc v. 6.4, The Document Foundation, Germany). Spatial data were manipulated, visualised, and analysed using QGIS (v. 3.26, the Open Source Geospatial Foundation, USA). Maps were made using ArcGIS 10.3.1 (ESRI, Redlands, CA, USA).
We combined locality data from the following sources, accessed in March 2023: 1) Personal observations of the authors and colleagues, either requested directly or kindly provided with permission from the users of the field data collection system https://SmartBirds.org �, comprising a mobile application (SmartBirds Pro) and a web-based interface; 2) Over 1000 publications (including dissertations, reports, and other grey literature) with information on the Bulgarian herpetofauna, spanning Historia naturalis bulgarica 45 (2023) from 1892-2022; 3) Data from GBIF.org, that aggregates a number of sources including iNaturalist, specialised Facebook groups, etc.; 4) Collection data from the Zoological Research Museum Alexander Koenig, Germany (ZFMK), and the California Academy of Sciences, USA (CAS). Records from the open data repositories (GBIF, ZFMK, CAS) were considered as "published".
We collated the records in a database, where a record usually represents an observation of one or more snakes of the same species in the same place at a given time. We carried out a quality assessment of the location accuracy, species identification, and other available details for each record manually, and removed the dubious data. We also manually identified and removed duplicate records (e.g., observations in SmartBirds.org also available in a publication; identical records from GBIF.org and CAS; observations with the same coordinates and time entered twice by two observers). Locality data that remained were thus assigned into one of three categories: 'exact locality' (accuracy within 50 m), 'approximate' (accuracy 51-250 m in GBIF, or only a UTM-grid identifiable), or 'unclear' (not used in spatial analyses; the accuracy was less than 250 m, or the locality could not be identified unambiguously, e.g., "Kresna gorge").
Records with exact geographic coordinates were either obtained by using handheld GPS units (including mobile phones) during the observation or were digitised at a later date using georeferenced aerial/satellite images and detailed locality description. Using QGIS, observations with exact coordinates were assigned an elevation based on a Digital Elevation Model (DEM) with a pixel size of 40 m, a 10-km UTM cell, and the Köppen-Geiger climatic class (available from Beck et al., 2018; due to missing data in the original raster from the Black Sea Coast, some observations were manually assigned one of the two classes present locally).
From published data lacking exact geographic coordinates but with a locality description, we considered as a separate record every original description of a location where one or more individuals were found and which was different (e.g., date, count, observer) from other descriptions in the publication. When names of settlements or geographic objects were given as reference points for the locations, we used digital and paper maps to assign these as best as possible to a UTM 10×10 km grid using the MGRS naming of cells (UTM zone 35N, datum WGS 1984). When possible, we used unpublished data from the original observers to clarify the locality. Some presence records were already provided as cells in the 10-, 5-, 2-, or 1-km UTM/MGRS grid; few records had precise coordinates, especially before ~2010, when handheld GPS units became more widely available and used locally. We also assigned, if possible, an elevation ("exact", or within a 100-and 500-m band) using maps.
The raw data varied greatly in quality and available details as they have been collected by multiple people, either in a non-systematic way or systematically, but within a limited geographic area and usually within a short time. Also, areas of special interest to herpetologists and tourists (e.g., the Black Sea Coast, Struma River Valley) and those near the major cities tended to be overrepresented. Although this biases the data, our field experience is that the data is partially reflective of the real situation (e.g., in terms of vertical distribution). We considered each record as a single individual and disregarded counts if provided. We have not removed records based on proximity to other records.
The increase in the amount of data positively correlates to an increased search effort by more experts, more unified and easily accessible data collection, and an increased number of publications, and is unlikely to signify increases in the range or population densities overall.
Altogether, we managed to place observations of N. natrix into 560 10-km UTM cells, based on 2320 records ( Fig. 1; Supplementary material 01  We stipulate that continued and more intensive sampling in areas currently lacking reports of the species (especially at lower elevations) would yield new localities. Several such regions with limited data for both species exist that seem potentially suitable based on expert knowledge and Maxent models of their potentially suitable habitats (Kornilev et al., in press): the Danubian Plain, the central and eastern Thracian Lowland, and the Ludogorie and Dobrudzha Regions. For example, Naumov et al. (2011) hypothesised that the lack of records for N. tessellata from Ludogorie and Dobrudzha Regions (square NJ; data for NJ08, NJ18, NJ80, NJ90; see op. cit. fig. 1 with distribution data and fig. 2 for names of major geographic objects referred to in the text) was most likely due to lack of sampling; here we provide data from two previously unpublished cells (NJ20, NJ44), supporting the need for additional sampling to reveal new localities in this part of the country. This is further reinforced when considering that the herpetofaunal data in SmartBirds.org for the area is lacking for a number of common species expected to be found there. Additionally, here we update the known distribution for both species from the vicinity of the Danube River, updating the recent mapping effort along the river . The geographic scope of the previous publication was limited to 10 km from the Danube, and most of the data were from projects targeted at a few protected areas and focusing on the river itself. We hypothesise that the species are distributed all along the Danube and its tributaries, supported by the potential distribution models (Kornilev et al., in press).
Overall, most observations of Natrix spp. were made close to big cities (e.g., Sofia, Plovdiv, Burgas), close to roads, and around herpetologically popular sites and ones where specific studies were made (e.g., the Struma River Valley, the Eastern Rhodopes, and the Black Sea Coast) and within Natura 2000 sites where specific surveys were made as part of scientific/conservation projects. Improving our knowledge on the distribution of common species can aid in identifying areas that are likely undersampled and can thus help obtaining further data on rare species. To collect a large database with observations, it is recommended to use modern methods for recording, managing, and sharing data, such as GBIF or SmartBirds.org. Since 2016, when SmartBirds.org became operational, over 33000 records have been submitted to it. Still, a way to fill knowledge gaps in distribution and minimise some of the collection bias is to increase the amount of data obtained through non-professional biologists. One such underutilised data source in Bulgaria is citizen science -which can be really helpful in obtaining distribution data, especially from locations outside of protected territories that usually are not sampled professionally. Online social networks, that could be used for citizen science, can provide new and interesting data (e.g., Naumov et al., 2020a).
In Bulgaria, both species were found predominantly at lower elevations, especially at 0-100 m, with observations rapidly decreasing with the increase of elevation (Fig. 3). Natrix tessellata is mainly observed below 500 m a.s.l. (over 92% of the records), closely matching the 85% reported in Naumov et al. (2011). Natrix natrix reaches higher elevations and multiple records are found up to 1500 m; only 16 records exist above that. For both species, the highest elevations continue being previously published records. For N. natrix, it is the observation in Naumov & Tomović (2005), reported there at around 2100 m; however, updated precise coordinates by BN and subsequent elevation estimation puts this record at ~2030 m. For N. tessellata, the highest observation was at 1420 m in the Rila Mountains . Six additional records exist above 1000 m: up to 1100 m on the SW slopes of the Pirin Mountains (Beschkov, 1961), 1034, 1055, 1056, Out of the 12 available Köppen-Geiger climatic classes in Bulgaria, 1627 N. natrix and 1129 N. tessellata records were attributed to territories belonging to all but three classes with very little areas, which also correspond to high elevations (Dsc, ET, Dsa) (Table 2; Fig. 4). Generally, the classes with larger area account for higher proportions of the records. Some clear exceptions exist. BSk represents only 7.4%, while harbouring 16.8% and 33.7% of the records of N. natrix and N. tessellata. These records are clustered along the Struma River valley in the SW and along the northern Black Sea Coast. Both of these areas are also highly popular herpetological spots and generally support high-density populations at many locations, which leads to detection bias.
Although the climatic preferences of both species remain understudied, environmental niche models for Bulgaria revealed limited impacts of temperature on their potential distribution (Kornilev et al., in press); although elevation explained 20-24% for both models and precipitation variables contributed ~15%; these likely relate to the presence of water bodies. Furthermore, the northernmost populations of the dice snake are likely the result of recent colonisation during the Holocene climatic optimum, demonstrating geologically rapid response to climatic changes (Marosi et al., 2012). Their extensive ranges further support that both species are environmental generalists that could thrive in diverse climatic zones.
Knowledge of the distribution of Natrix species is important not only for their long-term survival, but because of their ecological role. The diet of N. natrix For visualisation purposes, percents less than 0.3% were increased to 0.3%. on the Balkans is diverse, with adults feeding predominantly on amphibians and fish, while N. tessellata consumes mostly fish (Šukalo et al., 2014;Speybroeck et al., 2016). Locally, both species can reach potentially high densities -for example, up to 5800 N. tessellata/~18 ha island (Ajtić et al., 2013). Table 2. Percent distribution of records of Natrix natrix (N = 1627 records) and N. tessellata (N = 1129) within the available territory of each Köppen-Geiger (K-G) climatic class in Bulgaria (BG, %).   Table 2 for class' descriptions) in Bulgaria.

K-G class
On the other hand, the snakes are prey for a number of avian and mammalian predators, some of which might be protected. Therefore, coupled with their ectothermic biology, this makes Natrix an important component of food webs and provide an ecological service by aiding transfer of energy from aquatic to terrestrial environments.
In conclusion, the new data on the distribution of the two species, confirm our understanding that in Bulgaria the ranges of N. natrix and N. tessellata are continuous and largely overlap. A major limiting factor seems to be elevation, which generally correlates with local climate conditions. A more detailed comparison between the distribution of the two species could be obtained as a result of specific studies on their habitat preferences and requirements. life history traits in a very dense population of Dice snakes. Zoologischer Anzeiger 252 (3) Evolution and phylogeny of the genus Natrix