Evaluating taxonomic inflation: towards evidence-based species delimitation in Eurasian vipers (Serpentes: Viperinae)

The designation of taxonomic units has important implications for the understanding and conservation of biodiversity. Eurasian vipers are a monophyletic group of viperid snakes (Serpentes, Viperinae), currently comprising four genera (Daboia, Macrovipera, Montivipera and Vipera) and up to 40 species. Taxonomic units have been described using a wide variety of methods and criteria, and consequently, considerable controversy still surrounds the validity of some currently listed species. In order to promote a consensusand evidence-based taxonomy of Eurasian vipers, we analysed published mitochondrial and nuclear DNA sequences for this group to reconstruct phylogenetic relationships among currently recognized viper species. We also compiled information on external morphology to assess their morphological distinctiveness. Phylogenetic inference based on mtDNA sequences shows contrasting levels of divergence across genera and species and identifies several instances of non-monophyly in described species. Nuclear DNA sequences show extremely low levels of genetic variation, with a widespread pattern of allele sharing among distant species, and even among genera. Revision of morphological data shows that most species designations rely on scalation traits that overlap extensively among species of the same genus. Based on our combined assessment, we recognize 15 taxa as valid species, three taxa which likely represent species complexes, 17 taxa of doubtful validity as species, and five taxa for which species status is maintained but further research is highly recommended to assess taxonomic arrangements. We stress the need to implement integrative taxonomic approaches for the recognition of evidence-based taxonomic units in Eurasian vipers.


Introduction
The designation of taxonomic units has important implications for the way we study, describe and understand biodiversity, as well as for how we mobilize efforts and allocate resources to develop conservation strategies. Over the years, different criteria and tools have been used to define species, leading to a succession of species concepts that resulted in extended controversy within the research community (Mayden, 1997;de Queiroz, 2007). Nowadays, a species is often defined as a separately evolving metapopulation lineage that possesses relevant characteristics that allow assessing its distinctiveness from others (i.e., the unified species concept; de Queiroz, 2007, and its precursor, the evolutionary species concept; Simpson, 1961;Wiley, 1978;Frost and Hillis, 1990). This definition is linked to the integrative taxonomy framework, which is based on the combination of different lines of evidence (e.g., genetic, morphological, ecological) and methodologies (e.g., phylogenetic inference, ordination methods, ecological modelling) to objectively identify taxa (Dayrat, 2005) that -in an ideal case -would represent independently evolving species. Names of species are therefore intended to identify biologically cohesive populations with recent common ancestry rather than to recognize unusual patterns of distribution or morphology (Kaiser et al., 2013).
Eurasian vipers are a monophyletic group within the subfamily Viperinae (Serpentes, Viperidae), whose members are distributed primarily in the Palaearctic region, i.e., nontropical Eurasia and North Africa (Phelps, 2010). This group is phylogenetically sister to a clade of Middle Eastern vipers, constituted by the genera Eristicophis and Pseudocerastes (see Phelps, 2010;Zheng and Wiens, 2016), which are not considered in this work. At the time of writing, the most recent and comprehensive list of reptiles (i.e., The Reptile Database; Uetz, Freed and Hošek, 2019) lists four genera and 40 species within Eurasian vipers (table 1): Daboia, with 4 species; Macrovipera, with 3 species; Montivipera, with 8 species; and Vipera, with 25 species. However, Eurasian vipers have a long taxonomic history, and different authors have used a wide variety of methods and criteria to define taxonomic units, as reflected in previous species lists (e.g., Mallow, Ludwig and Nilson, 2003;Phelps, 2010).
At the genus level, the history of Eurasian vipers is relatively simple. Through most of the 20th century, all species considered here were included in the single genus Vipera (Boulenger, 1896(Boulenger, , 1913Schwarz, 1936;Klemmer, 1963;Minton, Dowling and Russel, 1968). The maverick German herpetologist Albert Franz Theodor Reuss described numerous genera within the Eurasian vipers (reviewed by Krecsák, 2007), but these gained little traction with subsequent authors, except where the names had priority for subsequently validated clades. Obst (1983) was the first author to challenge the monogeneric classification of Eurasian vipers, by separating the larger species into the genus Daboia, together with Pseudocerastes. This split however, was not adopted by most subsequent researchers. Herrmann, Joger and Nilson Table 1. Currently listed species of Eurasian vipers according to Uetz, Freed and Hošek (2019), depicting the criteria used for species designation, our recommended status, in some cases provisional, including justification, suggestions for further work and IUCN red list category (also for included taxa). Percentages of divergence from the closest sister species or clade are given for a small fragment of cyt b (196 bp; supplementary table S3). "DVAS" means Doubtful Validity As Species, may represent geographic variation, a subspecies or diverged population, hence, currently we recommend to decline its species status; "LSC" means Likely Species Complex (i.e., group of closely related taxa, possibly more than one species, but delimitations are not clear yet), hence, currently maintain species status; and "pending" indicates that, despite some incongruences with the divergence threshold delimitation, single-species recognition is currently maintained, but further research is required to assess taxonomic status. Daboia palaestinae (Werner, 1938) morphological species high genetic divergence (15%) LC Daboia russelii (Shaw and Nodder, 1797) morphological species high genetic divergence (9%) Daboia siamensis (Smith, 1917) (Reuss, 1933) morphological DVAS low to moderate genetic divergence to V. renardi (2-4%); it includes V. ebneri, which was previously synonymized and shows low genetic divergence to V. eriwanensis (1%); geographic isolation doubtful (Rajabizadeh et al., 2011, Tuniyev et al., 2018; private haplotypes in some nuDNA genes; a subspecific status must be appropriate.
nuclear inferences should be increased Vipera graeca (Nilson and Andrén, 1988) Nilson et al. (1999) described Montivipera as a new subgenus for the xanthina group. This was subsequently raised to full genus level by Joger (2005). Lenk et al. (2001), using mitochondrial DNA sequences, assigned the species mauritanica, deserti and palaestinae to Daboia, leading to the current generic arrangement of the group. The recent use of subgenus Pelias (Merrem, 1820) as a full genus for the berus and ursinii groups (e.g., Avcı et al., 2010;Tuniyev et al., 2012Tuniyev et al., , 2013Tuniyev et al., , 2018a, on the other hand, has remained a minority opinion in the literature. Certain changes proposed outside the peer-reviewed scientific literature are not considered here for reasons given in Kaiser et al. (2013). One third of the currently recognized species were described in the 18th and 19th centuries by recognized taxonomists and zoologists of the time (e.g., Carl Linnaeus described Coluber lebetinus, C. berus, C. aspis and C. ammodytes (now in Macrovipera and Vipera) in 1758; John Edward Gray described Clotho mauritanica and Daboia xanthina (now in Daboia and Montivipera, respectively) in 1849; Eduard Boscá described Vipera latastei in 1878; table 1). During the 20th century, eighteen of the currently recognized species were described; ten of them were described before or during the 1970s, and eight after that period. All species descriptions were based on morphological traits (i.e., scale counts, biometric measures, colour patterns), with the application of statistical analyses of these morphological and other phenotypic traits gradually becoming incorporated during more recent times (e.g., Herrmann, Joger and Nilson, 1992;Nilson and Andrén, 2001).
Recent phylogenetic and phylogeographic studies have transformed the taxonomic panorama in Eurasian vipers considerably, validating some taxa as species (e.g., Vipera graeca, Mizsei et al., 2017), rejecting or synonymising others (e.g., Vipera altaica with V. renardi, Zinenko et al., 2015 Zinenko et al., 2016), or modifying previously designated taxonomic units (e.g., assigning species to four genera; Garrigues et al., 2005). However, despite this multitude of studies, there is still considerable uncertainty regarding the validity of some species, which is hampering the development of optimized conservation strategy for the whole group.
Here, we apply an integrative approach to review the taxonomy of the Eurasian vipers, by bringing together and analysing existing molecular and morphological data. We compiled and analysed published and new mitochondrial and nuclear DNA sequences for this group to reconstruct phylogenetic relationships among currently recognized species. We also compiled information on external morphology, as well as on the criteria used for species delimitation, in order to assess the morphological distinctiveness of currently recognized species. Our objectives are: 1) to provide an updated taxonomy of Eurasian vipers; 2) to evaluate the validity of some of the newest species designations under the unified species concept and recommend an appropriate status reflecting our current data base; 3) to identify remaining knowledge gaps and the research required to achieve a robust, stable and evidence-based taxonomic framework for this group of vipers.

Taxonomic inference
Our taxonomic evaluation is built upon an integrative and evolutionary framework, based on the unified species concept (de Queiroz, 2007), under which species are defined as separately evolving lineages and their biological properties (e.g., monophyly, reproductive isolation, differentiated ecological niches, morphological distinctiveness) as "operational criteria" that ultimately provide evidence for their separation through a variety of methods. We employ the integrative approach of Padial et al. (2010) by using phylogenetic analysis of mitochondrial (mtDNA) sequences, the most widely available and standardised marker, to identify divergent lineages, that are then tested, to the extent that the available data allow, with additional data, in particular single-copy nuclear DNA (nuDNA) markers and morphological data. We set a sequence divergence percentage to propose a threshold for provisional taxonomic categorization (see phylogenetic inference section for more details). Above this threshold, evolutionary lineages may simply confirm established species, or if not described as such, they may be considered as candidate species and should be targeted in future studies to evaluate their taxonomic status. If currently recognized species are composed of several divergent lineages, we classify them as a Likely Species Complex (LSC). Below the cyt b threshold, currently recognized species are categorized as Doubtful Valid as Species (i.e., DVAS), until there is enough evidence to indicate species integrity. In addition, we used the "pending" category to indicate that single-species recognition is currently maintained, despite incongruences with the divergence threshold delimitation, but further research is recommended to assess taxonomic arrangements.

Phylogenetic inferences
We searched on GenBank for mtDNA and nuDNA gene sequences representing all the relevant lineages within Eurasian vipers (see supplementary table S1 and S2). Selection of sequences was based on published phylogenetic and phylogeographic studies (e.g., Ursenbacher et al., 2006Ursenbacher et al., , 2008aVelo-Antón et al., 2012;Zinenko et al., 2015;Stümpel et al., 2016;Mizsei et al., 2017;Freitas et al., 2018;Martínez-Freiría et al., 2020) to ensure that the selected sequences represent the genetic structure and units reported in those studies. Mitochondrial DNA sequences were available for a total of nine gene fragments, as well as the whole mitochondrial DNA genome for D. siamensis (mislabelled as D. russelii in GenBank). For phylogenetic analysis, we selected DNA fragments from a subset of seven mtDNA markers with a higher representation across species: CR (control region), COI (Cytochrome c oxidase subunit I), cyt b (cytochrome b), ND2 (NADH dehydrogenase subunit 2), ND4 (NADH dehydrogenase subunit 4), ND5 (NADH dehydrogenase subunit 5) and 16S (mitochondrial gene coding for 16S rRNA).
Additionally, we also provide unpublished DNA sequences generated for other studies to complement available data from GenBank. Sequences were concatenated when they originate from the same lineage or geographic locality, using SequenceMatrix software (Vaidya, Lohman and Meier, 2011). Since Eurasian viper species generally contain geographically cohesive, parapatric mitochondrial lineages, it is unlikely that they generate chimeras when concatenating sequences from multiple individuals from the same geographic locality. The final dataset included 97 units (representing several lineages for 39 of the 40 recognized species) with sequences ranging from 654 to 4621 base pairs (bp). Details of sequences used in mtDNA analyses are available in supplementary table S1. Sequences for nuDNA were available for 17 gene fragments, from which we selected six protein-coding genes, BDNF (Brain Derived Neurotrophic Factor), CMOS (Oocyte maturation factor Mos), MC1R (Melanocortin 1 Receptor), NT3 (Neurotrophin-3), PRLR (Prolactin Receptor), RAG1 (Recombination Activating protein 1) and one intron, B-fib (Beta-fibrinogen intron 7), as they allowed the most comprehensive taxonomic coverage (24 species, 61.5% of the currently recognized total). This set includes the most relevant nuclear genes used to support the designation of some of the most recently recognized species (e.g., RAG1, Ghielmi et al., 2016;NT3, Mizsei et al., 2017). Details on nuDNA sequences used are provided in supplementary table S2. Sequences were manually aligned and edited using Geneious v 4.8.5 (Kearse et al., 2012). For the nuclear genes, haplotype phases were produced by a coalescent-based Bayesian reconstruction implemented in PHASE (Stephens, Smith and Donnelly, 2001) available in DNAsp (Librado and Rozas, 2009).
Phylogenetic relationships and time of divergence between species were inferred using a Bayesian Inference (BI) method implemented in BEAST v 1.7.5  on the concatenated mtDNA dataset. An exhaustive search with PartitionFinder 1.1.1 (Lanfear et al., 2012) was conducted to select appropriate partitioning schemes and evolutionary models based on the Bayesian Information Criterion (BIC). The GTR + G + I model applied to all mtDNA fragments combined in a single partition was determined as the best-fit model and partitioning scheme.
Substitution rates were estimated under a strict molecular clock (Drummond et al., 2006) that assumes uniform rates across branches. A Yule model, most suitable for specieslevel phylogenies, was implemented as tree prior. Since the fossil record of Eurasian vipers is fairly incomplete and does not provide reliable and verified calibration dates (see Stümpel et al., 2016), our molecular dating strategy relied on secondary calibrations, including the splits of Vipera-Daboia and Macrovipera-Montivipera, dated at 26 Mya (Zheng and Wiens, 2016 -but see Šmíd and Tolley, 2019). We used a lognormal prior with a mean of 26.3 Mya and a standard deviation of 0.07 to constrain node ages. Three independent runs of 100 million generations were performed, sampling trees and parameter estimates every 10 000 generations with 10% of the trees discarded as burnin. Convergence was verified by looking at the effective sample sizes of all parameters (ESS > 300) using Tracer v1.7 (Rambaut et al., 2018). Trees obtained from multiple independent runs were then combined using LogCombiner v 1.7.5.  and summary trees were generated with TreeAnnotator v1.7.1 .
Haplotype networks for the cyt b gene (the most widely used marker across studies) were reconstructed with TCS v 1.21 (Clement, Posada and Crandall, 2000), with a 90% parsimony connection limit; and the graphical output was visualized in TCSBU (dos Santos et al., 2015). The initial alignment of 1141 bp was trimmed down to 196 bp with no missing data across 96 units. The short length of the single sequence available for V. shemakhensis precluded its inclusion in the haplotype network. However, it was directly compared to the remaining sequences (180 bp of overlap). Uncorrected p-distances between taxa were estimated based on the same cyt b fragment using MEGA ver. 5 (Tamura et al., 2011). We propose the use of this short standardized cyt b fragment as a candidate DNA barcode for the delimitation of evolutionary units in Eurasian vipers (see Hebert and Grégory, 2005). Furthermore, we recommend a value of 5% cyt b sequence divergence as a provisional threshold, below which an untested species-level designation appears inappropriate or premature, unless other lines of evidence would validate species classification. The cyt b-threshold results from two facts: 1) a 5% threshold provides a good delimitation of currently recognized species within Eurasian vipers, also recognizing deep evolutionary lineages within them; 2) cyt b divergence levels are consistently higher than 5% between closely related species that co-exist with restricted or no hybridization, e.g. sympatric species of the genus Vipera differ by more than 10% (Tarroso et al., 2014;Mebert et al., 2015b), sympatric watersnakes Nerodia fasciata and N. sipedon differ by 9% (Mebert, 2008), Natrix helvetica and N. natrix by 6.9% (Kindler et al., 2017) and 9% difference between Montivipera wagneri and M. raddei with no signs of mixing along a sharp contact line (Mebert et al., 2016;Stümpel et al., 2016).
Haplotype networks were drawn for each of the seven nuclear genes following the same procedure. Sequences with high proportion of missing data (>30% of the total length) were excluded from the dataset. For B-fib and RAG1, datasets were divided in two sets each based on sequence length (B-fib-1, B-fib-2 and RAG1-1, RAG1-2) and analysed independently to avoid excluding shorter sequences from the haplotype networks.

Morphological characterization
From the published literature, we compiled a list of the criteria used to identify each of the currently recognized species. In addition, we gathered information on the variability of 14 external morphological traits in each species, including one biometric, 11 pholidotic and two dorsal colour pattern characters. Information from the literature was collated with data from specimens measured by the authors in the field or in museum collections. This allowed us to establish the range of variation of different morphological traits, in some cases for each sex separately (see Results). In addition, whenever possible, and in order to more accurately represent morphological variation in some pholidotic traits, modal or mean values were retrieved for specific groups or subspecies within each species. Furthermore, we included verbal descriptions of colouration, which provide an idea of variation in visually striking qualitative traits.
Vipera is the most diverse genus and includes three subclades: (1)   Lineages are grouped and distinctively coloured considering a divergence equal or higher than 5% for a small fragment (196 bp) of cyt b (supplementary table S3). Black and light-grey dots represent posterior probabilities higher than 0.95 and between 0.9-0.95, respectively. Names of taxa and lineages are given accordingly to publications from where sequences were retrieved. See supplementary table S1 for details.
Haplotype networks based on the 196 bp cyt b fragment recover similar genetic relationships to the mtDNA phylogenetic tree (supplementary fig. S2). The 82 identified haplotypes and matrices of uncorrected genetic distances based on these cyt b fragments are provided in supplementary material (supplementary tables S1 and S3, respectively). The single sequence available for V. shemakhensis differed in one position from the sequences of V. ebneri and V. eriwanensis.

Phylogenetic inferences from nuDNA
Haplotype networks constructed for the seven nuclear genes show a pattern of wide haplotype sharing and few mutational steps among species; however, some species present distinct, well differentiated haplotypes ( fig. 2;  Haplotype networks for B-fib, PRLR and NT3 present higher variability than the other selected markers ( fig. 2). The B-fib shows unique, well-separated haplotypes (more than three mu-

Morphological characterization
We obtained data from 39 studies providing morphological descriptions or addressing the variability of species, 13 assessing phylogenetic relationships among species, and six using phylogenetic inferences and phenotypic characterization to describe taxa (see supplementary references S1). One book (Phelps, 2010) was used to extract maximum body size for some species for which other published data were lacking. Information retrieved from publications, complemented with data collected by the authors from museum specimens, is shown in supplementary table S5.
The most commonly used criteria to diagnose species relied on head (24 species) and body (26 species) pholidosis, and dorsal colouration (27 species). In contrast, phylogenetic analyses based on molecular data were initially considered for species description in only five cases (table 1, supplementary table S5).
Data compilation showed that species within Macrovipera and Daboia are the largest in body size, followed by Montivipera and with Vipera being the smallest. Sexual dimorphism in body size is reported in most species (supplementary  table S5).
Regarding pholidosis, some traits of head scalation (e.g., canthal, supralabial or infralabial scales) exhibit low variation, particularly within each genus (supplementary table S5). Modal or mean values of other head traits, however, show important variation within genera (e.g., apical and intercanthal + intrasupraocular scales within Vipera) and among them (e.g., loreal scales; fig. 3). Modal or mean values of body scalation exhibit variation among species (e.g., ventral and subcaudal scales; fig. 5) or genera (e.g., number of dorsal rows; fig. 4). Nevertheless, the ranges of variation for most pholidotic traits overlap extensively among species of the same genus and even among distinct genera (figs 3, 4). Sexual dimorphism in subcaudal scale counts is also mirrored in ventral scale variation in many species (fig. 4).
With respect to dorsal colouration, there is high variability at both the inter-and intraspecific levels, particularly in Vipera species, which display more dorsal marks and a higher number of distinct dorsal pattern types than species of the other genera (supplementary table S5).

Discussion
Eurasian vipers are a taxonomically challenging group due to the long-standing use of diverse species criteria and the prevalence of morphology-based classifications and species delimitations. Here, we inferred the molecular phylogeny of Eurasian vipers based on mitochondrial data, encompassing for the first time almost all the described species and evolutionary lineages within this group. We also gathered information on the morphological variability of currently recognized species and reviewed the criteria used to identify them. By assessing phylogenetic and morphological variability, we provide recommendations and future research directions for robust species delimitation, which will aid the advancement towards a more informed and coherent taxonomy for this group of vipers.

Phylogenetic inference from mtDNA
Phylogenetic analyses based on seven mtDNA fragments produced a mostly resolved topology that strongly supports the monophyly of the four genera ( fig. 1). Phylogenetic relationships and divergence dates estimated in this study mostly agree with those reported in previous works (e.g., Alencar et al., 2016;Stümpel et al., 2016;Zheng and Wiens, 2016), with the exception of divergence dates recently reported using alternative time calibration procedures (Šmíd and Tolley, 2019).
Overall, our inferences suggest distinct diversification dates and divergence levels for each genus. Diversification within Daboia is estimated in the early Miocene, followed by the diversification within Vipera in the middle Miocene; Macrovipera and Montivipera appear as the most recent genera, diverging from each other in the late Miocene ( fig. 1). Despite its old origin, Daboia shows low diversity levels, with only four species recognized: two tropical Asian species, D. russelii and D. siamensis, and two Mediterranean species, D. mauritanica and D. palaestinae. These species, however, are extremely divergent, with genetic distances based on a small (196 bp), variable fragment of cyt b ranging from 9%, between the Asian relatives D. russelii -D. siamensis, to 21%, between D. siamensis and the Mediterranean D. palaestinae (supplementary table S3). The old diversification and wide distributional range of this genus, together with the spatial gaps separating species ranges, suggest that this group likely experienced major extinctions along its evolutionary history, with only some of the representative taxa currently persisting. This suggestion is supported by the occurrence of Daboia-like vipers  With a more recent origin, Macrovipera also shows low diversity levels, comprising three described species ( fig. 1): M. schweizeri, from the western Cyclades, M. lebetina, widely distributed from central Asia to the Middle East, and M. razii from southern and central Iran. High levels of polymorphism, especially in colouration, led to the description of distinct subspecies within M. lebetina (Nilson and Andrén, 1988;Stümpel and Joger, 2009). Molecular studies have suggested the validity of four subspecies, i.e., lebetina, obtusa, turanica and cernovi (Stümpel and Joger, 2009), but rejected the species status of M. schweizeri, which is suggested to be a subspecies of M. lebetina with a wider distribution than previously considered (Lenk et al., 2001;Stümpel and Joger, 2009). Our phylogenetic analyses recover M. schweizeri and M. lebetina as a single unit (2% of cyt b genetic distance; supplementary table S3). Shared mitochondrial haplotypes between Macrovipera schweizeri and some southern Turkish populations of M. lebetina further corroborate conspecificity (Stümpel and Joger, 2009). Oraie et al. (2018) found cytochrome b divergences between moderate and high (up to 4.4%) among Iranian populations of M. lebetina obtusa and M.l. cernovi and between these and specimens from Uzbekistan and Azerbaijan, suggesting that further diversity may exist within this species. Moreover, increased geographical sampling may uncover additional phylogenetic diversity, as has been the case with the description of M. razii from Iran, previously allocated to M. lebetina (Oraie et al., 2018).
The genus Montivipera consists of two wellsupported species complexes, the xanthina and raddei clades. This genus was recently subjected to a comprehensive phylogenetic study based on a multilocus mitochondrial and nuclear dataset for all its constituent taxa (Stümpel et al., 2016). Here we confirm previous findings of a more recent origin of the raddei-complex (ca. 1.68 Mya, 2% genetic divergence within the group) and an older divergence with higher levels of genetic diversity within the xanthinacomplex, including the bornmuelleri-group (ca. 4.78 Mya, 6% genetic divergence within the entire complex). Major branches within the latter group are supported by high posterior probabilities, except for the monophyly of M. xanthina, for which low branch support was already shown in Stümpel et al. (2016).
Among all analysed genera, Vipera is the most diverse and well-studied group. Extensive phylogeographic and phylogenetic work has been conducted within this genus (e.g., Garrigues et al., 2005;Ursenbacher et al., 2006aUrsenbacher et al., , b, 2008Velo-Antón et al., 2012;Zinenko et al., 2015Zinenko et al., , 2016Martínez-Freiría et al., 2020). However, no previous study has included all known taxa (e.g., 12 taxa in Zheng and Wiens, 2016;19 taxa in Šmíd and Tolley, 2019). Although traditionally grouped in two subgenera (Garrigues et al., 2005), Vipera forms three well-supported monophyletic groups ( fig. 1): Pelias, Vipera 1 and Vipera 2 (see Nilson and Andrén, 1997 for a similar designation). Both Vipera 1 and Vipera 2 present deep phylogenetic structure and high divergence between and within taxa. Vipera 1 includes the western and Mediterranean V. aspis and V. latasteimonticola, and Vipera 2 is represented by V. ammodytes-transcaucasiana from the Balkans, Turkey and Asia Minor. V. latastei-monticola and V. ammodytes-transcaucasiana are species complexes that comprise highly divergent units (9% and 6% of cyt b genetic distance within each complex, respectively; supplementary table S3). Phylogeographic studies on these taxa show high levels of geographically structured genetic diversity and the oldest divergences among Eurasian vipers, with main diversification events likely occurring during late Miocene and early Pliocene (Ursenbacher et al., 2008 On the other hand, Pelias is the most diversified group in the phylogeny, with multiple described species. In general, genetic divergence within this group is shallower than that observed within the other Vipera subclades, with many taxa likely resulting from geographic splits during Pleistocene climatic oscillations (see Zinenko et al., 2015). Our phylogeny recovered two highly divergent subclades (TMRCA = 7 Mya, fig. 1; 9% genetic distance, supplementary table S3). One of the subclades includes V. berus (with V. barani nested within it) and V. seoanei. The other clade is highly diversified and includes three major groups: (1) V. renardi (with V. shemakhensis, V. lotievi and V. altaica nested within it) and V. eriwanensis (with V. ebneri) as sister group, together with V. ursinii and V. kaznakovi from Russia (V. orlovi, is an admixed population of V. kaznakovi and V. renardi, Zinenko et al., 2016; and V. dinniki is nested within V. kaznakovi from Russia), and V. graeca; (2) V. sakoi, V. kaznakovi from Georgia, V. darevskii (with V. olguni nested within) and V. walser; and (3) V. anatolica, which appears as a separate lineage at the root of this clade.
Our results support most taxa and recognized complexes (e.g., V. latastei-monticola) as monophyletic lineages, with the exception of V. kaznakovi and V. lotievi. As already shown by Zinenko et al. (2015) and Ghielmi et al. (2016), V. kaznakovi from Russia (Greater Caucasus) and Georgia (Lesser Caucasus) appear in the phylogeny as two polyphyletic lineages separated by a considerable genetic distance (5%), and V. lotievi is also polyphyletic, with several lineages included within V. renardi ( fig. 1, supplementary fig. S2). Multilocus RAD-sequencing data recovered low differentiation between both mitochondrial lineages within V. kaznakovi (Oleksandr Zinenko, unpublished data), a pattern that is further supported by low morphological differentiation and the existence of a continuous area of suitable habitats connecting these populations (Orlov and Tuniyev, 1990). Reasons for the discordance between mtDNA and nuDNA are still unclear but a historical occurrence of introgressive hybridization with asymmetric mitochondrial DNA capture could explain this pattern (e.g., Barbanera et al., 2009). The polyphyletic status of V. lotievi, on the other hand, is thought to be the result of possible confusion in the identification of species due to morphological convergence or hybridization and introgression leading to admixture of traits (Zinenko et al., 2015(Zinenko et al., , 2016.

Mito-nuclear discordance and low resolution of nuDNA
The nuclear data do not follow the same pattern observed for mtDNA. The haplotype networks constructed for each nuDNA marker show a pattern of widespread haplotype sharing among distant species, and even among genera, with extremely low levels of genetic variation ( fig.  2), suggesting incomplete lineage sorting of ancestral polymorphism (Wan et al., 2004) or extremely low levels of sequence evolution. Nuclear genes were already known to provide insufficient resolution to infer the phylogenetic variability within Eurasian viper taxa (e.g., V. latastei-monticola and V. aspis, Velo-Antón et al., 2012;Freitas et al., 2018;Martínez-Freiría et al., 2020;D. mauritanica, Martínez-Freíria et al., 2017a;Montivipera taxa, Stümpel et al., 2016). Yet, relying on mtDNA as the only source for phylogenetic inference could be problematic (see Ballard and Whitlock, 2004), and thus, the phylogenetic patterns retrieved here should be interpreted with some caution. The existence of local hybridization between highly differentiated species (e.g., Tarroso et al., 2014;Guiller, Lourdais and Ursenbacher, 2017) and the growing evidence for extensive gene flow in more recently differentiated taxa (Zinenko et al., 2015(Zinenko et al., , 2016 highlight the importance of nuDNA analyses. However, despite being commonly used in species delimitation studies, the sequenced nuDNA markers were too conserved to consistently differentiate between otherwise well supported species of Eurasian vipers. This applies even to introns and other fast-evolving single copy nuclear genes (e.g., NT3, PRLR, Townsend et al., 2008). Consequently, other molecular approaches such as even faster evolving markers (e.g., Kindler and Fritz, 2018;Pöschel et al., 2018), phylogenomic methods (e.g., Blair et al., 2019;Heinicke et al., 2018), and increased sampling to evaluate current or past gene flow between most-proximate populations or contact zones of two or more closely related species (Mebert, 2008(Mebert, , 2015a(Mebert, , b, 2020Hillis, 2019) should be favoured to infer evolutionary histories and resolve species limits among these taxa.

Taxonomic relevance of morphological traits
The value of morphological traits to define species in many groups has become questionable after the emergence of DNA-based methods in taxonomy. Morphological variation across populations often reflects local adaptation processes or phenotypic plasticity, rather than historical relationships (e.g., Kaliontzopoulou et al., 2011;Kaliontzopoulou, Carretero and Llorente, 2012;Alhajeri, Hunt and Steppan, 2015). Additionally, lack of morphological differentiation does not necessarily imply shared evolutionary histories, as is the case in cryptic species (e.g., Bickford et al., 2007;Kaliontzopoulou et al., 2011;Ghielmi et al., 2016) or taxa displaying convergent evolution (e.g,. Harmon et al., 2005).
The taxonomy of Eurasian vipers has been traditionally based on differences in pholidotic (head and body scalation) and colouration (dorsal pattern and colour) traits. Our data compilation shows that the ranges of variation of most pholidotic traits overlap extensively among species of the same genus (figs 4, 5), while colouration traits exhibit high variability at both the inter-and intraspecific levels (supplementary table S5). This apparent low prevalence of diagnostic traits for taxonomic purposes must be taken with caution given the limitations of our data compilation. While some species descriptions relied on few specimens from a small, not representative region and/or particular diagnostic traits (e.g., M. albizona Nilson, Andrén and Flärdh, 1990;M. kuhrangica Rajabizadeh, Nilson and Kami, 2011), robust studies using multivariate comparative analyses of morphological traits have shown morphological distinctiveness of some species (e.g., V. graeca, Nilson and Andrén, 2001;M. razii, Moradi, Rastegar-Pouyani and Rastegar-Pouyani, 2014), reinforcing the view that the combination of traits can in fact identify taxonomic units in some cases. On the other hand, the occurrence of morphologically-cryptic species, that have only been identified after molecular phylogenetic analyses (e.g., V. walser Ghielmi et al., 2016), suggests that other factors may be affecting the external morphological variability of Eurasian vipers.
Both pholidotic and colouration traits frequently display geographic variation associated with environmental gradients, reflecting adaptive processes (e.g., Shine, 2000;Sanders, Malhotra and Thorpe, 2004;Martínez-Freiría et al., 2009;Tomović, Crnobrnja-Isailović and Brito, 2010;Martínez-Freiría and Brito, 2013). The role of local adaptation in shaping intraspecific morphological differentiation has been highlighted in multiple studies on reptile species (e.g., Thorpe and Baez, 1993;Malhotra and Thorpe, 1997;Kaliontzopoulou, Pinho and Martínez-Freiría, 2018). In particular, traits related to fitness frequently present variation across different environmental and ecological conditions in order to meet the species-specific needs and enhance performance and fitness (Arnold, 1983;Kingsolver and Huey, 2003). In vipers, for instance, differences in dorsal pattern colouration can be an adaptive response to temperature gradients, enhancing thermoregulation capabilities, or to predation pressures, leading to aposematic signals or increased substratecrypsis (Wüster et al., 2004;Valkonen et al., 2011;Santos et al., 2014;Dubey et al., 2015;Martínez-Freiría et al., 2017). Increasing or decreasing scale numbers within species can influence water loss along environmental gradients (Malhotra and Thorpe, 1997;Sanders, Malhotra and Thorpe, 2004) or enhance locomotion over distinct substrates (Kelley, Arnold and Gladstone, 1997). These traits can also be plastic and depend on the thermal conditions experienced by the embryos during gestation. For example, both field and experimental studies indicate that thermal conditions at early embryonic stages influence the number of ventral scales, scale abnormalities, as well as dorsal colouration in vipers (Lourdais et al., 2004;Lorioux et al., 2013). Additionally, reproductive programs in low-effective-size populations have shown that inbreeding can lead to subsequent decrease of scale counts in viper offspring (Üveges et al., 2012). Altogether, these studies suggest that morphological traits might exhibit high interpopulation variability and little taxon-specific variation and thus, should no longer be used as the only source of data to inform taxonomic decisions in Eurasian vipers.

Taxonomic inflation and need for an integrative taxonomy
The taxonomy of Eurasian vipers has long been under intense debate, especially as it provides crucial underpinnings to the formulation of conservation management strategies and the allocation of economic resources for this purpose. Species have been described based on diverse criteria, and mostly using morphology and/or geographic isolation as the only source of inference. Not surprisingly, recent phylogenetic studies have often revealed major inconsistencies in relation to current taxonomic units (e.g., Ferchaud et al., 2012;Velo-Antón et al., 2012;Zinenko et al., 2015Zinenko et al., , 2016Stümpel et al., 2016;Martínez-Freiría et al., 2017), also unveiling the existence of morphologically cryptic taxa (e.g., Ghielmi et al., 2016). Yet, many recent species descriptions maintained the traditional methods of species delimitation, disregarding known limitations (e.g., Vipera altaica Tuniyev, Nilson and Andrén, 2010; Montivipera kuhrangica Rajabizadeh, Nilson and Kami, 2011;Vipera olguni Tuniyev et al., 2012;Vipera shemakhensis Tuniyev et al., 2013). In agreement with previous molecular studies, our phylogenetic reconstruction shows a clear mismatch between relevant evolutionary units and recognized species. This is particularly evident within Vipera, in which species complexes such as V. latastei-monticola or V. ammodytestranscaucasiana include highly divergent lineages from the Miocene, while other species are much younger (from late Pleistocene), polyphyletic (e.g., V. lotievi) or are nested within others (e.g., V. altaica, V. shemakhensis). This puzzling scenario mainly results from the indiscriminate use of morphological traits and geographic isolation as exclusive criteria for species delimitation (see table 1). Polymorphism as a result of local adaptation and plasticity may lead to taxonomic inflation, whereas low morphological variability can hamper the identification of cryptic species, such as in the case of the V. latastei-monticola and V. ammodytestranscaucasiana complexes. One example that illustrates well the problem of high morphological variation within species is the case of V. lotievi, for which morphological convergence across similar environments and confusion over species identification are highlighted as a possible explanation for its polyphyly (Zinenko et al., 2015). Similarly, the occurrence of ecotypes can lead to the designation of taxonomic units which are not concordant with evolutionary history (e.g., V. aspis atra, Ursenbacher et al., 2006b; V. aspis montecristi, Barbanera et al., 2009;Luiselli et al., 2015;V. monticola, Velo-Antón et al., 2012). Additionally, genetic introgression was also shown to be a confounding factor on species classification within the Pelias subgenus, leading to intermediate phenotypes in admixed populations (e.g., the description of V. magnifica and V. orlovi; Zinenko et al., 2016).
Due to these limitations, an increasing number of studies have applied molecular methods to inform taxonomic decisions. The amount of sequence divergence among groups has been extensively used as a criterion to delimit taxa (e.g., in small mammals, Bradley and Baker, 2001;in rat snakes, Hofmann et al., 2018). Although the use of a standard percentage of sequence divergence to separate species is debatable, in this work, we draw an arbitrary but conservative (in terms of recognizing most described species) threshold for the delimitation of evolutionary units in Eurasian vipers corresponding to an uncorrected genetic distance equal or higher to 5% for a 196 bp cyt b fragment ( fig. 1; supplementary table S3). This strategy highlights the distinct levels of divergence for the currently recognized species, allowing the formulation of recommendations for specific status of already described species (table 1). Furthermore, our barcoding approach, using a short cyt b fragment that fully recovers all phylogenetic units (supplementary fig. S2), may be a useful tool to assess the distinctiveness of newly discovered populations that are suspected to represent new taxa (Hebert and Gregory, 2005).
Ultimately, the taxonomic status of candidate species should be best addressed in an integrative fashion, that is, by searching for concordant differences in genetic, morphological and ecological traits (Padial et al., 2010). Since mtDNA could be particularly misleading in vipers, both by inflating the number of taxa or by missing lineages that have lost their mtDNA due to introgression (e.g., Zinenko et al., 2016), we encourage the use of multilocus genetic data (e.g., UCE loci, Blair et al., 2019). Until now, among all Eurasian vipers, only three species have been described using integrative approaches (i.e., by addressing phylogenetic divergence and characterizing phenotypic variability; Vipera walser Ghielmi et al., 2016; Macrovipera razii Oraie et al., 2018; Vipera sakoi Tuniyev et al., 2018). In the absence of more integrative studies, which are sometimes constrained by low financial support and by the scarcity or remoteness of populations (e.g., Göçmen et al., 2014aGöçmen et al., , b, 2017Freitas et al., 2018), the systematic situation of the group will likely remain unresolved in the immediate future.
Species classification is not the sole controversial issue in Eurasian viper taxonomy; allocation below (i.e., subspecies) or above (i.e., genera) species level is even more problematic due to the lack of objective, operational concepts defining these ranks. While our phylogenetic inferences are directed towards assessing the validity of species as independently evolving lineages, some provisional consensus can be reached at higher hierarchical levels. For instance, if Pelias is used to refer to Euro-Siberian Vipera species (as proposed in several studies, see Tuniyev et al., 2009Tuniyev et al., , 2012Tuniyev et al., , 2013Tuniyev et al., , 2018aAvcı et al., 2010), the designation of an additional genus within the current Vipera and the separation of Daboia into two or three genera would be required to reflect equivalent phylogenetic distances (see fig. 1; supplementary table S3). However, in order to avoid illfounded splitting procedures and even further confusion, we advise against taking such steps (see also Vences et al., 2013), and argue that the use of subgenera may be a better way of providing names for clades without disrupting the binomial nomenclature (Wallach, Wüster and Broadley, 2009).

Concluding remarks
In this work, we integrate currently available information on the phylogenetic and morphological variability of Eurasian vipers to advance into a more coherent and objective taxonomy for this group. Based on our integrative assessment, we provide recommendations on the specific status of 40 described species and propose some guidelines to clarify the taxonomic status of some of them (table 1). Species complexes such as V. latastei-monticola or V. ammodytestranscaucasiana require further analyses on the extent of gene flow among distinct lineages to delineate taxonomic units, while some of the currently recognized species, described using morphological data only (e.g., V. altaica, V. magnifica, V. shemakhensis), must be explicitly regarded as non-valid species due to low genetic differentiation in relation to other previously recognized species. Again, other taxa, such as eastern M. xanthina and V. sakoi require more extensive geographic sampling to arrive at robust conclusions. As discussed above, integrative taxonomic approaches bringing together independent evidence and using different methodological approaches, particularly incorporating genomic data instead of relying solely on mitochondrial data, will allow the robust delineation of coherent evolutionary units in a unified taxonomic framework.
Previous attempts to propose priorities for the conservation of vipers (i.e., Maritz et al., 2016) suffered from taxonomic inflation and lack of a geographic comprehensive sampling scheme (see table 1). It is striking that more than half of the Eurasian viper species listed as globally endangered in our assessment (nine of 16 species with categories CR, EN and VU; table 1) were classed as of doubtful validity as species. Advancing in a robust, evidence-based designation of taxonomic units is therefore essential for the future development of conservation strategies aimed to anticipate threats related to anthropogenic factors, while slowing down or even stopping the rapid decline of many populations of Eurasian vipers.