Contribution to the phylogeography of the nose-horned viper (Vipera ammodytes, Linnaeus, 1758) in the Central Balkan Peninsula

Seven genetic clades have been recognized within the species Vipera ammodytes (nose-horned viper); however, the precise phylogenetic position of many Balkan populations is unknown. We used Bayesian analysis of the mtDNA sequences from the 16S rRNA mtDNA gene obtained from 47 individuals (26 novel samples sequenced in this study and 21 sequences available from GenBank). Our results show that sampled nose-horned vipers from localities in Serbia are clustered within three clades: the northeastern (23 individuals), the northwestern (two individuals) and the southeastern (one individual). Results revealed an overlapping distribution of the northeastern and the northwestern clades in two populations. We have revealed that the northeastern clade extends further south than previously suggested, to the Ohrid/Prespa lakes in North Macedonia. Our findings contribute to the knowledge of the genetic diversity of this species in Serbia and help to clarify the geographical distributions of mtDNA-defined clades.


INTRODUCTION
Phylogeographic studies serve as a valuable tool in conservation biology [1][2][3][4], not only to help understand the evolutionary history of a species but for providing insight into a species' genetic diversity, which further enables proper assessment of the impact of various threats [5,6].The use of mitochondrial DNA (mtDNA) has been applied in molecular analyses to aid conservation studies to investigate the genetic diversity and genetic structure of a species and identify any evolutionarily significant units [7][8][9][10][11][12].Furthermore, phylogeographic studies are used to understand the role of paleogeographical events and demographic processes in shaping diversification within a species and the geographical distribution of the lineages [13].
Vipera ammodytes (Squamata: Serpentes: Viperidae: Viperinae) is a venomous viper with a distribu-tion ranging from southern Austria and northern Italy, through the Balkans to Asia Minor [14].This snake inhabits dry rocky habitats including woodland and scrub, sand dunes, hillsides, screes, stone walls, and sometimes cultivated areas with an elevation range limit of up to 2500 m a.s.l.[14,16,17].Seven subspecies of V. ammodytes were recognized in the past on the basis of phenotypic characters [14,18] before contemporary multivariate analyses of morphology synonymized three, and suggested the validity of only four subspecies, V. a. ammodytes, V. a. meridionalis, V. a. montandoni and V. a. transcaucasiana [19].
Using a phylogeographic approach, seven main genetic clades within V. ammodytes were defined [15] as follows: (1) a northwestern clade that inhabits Italy, Austria, Slovenia, Croatia and Bosnia and Herzegovina; (2) a northeastern clade from western, central and eastern Serbia and western Bulgaria; (3) a Montenegrin clade from Montenegro; (4) a southwestern clade from Albania and northwestern Greece; (5) a southeastern clade from the southernmost part of Serbia, most of the Republic of North Macedonia, central and eastern Bulgaria, northern and central Greece, Turkey, Armenia and Georgia; (6) a Peloponnese clade from the Peloponnesian peninsula; (7) a Cyclades clade from the Cyclades islands.However, the distribution boundaries of most of these clades have not been clearly defined due to the lack of samples throughout the range and from boundary zones of the clades, and because of the absence of nuclear markers.More information on the contact zones would be of great importance for our understanding of the current distribution and evolutionary integrity of these lineages.
The main aim of this study was to identify the position of some populations of V. ammodytes in Serbia from localities that have not been included in previous phylogeographic studies and are geographically positioned in proposed delineation zones of the main clades [15].For analysis, we used mtDNA sequences from the 16S rRNA gene as several sequences of the same gene were already available for this species [15]; this is a common and reliable universal barcoding marker used for vertebrate studies [20].

Ethics statement
All samples were collected with permits obtained from the Ministry of Agriculture and Nature Protection (No. 353-01-170/2016-17 and No.353-01-2666/2016-17).At present, vipers are not protected by law in Montenegro, while the sample from the border between Albania and the Republic of North Macedonia was collected opportunistically from a road-kill specimen.

Sample collection
This study was based on a combination of 16S rRNA sequences from novel samples and sequences obtained from GenBank.In total, we analyzed 47 16S rRNA sequences.Tissues were collected from tail tips (preserved in 96% ethanol) from 26 individuals sampled at 15 localities in Serbia, one locality in Montenegro and one locality in the Republic of North Macedonia directly bordering Albania.The samples in Serbia were mostly collected during 2016 and 2017, whereas samples from Montenegro and Albania/North Macedonia were collected in 2010.Twenty-one 16S rRNA sequences submitted to GenBank by Ursenbacher et al. [15] were included.GenBank accession numbers for previously published sequences and from this study with the localities for collected samples are given in Supplementary Table S1.Four additional 16S rRNA sequences from GenBank of four different species of vipers (Vipera latastei, Vipera aspis, Montivipera xanthina and Vipera berus) from GenBank were used as outgroups in the phylogenetic analysis.

Phylogenetic analysis
Sequences obtained in this study were visually inspected by FinchTV (Geospiza, Inc., Seattle, WA) and aligned with published sequences of the same mtDNA region using ClustalW [22].Bayesian phylogenetic inference was performed with MrBayes software v.3.2 [23], using TPM3uf (three parameter model) [24] + I model, which was selected by the JmodelTest program v.2.1.4[25].A Markov Chain Monte Carlo (MCMC) search was made for 3x10 6 generations using tree sampling every 100 th generation; 25% of the trees were discarded as 'burn-in' before generating a consensus tree.Consensus trees were viewed and edited in FigTree v. 1.4.0.Genetic distances between novel sequences were calculated using neighbor-joining method implemented in PAUP v. 4.0b10.

RESULTS
The analysis of 26 sequences of the 16S rRNA gene revealed eight unique haplotypes.In Serbia, six haplotypes were detected: the first was shared by 15 samples, the second was shared by two (15KRP_RS and 3SVL_RS), and each of the remaining detected haplotypes were unique to one sample each (19RKO_RS, 24POZ_RS, 26LJU_RS and 27MGT_RS).The remaining two haplotypes represented samples from Albania/ North Macedonia (21ALM_MK) and another from Montenegro (18MON_ME).All samples sequenced in this work had a deletion at position 2206 nt of mtDNA in the reference sequence NC 036956, except for the sample from Montenegro (18_MON_ME), which had thymine at the position, as did a few other Gen-Bank sequences (MO6_ME, MO2_ME, SW3_GR, SW2_GR, SW1_AL, SE5_BG and SE1_RS, SE3_AR, SE4_TR, CY4_GR, CY3_GR, CY2_GR, PE2_GR, PE1_GR, PE3_GR).
The phylogenetic tree (Fig. 1) was concordant with the previous study on V. ammodytes [15] regarding the main genetic clades.In two localities (in western and central Serbia) we confirmed the overlapping distribution of two distinctive clades (northeastern and northwestern).One of the two individuals belonging to the same local population in western Serbia (code 4KRP_RS in Fig. 1) clustered with the northeastern clade, the other (code 15KRP_RS) clustered with the northwestern clade.The genetic distance between these two vipers was estimated as 0.01 (1%).At the other locality from central Serbia, six of the seven analyzed individuals from the same local population clustered with the northeastern clade, while individual 3SVL_RS clustered with the northwestern clade (see Fig. 1).The genetic distance between the "northwestern" individual (code 3SVL_RS) and another six from that same population was estimated as 0.01 (1%).All the other novel samples from Serbia clustered with the northeastern clade (Fig. 1).The specimen from Montenegro (code 18MON_ME) clustered with the Montenegrin branch of the northeastern clade.The nose-horned viper collected on the Albanian-North Macedonian border (code 21ALM_MK) clustered with the northeastern clade.We have presented the geographical positions of the identified genetic clades for the samples analyzed in this study in Fig. 2, together with the surrounding sequences from GenBank showing their affiliation to the specific clade.

DISCUSSION
We recognize the presence of three main V. ammodytes mtDNA-defined clades in Serbia instead of the two that were previously identified [15]: (1) populations from most of Serbia are part of the northeastern clade; (2) in the extreme southeast of Serbia, the populations are part of the southeastern clade; (3) we have also identified a couple from the northwestern clade in western and central Serbia, where they overlap with the northeastern clade.Western Serbia has been mentioned [26] as part of the eastern border for the northwestern clade, and it was assumed that specific orographic formations and aquatic barriers prevented the further spread of the northwestern clade east and southeast into the Balkans.It is still not clear which barriers contributed to the distribution of these two clades in Serbia.The detec- tion of populations in the westernmost part of Serbia that have individuals clustering with the northwestern clade suggests that the Drina River, which represents the natural border between Bosnia and Herzegovina and Serbia, has not been a strong enough barrier to isolate the northwestern clade in the west.
We assume that the presence of the northwestern clade in western and central Serbia is probably due to a natural transition zone between the two clades.However, we must take into consideration information obtained from our interviews with local inhabitants at the location in central Serbia where one of the seven samples clustered with the northwestern clade (3_SVL_RS); two farmers stated that 10 to 20 years ago there was a deliberate release of V. ammodytes collected from different localities in former Yugoslavia for the purpose of venom milking and production of anti-venom components.Therefore, the presence of genes corresponding to the northwestern clade could be the result of persisting haplotypes originating from individuals released outside of their place of origin.
Another contribution of this study is the more detailed understanding of the geographic distribution of the main genetic clades of V. ammodytes in the Central Balkans.Our results suggest that the northeastern clade extends further southward into the Balkans than previously defined [15].According to the molecular clock [15], the first splits within V. ammodytes occurred during the early Pliocene (3.4-4.9Mya), separating the Montenegrin and northeastern clade from the northwestern, the southwestern from the Peloponnese, Cyclades and southeastern clades.However, the authors could not determine where the split occurred.There were extensive inflows of the sea during the Pliocene at the area of the Axios and Strymon river basins that reached northwards almost to the Central Balkans along the Vardar and perhaps the Morava River valleys [27].The presence of the Dinaric, Pindus and Rhodope massifs could explain the distribution of the northeastern and the southeastern clades on the territory of North Macedonia.These orographic formations and aquatic barriers isolated the southeastern, the southwestern and the northeastern clades in this area.As postulated for V. ammodytes [15], Ichthyosaura alpestris [28] and Lissotriton vulgaris [29], the distribution of the southwestern clade corresponds to regions of the "refugia-within-refugia" isolated by the Dinaric and Pindus Mts.Nevertheless, as there is a lack of samples from most of the territory of Kosovo, North Macedonia and Albania, in order to clearly define the boundary zones between the northeastern and the southeastern clades, further analyses of additional populations from these areas are needed.
The precise shape of the western boundary of the Montenegrin clade has not been defined [15] and implies the necessity of analyzing populations in the area extending from Dubrovnik in the north to the border between Montenegro and Croatia in the south.Sampling should also include populations distributed from the border zone between Bosnia and Herzegovina and Montenegro north of the Tara River canyon and further west to the Neretva River, which is known to be one of the biogeographic barriers in this part of the Balkan Peninsula [30].
Further investigations include an analysis of other nuclear genetic markers apart from mtDNA, in order to identify the level of gene flow between these clades as it is possible that nuclear gene exchange does occur [31].This species has low dispersal capability.Together with the results of the investigation of Ursenbacher et al. [15] based only on mtDNA, it would appear that sex-biased dispersal contributed to the obtained pattern of haplotype distribution.If males tend to disperse further than females during the mating season, which is expected and confirmed for most snake species [32][33][34], genetic analysis of nuclear genes would be necessary for further clarification of the contact zones between clades.To explain the distribution of nose-horned viper haplotypes in the studied area, we should not ignore the possible impact of occasional human-mediated transport.
The results of this study support future conservation plans for the nose-horned viper in the Central Balkans.This species is legally protected in Serbia because of uncontrolled overexploitation in the past [35], and it is still subject to illegal collecting at one locality (personal communication from a local collector).Even sustainable exploitation can lead to alteration and possible loss of genetic variation [36], and maintaining a species' genetic diversity in areas where it is exposed to illegal collecting is crucial [2].
To conclude, the results of this study provide insight into the phylogenetic positions of nose-horned viper populations that inhabit the boundary areas of mtDNA-identified clades in the Central Balkans.Our findings underline the importance of analyzing multiple samples from the same local population in genetic studies, especially if those populations are located in the supposed transition zones between distinctive genetic clades.To obtain a more complete insight into species' or population genetic structures, it is necessary to conduct comprehensive analyses using both mitochondrial and nuclear genes and to include a range of samples from locations in the supposed contact zones of the main genetic clades.

Fig. 2 .
Fig. 2. Map of the sampling localities of identified genetic clades (this study) with surrounding samples used from GenBank.Multiple samples from the same location that belong to the same clade are represented with a single mark for the following: NE2 (6), NE5 (2), NE6 (2).Novel samples are marked with stars while the used GenBank samples are marked with squares.NE corresponds to the northeastern clade, NW corresponds to the northwestern clade, MO corresponds to the Montenegrin clade, SE corresponds to the southeastern clade, SW corresponds to the southwestern clade; these abbreviations correspond to the abbreviations of clades in Fig. 1.Shaded areas suggest the updated potential geographic distribution of the genetic clades in the Central Balkan Peninsula.

Funding:
The fieldwork was kindly funded by The Rufford Foundation (Grant Nos.19578-1 and 23392-2) for TČ, and by the Ministry of Education, Science and Technological Development of the Republic of Serbia, Grant No. 173025 for both JCI and TČ.Genetic analyses were supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia, Grant no TR37009.