New data on the distribution of the Vipera ammodytes (Linnaeus, 1758) mitochondrial lineages place their contact zone in western Bulgaria

Molecular studies have outlined several mitochondrial lineages of Vipera ammodytes , but the limits of their distribution ranges remain unclear due to limited sampling throughout the vast range of the species. One particularly understudied area is the Eastern Balkan Peninsula where at least three lineages occur, showing ranges that may be partly overlapping. We used two mitochondrial markers (cytb and ND2) to unveil mtDNA distribution patterns based on vipers from 31 localities across Bulgaria. Our results confirmed the presence of the north-eastern (NE) and the south-eastern (SE) mitochondrial clade in Bulgaria, the latter being represented by its southern (S) and eastern (E) subclades. Among the sampled localities, two were suspected to be potential contact zones be-tween these mtDNA lineages based on old morphology-derived distributional records. The NE clade was absent from both potential contact zones. However, our data showed that in western Bulgaria, populations of this clade establish contact with populations of the E subclade in at least one area, and also come close to contact with populations of the S subclade. These results indicate the need for more detailed research in the potential areas of contact in western Bulgaria, integrating morphological data with extensive mitochondrial and nuclear DNA-sampling to better understand the phylogeographic patterns of intraspecific differentiation in V. ammodytes .


Introduction
Despite its relatively wide range, the Nose-horned Viper, Vipera ammodytes (Linnaeus, 1758), is one of the poorly studied European viper species.The species is distributed from north-eastern Italy and southern Austria, through most of the Balkan Peninsula and many Aegean islands, to north-western and northern Asia Minor and the Lesser Caucasus (Speybroeck et al. 2016).Seven different subspecies have been recognized in the past, based on various morphological traits (Boulenger 1913;Bruno 1968;Sochurek 1976;Biella 1983;Biella and Blättler 1989).Using multivariate morphological analyses, Tomović (2006) suggested only four valid subspecies: V. a. ammodytes, inhabiting northern Italy, southern Austria, Slovenia, Croatia, Bosnia and Herzegovina, Montenegro, northern Albania, most of Serbia, north-western North Macedonia, north-western Bulgaria, and western Romania; V. a. montandoni, inhabiting south-eastern Romania, most of Bulgaria, western Turkey, north-western, northern and north-eastern Greece, most of North Macedonia, the southernmost parts of Serbia, and south and central Albania; V. a. meridionalis, inhabiting central Greece, Peloponnese and the Cyclades; V. a. transcaucasiana, inhabiting the eastern parts of Turkey, Georgia and Armenia.
Recent studies based on mitochondrial DNA markers, however, showed that the current taxonomy needs revision because the species exhibits a high genetic diversity throughout its range, and especially in the Balkan Peninsula, which the morphology-based approach fails to reflect (Ursenbacher et al. 2008;Čubrić et al. 2019;Freitas et al. 2020).Seven major mtDNA clades were found within V. ammodytes: NW, NE, MO, PE, CY, SW, and SE, the latter being comprised of the E, S, T, and AM subclades (Ursenbacher et al. 2008).Very recently, a study based on a combination of genomic and mtDNA phylogenies confirmed the presence of all clades described by Ursenbacher et al. (2008) and reported the presence of an additional clade from north-eastern Peloponnese (NEP) (Thanou et al. 2023).Genomic-based analyses revealed two distinct clusters, corresponding to the north-western ("North Balkan Clade" or NBC) and south-eastern ("South Balkan Clade" or SBC) Balkan Peninsula, respectively.The first cluster includes the NW, NE and MO mtDNA clades (all of which belong to the V. a. ammodytes morphotype), while the second cluster consists of the rest of the clades (belonging to V. a. montandoni, V. a. meridionalis and V. a. transcaucasiana morphotypes) (Thanou et al. 2023).The results of Thanou et al. (2023) also indicated very limited or no gene flow between the two clusters.The authors proposed that the NBC and SBC clusters should be considered different species, namely V. ammodytes and V. meridionalis.
In Bulgaria, the Nose-horned Viper is widespread throughout the country, except in the high parts of the mountains or urbanized and intensively cultivated agricultural land (Stojanov et al. 2011).Both V. a. ammodytes and V. a. montandoni morphotypes have been reported from Bulgaria.They are considered to have a parapatric distribution, possibly coming into contact in the western parts of the country.Specifically, a narrow area of contact has been suspected to exist along the Iskar River, northeast from the northern end of the Iskar Gorge (Fig. 1) (Stojanov et al. 2011).Another putative contact zone might exist along the southern slopes of Stara Planina Mtn., the eastern slopes of Ihtimanska Sredna Gora Mtn. and the western slopes of Sashtinska Sredna Gora Mtn., within the "triangle" locked between the towns of Chelopech, Ihtiman and Panagyurishte (Fig. 1) (A.Stojanov personal communication to A. Dyugmedzhiev).A third contact zone might also exist in the valley along the upper reaches of the Struma River, situated on the south-western slopes of Vitosha Mtn.(the area between the villages of Bosnek and Chuipetlovo, Fig. 1) (Tzankov et al. 2014).However, molecular data are not available for these populations and the presence, location and extent of these contact zones have not been confirmed.
On the other hand, the north-eastern (NE) clade as well as the eastern (E) subclade of the SE clade (hereafter referred to as the SE/E subclade) inhabit the country (Ursenbacher et al. 2008;Thanou et al. 2023).The NE clade is confirmed to be present in north-western Bulgaria (Ursenbacher et al. 2008) and its distribution in the country might correspond to the distribution range of the V. a. ammodytes morphotype in Bulgaria, depicted by Stojanov et al. (2011).The SE/E subclade is confirmed to be present along the latitudinal gradient in the central and eastern parts of Bulgaria (Ursenbacher et al. 2008;Thanou et al. 2023).Its distribution in the country seems to correspond to a large degree to the distribution range of the V. a. montandoni morphotype in Bulgaria, depicted by Stojanov et al. (2011), with the exception of south-western Bulgaria, where the southern (S) subclade of the SE clade (hereafter referred to as the SE/S subclade) might be distributed (Thanou et al. 2023).The latter is found in immediate proximity to the south-western border of Bulgaria in northern Greece and western parts of North Macedonia (Ursenbacher et al. 2008;Thanou et al. 2023).However, its presence has not been confirmed yet with molecular studies, since neither Ursenbacher et al. (2008) nor Thanou et al. (2023) sampled populations from south-western Bulgaria.
The goals of the current study are to 1) evaluate the mitochondrial diversity of the Bulgarian V. ammodytes populations, 2) assess the geographical distribution of the mtDNA lineages across Bulgaria, and 3) reveal and define the spatial distribution of putative contact zones between these lineages.In order to achieve these goals, we conducted an extensive sampling, covering most of the territory of Bulgaria, and analyzed information derived from two mitochondrial markers.The widely used gene cytochrome b (cytb) was chosen to allow the comparison of Bulgarian populations with previously published data, while a second gene, the NADH dehydrogenase subunit 2 (ND2), was added to provide further resolution between Bulgarian populations.

DNA sampling
Sampling was carried out between 2014 and 2023 in Bulgaria and included 62 individuals from 31 localities covering different parts of the distribution range of V. ammodytes in the country.Among these localities, we included two of the three potential contact zones of the different V. ammodytes morphotypes in western Bulgaria: along the Iskar River near Karlukovo (43.17°N, 24.06°E) and along the upper reaches of Struma River between the villages of Bosnek and Chuipetlovo (42.49°N, 23.21°E) (Fig. 1).In each locality, ventral scale samples were taken in the field from living individuals (n = 52), which were released after the procedure.Additionally, tissue samples from specimens found dead were collected (n = 10).The material was stored in absolute ethanol for long-term preservation.

Molecular procedures
Total DNA was extracted from tissue samples using DNeasy Blood & Tissue Kit (Qiagen), following the instructions of the manufacturer.Two mitochondrial markers were amplified through polymerase chain reaction (PCR) with the Thermo Scientific DreamTaq Green PCR Master Mix (Thermo Fisher Scientific) and according to the instructions of the manufacturer.Specifically, cytb and ND2 were respectively amplified with the primer pairs L14724Vb / H15914Vb (Ursenbacher et al. 2006 and references therein) and L4437b / tRNA-trp (Kumazawa et al. 1996;Ashton and de Quieroz 2001) and the cycling protocols of Ursenbacher et al. (2006) and Garrigues et al. (2005).Sanger sequencing was performed by Macrogen Europe BV as an external service.
The phylogenetic analysis was based on the cytochrome b fragment, which is commonly used in molecular studies of vipers.DNA sequences of V. ammodytes from a previous European study (Ursenbacher et al. 2008) were downloaded from Genbank and added to the newly obtained data (see Suppl.material 1) to provide context.Sequences were aligned using Mega X (Kumar et al. 2018).The best substitution models for each coding position were estimated with PartitionFinder ver.2.1.1 (Lanfear et al. 2017).Phylogeny was reconstructed through Bayesian inference (BI) accomplished in MrBayes v. 3.2.7 (Ronquist et al. 2012).Four simulations of Markov chains and 4 × 10 6 generations were run, sampling one of every 100 trees.The chain parameters were examined in Tracer ver.1.7.1 (Rambaut et al. 2018).The first 25% of trees were discarded as burn-in.

Genetic diversity of Bulgarian clades
In order to further study Bulgarian populations, we concatenated cytb and ND2 fragments (ND2+cytb hereafter), thus obtaining a larger mitochondrial dataset.Haplotype networks were based on the concatenated mitochondrial alignment using the TCS method (Clement et al. 2002) implemented in PopART 1.7 (Leigh and Bryant 2015).

Genetic diversity and distribution patterns of V. ammodytes in Bulgaria
New molecular data from Bulgaria included 62 DNA sequences from two genes; from one sample only ND2 was obtained.No stop-codons were detected in the protein coding genes.Localities, geographical coordinates and Genbank accession numbers are given in Suppl.material 1.
The cytb alignment, used to infer the general phylogenetic tree of V. ammodytes, consisted of 103 in-group sequences (61 own and 41 previously published data from the species, and one outgroup Vipera berus (Linnaeus, 1758)).The final matrix was 927 bp long.The ND2 fragment used for haplotype networks included new data only and was 864 bp long.
The arrangement of the Bulgarian samples in the Bayesian phylogenetic tree confirmed the presence of the north-eastern clade (NE) and the south-eastern clade (SE) in the country, the latter being represented by its southern (S) and eastern (E) subclades (Fig. 2).The samples collected in the central and northern parts of western Bulgaria (n = 18) belonged to the NE clade.The samples from the southernmost parts of central-western and south-western Bulgaria (n = 14) were arranged within the SE/S subclade.Finally, those collected along the latitudinal gradient from the easternmost parts of western Bulgaria, throughout central to eastern Bulgaria (n = 30) belonged to the SE/E subclade (Fig. 1 and Suppl.material 1).Two individuals from a single locality (Brusen, 42.89°N, 24.12°E) were arranged within separate mtD-NA clusters (NE and SE/E, see below).
The (ND2+cytb) dataset used for haplotype network reconstruction was 1788 bp long.The haplotype network confirmed the presence of three distinct clusters in Bulgaria (Fig. 3).Based on our concatenated dataset, the estimated genetic distance was ca. 5% between NE and E, and ca.2.8% between S and E clusters.

Analysis of individuals from contact zones
In order to investigate the potential contact zones, we sampled 10 individuals from the area along the Iskar River and seven from the area along the upper reaches of Struma River.Additional samples were also collected from two individuals captured near the village of Reselets, along the Iskar River (43.24°N, 24.02°E).All individuals captured from the first area and Reselets belonged to the SE/E subclade, and those from the second area belonged to the SE/S subclade.No individuals belonging to the NE clade were discovered in either of the two sites.
The closest sampled populations from the NE clade and the SE/S subclade occur in the Krayshte region, western Bulgaria, approximately 15 km from each other (Fig. 1).More interestingly, representatives of both the NE clade and the SE/E subclade were discovered within a narrow perimeter of less than 100 m, in the area between the valleys of Malki Iskar and Cherni Vit rivers, near Brusen, north-western Bulgaria (42.89°N, 24.12°E) (Fig. 1).Four neonate V. ammodytes (three females and one male) were captured there on 18.09.2022and DNA was collected from one of the females.The site was visited again on 28.09.2023 when the remains of a dead subadult female were found 7 to 60 m away from where the four neonate vipers were found the previous year.Our analysis grouped the first female within the SE/E subclade and the second one within the NE clade (Fig. 3).

Genetic diversity and distribution patterns of V. ammodytes in Bulgaria
The Bayesian phylogenetic tree constructed in our study was in general agreement with earlier mitochondrial phylogenies (Ursenbacher et al. 2008;Thanou et al. 2023).Our results confirm the presence of both the NE clade and the SE/E subclade in Bulgaria, as shown by Ursenbacher et al. (2008) and Thanou et al. (2023) and reveal that the distribution range of the former reaches further into the east than previously known.Furthermore, we confirm that the SE/S subclade is also distributed in Bulgaria, as depicted by Thanou et al. (2023), extending its distribution range further into the north-east.The NE clade inhabits the central and northern parts of western Bulgaria, the SE/S subclade is found in the southernmost parts of central-western Bulgaria and south-western Bulgaria, while the SE/E subclade occurs throughout the rest of the country (Fig. 1).During the course of this study, we were not able to collect any samples from nose-horned vipers along the Mesta River Valley in south-western Bulgaria (the separate purple triangle area with question mark on Fig. 1).Genetically admixed individuals between the SE/S and SE/E subclades are found nearby, to the south of this area in northern Greece, where their respective distributional ranges might overlap as a result of secondary contact following periods of population expansion (Thanou et al. 2023).So, further DNA sampling is needed in order to clarify the clade-assignment of individuals from this part of Bulgaria.

Distribution within the potential contact zones
According to Thanou et al. (2023), gene flow and hybridization between the NE clade (comprising the NBC genomic cluster) and both the SE/E and SE/S subclades (comprising the SBC genomic cluster) is absent or very limited.The authors considered the two clusters allopatric since the closest observed distances between them were at 50 and 70 km in Serbia and North Macedonia, respectively.Based on the geographical isolation and the deep divergence time, Thanou et al. (2023) suggested that speciation between the two clusters is complete and proposed that NBC and SBC could be considered different species, Vipera ammodytes and V. meridionalis, respectively.The authors, however, did not include samples from Bulgaria (with the exception of only one locality for the SE/E subclade from the south-eastern part of the country); the easternmost samples from the NE clade in their study being from western North Macedonia and central Serbia.Our results show that the NE clade and the SE/S subclade come very close to each other in the Krayshte region (Fig. 1), a pattern that is also suggested by the distributional ranges depicted in Ursenbacher et al. (2008) and Čubrić et al. (2019).During the course of the current study, the SE/S subclade was newly discovered in the southern parts of Krayshte, near the village of Sirishtnik (42.57°N, 22.78°E), and the NE clade was newly discovered with several sampled populations in the northern and central parts of Krayshte, the southernmost of them being located 15 km to the north-west from the SE/S subclade's population, near the village of Stanyovtsi (42.69°N, 22.70°E).Even more interestingly, there is evidence that the NE clade and the SE/E subclade, although parapatric, have reached contact at least in one area between the valley of the rivers Malki Iskar and Cherni Vit in western Bulgaria.This contact zone is situated between several closely sampled localities for both of the two mtDNA lineages within the area between the "triangle" locked between the Malki Iskar and Iskar Rivers (Fig. 1).For the NE clade, these localities are near the village of Bozhenitsa (43.01°N, 23.83°E), along the town of Etropole (42.84°N,24.00°E and 42.80°N,23.95°E),between the villages of Lopyan and Yamna (42.84°N, 24.08°E), and near the village Divchovoto (42.83°N, 24.19°E).For the SE/E subclade, they are near the villages Reselets, Karlukovo, Babintsi (42.93°N, 24.26°E and 42.94°N, 24.27°E), and Goliam Izvor (42.95°N, 24.10°E) (Fig. 1).Several small and medium-sized river valleys pass through this entire area, providing natural and easily accessible bio-corridors for dispersal of individuals (Fig. 1, small picture).Therefore, it could be expected that at least some mixed populations could be found there.
Although two of the localities sampled in this study (labelled as 1 and 2 in Fig. 1), were considered areas of contact between the V. a. ammodytes and V. a. montandoni morphotypes (Stojanov et al. 2011;Tzankov et al. 2014), we didn't find any genetic evidence in either of them that would support the occurrence of such a contact zone.Unfortunately, all potential and existing contact zones remain very poorly studied.Future studies should integrate morphological data with extensive DNA sampling in order to fully understand the ecological drivers of coexistence and segregation and estimate possible gene flow between evolutionary units of V. ammodytes.Molecular analyses should involve nuclear markers in order to better assign individuals to evolutionary units and to shed additional light on gene flow and hybridization.Such studies should also focus on the area locked between the southern slopes of Stara Planina Mtn., the eastern slopes of Ihtimanska Sredna Gora Mtn. and the western slopes of Sashtinska Sredna Gora Mtn.(labelled as 5 in Fig 1), which has also been considered a potential contact zone between the morphotypes.

Figure 1 .
Figure 1.Distribution of the different mtDNA lineages of Vipera ammodytes in Bulgaria.Red, green and blue circles indicate the location of populations from which DNA samples were collected and their assignment to the respective clade; the yellow circle indicates the location of the population near Brusen.The red, green and blue areas indicate the potential distribution of the different clades in Bulgaria based on the distribution ranges depicted by Stojanov et al. (2011) and the authors' personal data; purple area with a question mark indicates the area around the valley of Mesta River, in which the clade-assignment of individuals needs to be clarified by DNA sampling.Numbers 1 and 2 indicate the suspected contact zones based on literature data (along the Iskar and Struma rivers, respectively) which were sampled during this study.The potential contact zones identified in the study are indicated with contours and numbered as follows: 3 -Iskar and Malki Iskar rivers (contact NE -SE/E), 4 -Krayste region (contact NE -SE/S), 5 -Stara Planina, Ihtimanska Sredna Gora and Sashtinska Sredna Gora Mtns (contact NE -SE/E, based on personal morphology-based data, which need confirmation by DNA sampling).

Figure 2 .
Figure 2. Bayesian phylogenetic tree of Vipera ammodytes inferred from the mitochondrial cytochrome b gene.Colored clades/ subclades are found in Bulgaria, each color corresponding to the red, green and blue sampling localities shown in Fig. 1.

Figure 3 .
Figure 3. TCS haplotype network based on a 1788 bp.mitochondrial fragment (ND2+cytb) including all samples from the NE clade, the SE/E and the SE/S subclades, collected from Bulgaria.The size of the circle is proportional to the number of samples showing the respective haplotype.The yellow circles indicate the placement of the individuals from Brusen within the NE clade and the SE/E subclade, respectively.