The biogeography of Elaphe sauromates (Pallas, 1814), with a description of a new rat snake species

Background The rat snake genus Elaphe once comprised several dozens of species distributed in temperate through tropical zones of the New and Old World. Based on molecular-genetic analyses in early 2000s, the genus was split into several separate genera, leaving only 15 Palearctic and Oriental species as its members. One of the three species also occurring in Europe is Elaphe sauromates, a robust snake from the Balkans, Anatolia, Caucasus, Ponto-Caspian steppes, and Levant that has been suspected to be composed of two or more genetically diverse populations. Here, we studied the genetic structure and morphological variation of E. sauromates, aiming to better understand its inter-population relationships and biogeography, and subsequently revise its taxonomy. Methods We reconstructed the phylogeography and analyzed the genetic structure of E. sauromates populations originating from most of its geographic range using both mitochondrial (COI, ND4) and nuclear (C-MOS, MC1R, PRLR, RAG1) DNA gene fragments. We employed Maximum likelihood and Bayesian inference methods for the phylogenetic tree reconstructions, supplemented with species delimitation methods, analysis of haplotype networks, and calculation of uncorrected p-distances. Morphological variation in 15 metric and 18 meristic characters was studied using parametric univariate tests as well as multivariate general linearized models. In total, we analyzed sequences originating from 63 specimens and morphological data from 95 specimens of E. sauromates sensu lato. Results The molecular phylogeny identified two clearly divergent sister lineages within E. sauromates, with both forming a lineage sister to E. quatuorlineata. The genetic distance between them (5.80–8.24% in mtDNA) is similar to the distances among several other species of the genus Elaphe. Both lineages are also moderately morphologically differentiated and, while none of the characters are exclusively diagnostic, their combination can be used for confident lineage identification. Here, following the criteria of genetic and evolutionary species concepts, we describe the lineage from eastern Anatolia and parts of the Lesser and Great Caucasus as a new species E. urartica sp. nov. Discussion Elaphe urartica sp. nov. represents a cryptic species whose ancestors presumably diverged from their common ancestor with E. sauromates around the Miocene-Pliocene boundary. The intraspecific genetic structure indicates that the recent diversity of both species has been predominantly shaped by Pleistocene climatic oscillations, with glacial refugia mainly located in the Balkans, Crimea, and/or Anatolia in E. sauromates and Anatolia and/or the Caucasus in E. urartica sp. nov.


INTRODUCTION
The Western Palearctic region is home to great diversity of herpetofauna (Sindaco & Jeremčenko, 2008;Sindaco, Venchi & Grieco, 2013). Among the most emblematic members are the rat snakes of the genus Elaphe Fitzinger in Wagler, 1833, once a large genus with distribution in temperate, subtropical, and tropical zones of both eastern and western hemisphere (Schulz, 1996). With the rise of molecular taxonomy in the last decades of the 20th century, the complicated relationships among dozens of species have been studied with a new rigor, resulting in the New-World as well as most Old-World species being identified as new or resurrected genera Lenk, Joger & Wink, 2001;Chen et al., 2017). As a consequence, the genus Elaphe is comprised of 15 species (Wallach, Williams & Boundy, 2014;Chen et al., 2017), with most occurring exclusively in Asia. The high diversity in Asia supports the hypothesis of the origin of the genus in this region (Lenk, Joger & Wink, 2001;Chen et al., 2017). Outside of Asia, three Elaphe species are also found in southern-eastern and eastern Europe: Elaphe sauromates (Pallas, 1814), E. quatuorlineata (Bonnaterre, 1790), the species with the western-most distribution of all Elaphe, and E. dione (Pallas, 1773) with mainly Asian distribution (Sindaco, Venchi & Grieco, 2013). Elaphe sauromates is also the type species of the genus Elaphe.
Elaphe sauromates is a robustly built snake whose biology has been systematically understudied and, despite its relatively large size, is only rarely encountered in the wild (Böhme & Ščerbak, 1993). It was described as a separate species under the name Elaphis sauromates (see : Strauch, 1873), followed by being synonymized with Coluber quatuorlineatus Bonnaterre, 1790, and subsequently treated as its subspecies C. quatuorlineatus sauromates (Boulenger, 1894). At the beginning of the 21st century, morphological and molecular-phylogenetic studies justified the species status of E. sauromates Lenk, Joger & Wink, 2001; and this view has been largely accepted to date (Sindaco, Venchi & Grieco, 2013;Sillero et al., 2014;Wallach, Williams & Boundy, 2014;Jablonski, Soltys & Simonov, 2019). Aside from confirming the separate species status of E. sauromates, early genetic studies found a surprisingly high genetic distance between populations from Kazakhstan and eastern Turkey (5.0-6.5% in Lenk, Joger & Wink, 2001 or even 7.6% in Huang et al., 2012), which was quite similar to the distances between E. sauromates and E. quatuorlineata, E. dione and E. bimaculata Schmidt, 1925, E. quadrivirgata (Boie, 1826 and E. carinata (Günther, 1864), or between E. quadrivirgata and E. schrenckii (Strauch, 1873) Huang et al., 2012). In addition to genetic divergence, some morphological and color variation seems to be geographically correlated, indicating possible phenotypic differentiation of at least two genetic forms (Zinner, 1972;Tiedemann & Häupl, 1978;Schulz, 1996). Thus, there is a clear indication that the taxonomy of E. sauromates requires a revision based on a detailed study of both genetic and morphological variation.
Here, we studied the genetic structure and morphology of E. sauromates in biogeographic framework using robust geographic sampling, with the aim of elucidating the levels of divergence and phylogenetic relationships among the populations originating from most areas of the known taxon range. We analyzed sequences of mitochondrial and nuclear DNA fragments, snake measurements and pholidosis and have revealed the existence of cryptic diversity within E. sauromates. As a result, we describe a new species of a rat snake of the genus Elaphe.

MATERIAL AND METHODS
Tissue sampling, DNA extraction, and PCR amplification Between 2000 and 2017 we collected tissue samples of E. sauromates sensu lato (s.l.) from an area extending from the eastern Balkans to western Kazakhstan. We covered almost the entire distribution range except for the central-southern Turkey, northern Iran, Turkmenistan, Uzbekistan, and the Levant (Israel, Lebanon, Syria; Fig. 1). In total we obtained 90 tissue samples of E. sauromates s. l., and 63 of them were subsequently used in our molecular analyses. All sample information (DNA working codes, sampling localities, GenBank Accession Numbers, and comments) is presented in Table 1. We collected blood samples, ventral scales, and saliva from living specimens, while muscle and liver samples were collected from the dead (road-kills) or museum deposited individuals. All tissue samples were preserved in 96% ethanol or frozen and stored at -25 or -80 C.
Total genomic DNA was extracted using the NucleoSpin Tissue kit (MACHEREY-NAGEL, Duren, Germany), following the manufacturer's instructions. Sequences of two mitochondrial (mtDNA) and four nuclear DNA (nDNA) genes were targeted: the DNA barcoding segment of cytochrome c oxidase subunit 1 (COI), the mitochondrial protein-coding segment of NADH dehydrogenase subunit 4 (ND4) (including the flanking tRNAs Serine, Histidine, and part of Leucine), and the nuclear protein-coding genes for the prolactin receptor (PRLR), the oocyte maturation factor Mos (C-MOS), recombination activation gene 1 (RAG1), and melanocortin 1 receptor (MC1R). Primers used for PCR amplifications are listed in Table S1. All amplification procedures involved an initial cycle of denaturation at 94 C for 3-7 min, 35-40 subsequent cycles at 94,  with the ABI PRISM 3100 capillary sequencer using the PCR amplification primers. All newly obtained sequences have been deposited in GenBank (Table 1).
Alignment, genetic divergence and model selection DNA sequences were manually checked, aligned, and inspected using BioEdit 7.0.9.0 (Hall, 1999). No stop codons were detected with into amino acid translated sequences using the vertebrate genetic code in the program DnaSP 5.10 (Librado & Rozas, 2009). The same program was used to calculate uncorrected p-distances among the main clades and to estimate the haplotype diversity (Hd), number of segregating sites (S), variables (V), nucleotide diversity (π), and parsimony informative (Pi) sites for the selected groups. Heterozygous positions in nuclear genes were manually identified based on the presence of double peaks in chromatograms. Identified heterozygous loci were coded according to the IUPAC ambiguity codes. The best-fit codon-partitioning schemes and the best-fit substitution models for phylogenetic analyses with concatenated dataset (COI 581 bp + ND4 810 bp) were selected using PartitionFinder v1.1.1. (Lanfear et al., 2012) with the following parameters: Bayesian approach (BA)-linked branch length; all models; BIC model selection; greedy schemes search; data blocks by codons for each used marker. The best partitioning scheme and models of nucleotide substitutions were: first and second positions (HKY+G), third position (K80). The Maximum likelihood (ML) analysis followed the same approach; the best substitution model in this case was GTR+G+I with a single partition.

Phylogenetic and haplotype network analysis
Concatenated (COI + ND4) mitochondrial phylogenetic trees were inferred using the BA performed with MrBayes 3.2.1 (Ronquist et al., 2012) and ML analysis performed with RAxML 8.0 (Stamatakis, 2014). The BA analysis was set as follows; two separate runs, with four chains for each run, 10 million generations with trees sampled every 100th generation. The first 20% of trees were discarded as the burn-in after inspection for stationarity of log-likelihood scores of sampled trees in Tracer 1.6 (Rambaut et al., 2013; all parameters had effective sample sizes (ESSs) of >200). A majority-rule consensus tree was drawn from the post-burn-in samples and posterior probabilities were calculated as the frequency of samples recovering any particular clade. Nodes with posterior probability values ! 0.95 were considered as strongly supported. The ML clade support was assessed by 1,000 bootstrap pseudoreplicates. Sequences of other, closely related Eurasian rat snakes, Oocatochus rufodorsatus (Cantor, 1842) and Zamenis lineatus (Camerano, 1891) (for GenBank accession numbers see Table 1), were used in the analyses included as outgroups Chen et al., 2017).
In order to obtain better support for putative species boundaries of the divergent lineages of Elaphe, we employed three species delimitation methods: (i) Bayesian implementation of the Poisson tree processes model (bPTP; Zhang et al., 2013, https://species.h-its.org/), (ii) multi-rate Poisson Tree Processes (mPTP; Kapli et al., 2016) model, using the webserver (http://mptp.h-its.org/), and (iii) the general mixed yule-coalescent model (GMYC; Pons et al., 2006). For species delimitations analysis we used as input an ultrametric tree of mtDNA haplotypes constructed with BEAST 2.1 (Bouckaert et al., 2014). BEAST analyses were run under the uncorrelated log-normal relaxed clock approach with a Yule tree prior. Two independent runs were conducted with a chain length of 5 Â 10 7 iterations. Tracer 1.6 (Rambaut et al., 2013) was used to check for convergence of the chains and adequate ESSs. Independent runs were combined using LogCombiner (part of the BEAST package), discarding the first 25% of each run as burn-in. The maximum clade credibility tree was summarized with TreeAnnotator (part of the BEAST package) and visualized with Fig-Tree 1.4 (http://tree.bio.ed.ac.uk/software/ figtree). GMYC species delimitation was conducted using the "splits" package (Ezard, Fujisawa & Barraclough, 2009) in R (R Development Core Team, 2019) under the single-threshold method.
The genealogical relationships between haplotypes of mtDNA and each nDNA marker were separately assessed with haplotype networks. For the purpose of allele network construction, sequences with more than one heterozygous site were resolved in PHASE 2.1.1 (Stephens, Smith & Donnelly, 2001) for which the input data were prepared in SeqPHASE (Flot, 2010). PHASE was run under default settings except for the probability threshold, which was set to 0.7. Haplotype networks of all analyzed markers were examined and drawn using PopArt (http://popart.otago.ac.nz) and the implemented parsimony network algorithm of TCS (Clement, Posada & Crandall, 2000), with 95% connection limit. Independent networks are considered distinct evolutionarily significant units, following Fraser & Bernatchez (2001), thus this analysis was also used to infer genetic structure within the studied taxa.

Material for morphological analyses and collection of specimens
In total, we examined external morphology of 95 specimens (46 males, 29 females, 20 individuals of unidentified sex) of E. sauromates s. l. from the Crimean Peninsula, Turkey, Armenia, and Azerbaijan ( Fig. 1D; Table 2). Only snakes with SVL ! 650 mm were used in the analyses of metric characters (75 adult animals in total), while all available specimens with known sex were used in scale count descriptive statistics and comparisons. The material was obtained either directly in the field or from collections of the following institutions: Zoology Department of Ege University, Turkey (ZDEU); Institute of Zoology, Azerbaijan, National Academy of Sciences, Baku, Republic of Azerbaijan (IZANAS); Institute of Zoology, Russian Academy of Sciences, Saint Petersburg, Russia (ZISP); Zoological Department of Tula State Regional Exotarium, Ministry of Culture of Tula Region, Tula, Russia (TuRE).
Field permits for this study were issued by the Republic of Turkey Ministry of Forest and Water Affairs (No. B.23.0.DMP.0.15.01-510-38417) and by the Crimean authority under the title "Study of Biodiversity and Landscape Structure of the Southeastern  Crimea, monitoring of Biotic and Abiotic Components of the Regional Ecosystems" (FASO: AAA-A16-116022510087-5) and by "Study of the structure and dynamics of land ecosystems in various climatic zones" No. AAAA-A19-119012490044-3). All efforts were made to minimize animal suffering.

Analysis of the external morphology
Body and tail lengths were measured with a ruler to the nearest one mm, while the head and individual skin plates were measured with a caliper with accuracy to the nearest 0.01 mm. For the purpose of obtaining basic information, the head scales of the juvenile snakes from Armenia were measured using their photographs (see Kropachev et al., 2015), but they were not used in further statistical analyses. In total, we targeted 15 metric and 18 meristic (scale count) characters, however not all measurements or counts were available for each individual (see Sahlean et al., 2016; Table S2 for explanation of the characters). In addition to the measurements and scale counts, we calculated two main morphometric indices: the snout-vent length to tail length (SVL/TL) and the snout-vent length to head length (SVL/HL) ratio.
Prior to morphological analyses, the snakes were split into two groups based on the results of molecular-genetic analyses and their origin: E. sauromates sensu stricto (s. s.; 63 specimens) and a new species that is formally described and named below (32 specimens). All measures were a priori log-transformed and the homoscedasticity (equal variances) was checked by Levene's test before the comparisons were carried out. We used independent sample t-tests to compare the length (SVL and total) of both snake groups. The tail and head lengths were compared using multivariate analysis of variance with SVL as a co-variate (MANCOVA). The overall head size was further compared with MANCOVA with head length used as a co-variate. The scale counts, which are not dependent on the snake length, were compared with multivariate analysis of variance (MANOVA). In all tests, the sexes were treated separately. We also examined the sexual differences within both groups and employed the same statistical strategy as in the species comparisons. The significance level was set to 0.05 and Bonferroni correction was applied in multiple-comparison analyses. All analyses were performed in SPSS Statistics 17.0.0. (SPSS Inc., Chicago, IL, USA).

The type material note
While the holotype of the newly described species (see below) as well as eight paratypes are deposited in public museums or scientific institutions (see Tables 1 and 2), four paratypes are currently still alive at TuRE (as of March 2019; Tula, Russia). These individuals have been described in detail, photographed, and assigned identification numbers. Upon death, the animals will be fixed and deposited in one of the central zoological museums in Russia (ZISP and/or ZM MSU). While the use of live animals as type material is not recommended by some authors (Dubois & Nemésio, 2007), such practice is not prohibited by International Commission on Zoological Nomenclature (ICZN). We decided to include these specimens (originating from Armenia) among the type series to provide a better representation of the material that was predominantly used as tissue sample source for molecular-genetic analyses.

Nomenclatural note
The electronic version of this article in portable document format will represent a published work according to the ICZN, and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank Life Science Identifiers (LSIDs) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix http://zoobank.org/. The LSID for this publication is: urn:lsid:zoobank.org:pub:F0BF1D63-7BD5-4340-85A0-E044CC56CD5B. The online version of this work is archived and available from the following digital repositories: PeerJ, PubMed Central and CLOCKSS.

Phylogeography and genetic diversity
The final nucleotide dataset of 4,068 bp comprised sequences of two mtDNA genes: COI (581 bp), ND4 (810 bp), and fragments of four nuclear genes: C-MOS (459 bp), MC1R (650 bp), PRLR (552 bp), and RAG1 (1,016 bp). The phylogenetic trees (ML and BA) of mtDNA concatenated dataset showed the same topologies with very similar, high statistical supports for each node ( Fig. 2; Figs. S1 and S2). Therefore, we only present the BA tree here to show the topology and interrelationships among the studied snake lineages. We identified two major, deeply divergent clades within E. sauromates s. l. Together they form a sister clade to that of E. quatuorlineata sequences. Uncorrected p-distances (mtDNA) among the clades ranged between 5.80% (ND4) and 8.24% (COI; Table 3). One of these clades corresponds to E. sauromates s. s., while the other one represents a distinct phylogenetic lineage, which we consider a new, separate species. The clade corresponding to E. sauromates has been detected (based on mtDNA) from western and central Anatolia (Turkey), the Balkans (Bulgaria, Turkey), southern Ukraine, Crimean Peninsula, southern Russia, and western Kazakhstan. Samples from the new species originated from eastern Anatolia (Turkey), Armenia, Azerbaijan, Georgia, and Dagestan (southern Russia). Both species seem to be geographically separated, presumably by the Anatolian Diagonal, and partially, also by the Great Caucasus (Fig. 1).
All three methods of species delimitation, that is, bPTP, mPTP, and GMYC recognized the three divergent lineages as distinct entities among the analyzed members of Elaphe ( Fig. 1): E. quatuorlineata, E. sauromates, and the new species described below that has a sister relationship to E. sauromates (Fig. 2).

Comparison of Elaphe sauromates s. s. and the new species
For all morphological analyses, we split the snakes into two taxonomic groups corresponding to E. sauromates s. s. (abbreviated as ES in this and the following section) and the new species (EU; see Methods; Tables 1 and 2), and males and females were treated separately for the purpose of both descriptive statistics (Tables 4 and 5; Tables S3  and S4) and comparisons (Table S5).

Note:
Plates on the left side were measured.
Morphological data originating from the type material of the new species can be found in Table 6 and data on a few juvenile specimens originating from Armenia in Table S4. Due to very low numbers, these were not subjected to statistical comparisons.
There are also no sexual differences in scalation in ES (F(1,15) = 2.575, p = 0.237, Wilks' L = 0.082, n(M) = 13, n(F) = 4), but according to the between-subject comparisons some differences may exist in numbers of ventral and subcaudal scales. However, model showed differences in scalation between males and females of EU (F(1,20) = 3.995, p = 0.036, Wilks' L = 0.111, n(M) = 12, n(F) = 10), which was mainly driven by differences in numbers of subcaudal scale pairs (70 ± 3 in males, 64 ± 4 in females) and preocular (1-3 in males, 1-2 in females) scales. More thorough morphological comparison of larger and more complete datasets will be necessary to obtain a more robust picture of both interspecific and intersexual differences in both species.

Systematic account
Our findings indicate that the taxon E. sauromates s. l. (type locality of E. sauromates: Pre-Sivash area of the Crimean peninsula, the Perekop Isthmus and adjacent territories of the Lower Dnieper region; related sample nos. 1,178 and 1,179 in this study; see Table 1) is composed of two clearly genetically differentiated allopatric populations. The populations from Transcaucasia and eastern Anatolia, which are also morphologically differentiated from E. sauromates s. s., represent a cryptic phylogenetic lineage. In accordance with the definition of the genetic species concept (i.e., genetic species is a group of genetically compatible interbreeding natural populations that is genetically isolated from other such groups; Baker & Bradley, 2006) and evolutionary species concept (i.e., populations with a long independent evolutionary history, represents a lineage of ancestral descendent populations, and maintains its identity from other lineages both on genetic (mtDNA and nDNA) and morphological level; Simpson, 1951;Wiley, 1978), we describe this lineage as a new species: Holotype. (Fig. 4) Tuniyev et al., 2014).
Full-size  DOI: 10.7717/peerj.6944/ fig-6 Abovyan, Armenia,40.37 N,44.69 E (Fig. 6A) Etymology. The specific epithet is a feminine adjective derived from the name of the ancient kingdom of Urartu that flourished in the Armenian Highlands and around lake Van, an area of recent distribution of E. urartica sp. nov., in the 9th-6th century BCE (Asher & Asher, 2009). We are choosing this name out of respect for Peter Simon Pallas, who proposed the name for E. sauromates, now the sister species of E. urartica, which most likely refers to Sarmatians (Sauromatae; Sayρομasaι in Greek), a confederation of nomadic peoples inhabiting vast portions of the recent range of E. sauromates between the 5th century BCE and 4th century CE.

Diagnosis.
A new species of western Palearctic genus Elaphe, very similar to E. sauromates (Pallas, 1814), characterized by the combination of the following characters: total length usually does not exceed 1,200 mm (796-1,205 mm), snout-vent (SVL) length usually less than 1,000 mm (650-970 mm), tail length less than 250 mm (146-245 mm) (see Tables 4 and 6). Tail forms about 25% of the SVL in males and about 21% in females. Head relatively large, distinguished from the body. Snout in prefrontal and internasal area is conspicuously convex which usually forms a hook-nosed head profile. Pileus length on average 1.8-1.9 times larger than its width. Frontal plate 1.2-1.3 times longer than wide. Anterior inframaxillar scute relatively large and wide, 1.2-1.3 times longer than the narrow posterior inframaxillar scute. One or two preocular scales, one loreal, two postoculars, two temporals, three or four posttemporals, eight labials, 10-11 sublabials on each side of the head. Eye in contact with fourth and fifth labials (Table 5; Table S3). Variation in head scale counts is relatively low (see Table S3). Usually two gulars located the anterior inframaxillars. The total number of gulars between inframaxillars and first preventral scale exceeds 12. Number of ventrals is 154-211 (154-206 in males, 194-211  Pileus is dark, often almost black, slightly lighter on the tip of the snout. Upper preoculars and temporals are dark forming a postocular stripe extending toward the mouth corner. This stripe blends with the dark dorsolateral head coloration anterior to the eye. Pale spots on the labials, only barely visible or lacking on sublabials. Ventral side of the body is whitish to pale yellow, sometimes with pinkish tint. There are marbled patterns of numerous small irregular dark brown and light gray spots with reddish contours that are more pronounced on the lateral sides of ventral plates. Throat is light, with numerous reddish-orange or brownish speckles on the lower jaws and anterior ventral plates. Iris is dark brown or almost black with a thin light rim around the pupil. Differential diagnosis. Elaphe urartica sp. nov. is closely related to E. sauromates and E. quatuorlineata. The genetic distance between E. urartica sp. nov. and E. sauromates is 7.20% and 6.91% in COI and ND4, respectively, and 0.15%, 0.15%, 0.73%, and 0.24%, in C-MOS, MC1R, PRLR, and RAG1, respectively. The genetic distance between E. urartica sp. nov. and E. quatuorlineata is 8.24% and 5.80% in COI and ND4, respectively, and 0.15%, 0.15%, and 0.20%, in C-MOS, MC1R, and RAG1, respectively. E. urartica sp. nov. is also morphologically very similar to E. sauromates. E. urartica sp. nov. attains shorter lengths than E. sauromates (795 ± 80 mm vs. 937 ± 152 mm in E. sauromates males, 861 ± 97 mm vs. 929 ± 160 mm in E. sauromates females; though E. sauromates from other parts of the range, for example, the Balkans, can be even larger than specimens in our dataset, see Sahlean et al., 2016). Both taxa also differ in some relative head dimensions. E. urartica sp. nov males have relatively (in comparison to the head length) longer pileus, higher rostrum, but shorter frontal plate and anterior inframaxillary scute. Females only differ in the pileus length (Table 4). The upper surface of the head is more convex near orbits, prefrontals, and internasals and the rostrum is more anteriorly pronounced in E. urartica sp. nov. than in E. sauromates. Males of both species also slightly differ in their scalation: E. urartica sp. nov. males have fewer subcaudal pairs (64 ± 4 vs. 75 ± 3 in E. sauromates) and loreal scales (1-2 vs. 1-3 in E. sauromates). Other differences in metric and meristic characters were not found statistically significant. The coloration of E. urartica sp. nov. is generally darker than that of E. sauromates (Figs. 4 and 6-10). The dorsal side of the head is very dark, sometimes almost black and without the whitish area separating two blotches just posterior to the head as seen in E. sauromates (Figs. 4 and 6-10). On the lateral side of the head the dark stripe running from behind the eye toward the corner of the mouth is also less distinguished in E. urartica sp. nov. compared to E. sauromates, in which it is clearly separated by lighter color from the darker head coloration. E. urartica sp. nov. has more conspicuous dorsal body spots that are more rounded than in E. sauromates, in which transverse elongation of the spots is common. These dorsal spots are typically lined with whitish color in E. urartica sp. nov. rather than yellow or yellowish as in E. sauromates. Description of the holotype. Adult male (ZDEU 26/2012; Fig. 4). Body cylindrical, snout-vent length 803 mm, tail length 210.00 mm. Head big and clearly distinct from the neck, pileus length 24.90 mm, width 16.20 mm. Head and body scales keeled. Rostral slightly curved toward the top of the head, indistinctly wedged between the internasals. Rostrum height 5.3 mm and width 6.6 mm, in contact with two labials, two nasals, and two internasals. Nostrils located within the nasal scales, inter-nostril width is 8.5 mm. Loreal on either side of the head in contact with second and third labials. Two preocular and postocular plates on each side of the head. Eyes circular with circular pupil of 4.8 mm diameter. Length of the narrow frontal 8.5 mm, width 6.6 mm. Eight labials, with fourth and fifth in contact with the eyes on each side. A total of 11 sublabials on each side, five sublabials in contact with the anterior chin shields on each side. Two temporals and two posttemporals on each side. Nine dorsal and temporal scales surrounding the posterior margin of the parietals. Two gular scales in contact with the anterior chin shield. A total of 23 dorsal scale rows at midbody, 25 at the level of one head length posterior to the head, and 19 at one head length anterior to the cloaca level. One preventral, 202 ventral plates, two anal plates, and 72 and 73 subcaudals on each side, respectively (Table 4). The dorsal side of the body is yellowish with round black spots. Flanks with two rows of dark blotches. Dorsal and lateral spots form striped pattern on the tail. The top of the head is black. Temporal stripe distinct, blending with the dark dorsolateral head coloration anterior to the eye. The ventral side of the body is yellowish white, with black markings (Fig. 4).
Variation. Details on variation among the type specimens of E. urartica sp. nov. are presented in Table 6. The coloration of the paratypes is very similar to that of the holotype.
Distribution and habitat. The geographic range of E. urartica sp. nov. is bordered by the Armenian Plateau, south-eastern foothills of the Great Caucasus, Alazan Valley, Kur-Aras, Lenkoran Lowlands, and the area of Qobustan. The species is distributed in Turkey, Georgia, Armenia, Azerbaijan, Nagorno-Karabakh, Iran, and Russia. In Turkey, it can be found east of the Anatolian Diagonal with reliable records from Kars, Bitlis, Diyarbakır, and Van Provinces, presumably also in Erzurum, Iğdır, and Ağrı Provinces (Baran et al., 2012). In eastern Transcaucasia E. urartica sp. nov. is distributed from south-eastern Georgia to the Zalka Plateau or to Suramskyi Ridge in Southern Ossetia in the West, throughout most of the Armenian territory, Nagorno-Karabakh, and Azerbaijan with the exception of the Abşeron Peninsula. The eastern part of the range lies in northern Iran to the Golestan Province to the East, and Kermanshah and Semnan Provinces to the South (Alekperov & Loginov, 1953;Muskhelishvili, 1970;Flärdh, 1983;Schulz, 1996;Sindaco et al., 2000;Arakelyan et al., 2011;Bunyatova, Akhmedov & Dzhafarov, 2012;Bunyatova, 2013;Najafov, Hashimov & Isgenderov, 2013;Safaei-Mahroo et al., 2015).
In the Russian Federation, E. urartica sp. nov. occurs in Samur-Devichi Lowlands of southern Dagestan and probably in the Dagestan Intermontane Region as well (Ananjeva et al., 2006;Mazanaeva & Askenderov, 2014). The species could also occur in the extreme northern regions of Iraq (Sindaco, Venchi & Grieco, 2013). The snake occurs in a wide range of altitudes-from ca. 25 m below sea level in the Lenkoran foredeep to about 2,600 m a.s.l. in the Shirak Province in Armenia (Arakelyan et al., 2011). It is an eurytopic species inhabiting a wide variety of landscapes: mountain and lowland semi deserts, different types of the steppe, semi subtropical savannah-like forest-steppes with oreoxerophytes, sparse juniper forests, montane broad-leaved forests, and alpine meadows (Fig. 5). The climate within the E. urartica sp. nov. range varies from the subtropical in Lenkoran and piedmont area of eastern Transcaucasia to cold mountain climate in Armenia and north-eastern Anatolia. Humidity varies from highly arid (with annual precipitation of less than 200 mm) to moderately humid (1,400-1,600 mm per year; Clark & Clark, 1973;Arakelyan et al., 2011;Bunyatova, Akhmedov & Dzhafarov, 2012;Şensoy et al., 2016).
Elaphe urartica sp. nov. is sympatric with E. dione in Dagestan, central-eastern Azerbaijan, eastern Georgia, and presumably in north-eastern Turkey, southern Armenia, and northern Iran. All other species of the genus Elaphe have allopatric distribution relative to E. urartica sp. nov. Since the species occurs in a region of southern Russia (Dagestan), north of the Caucasus, that is geographically and politically considered a part of Europe (Sillero et al., 2014), E. urartica sp. nov. is considered another member of the European herpetofauna.
Conservation status. Elaphe sauromates is a species with declining population numbers and is listed (as E. quatuorlineata) in the Red Data Books of: Ukraine (1994)-category 3, Kazakhstan (1996)-category 4-a little-known species, and Turkmenistan (1999)category 3-a rare species on the periphery of the distribution range. It is also listed among the taxa recommended to be included in the new edition of the Red Data Book of the Russian Federation as vulnerable species (requiring priority of protection measures-III) (Ananjeva et al., 2006;Iljashenko et al., 2018). The current global IUCN Red List status for E. sauromates s.l. is of Least Concern (Aghasyan et al., 2017). Since E. urartica sp. nov. occurs in the parts of the territory from which the data for these lists was derived, similar conservation concerns might apply for it as well. However, the consequence of splitting E. sauromates s. l. into two separate species is that the ranges of both species have significantly decreased, and this should be considered in future conservation measures as well.
Proposal of common names. We propose the English name "Urartian Rat Snake" for E. urartica sp. nov. Along with the name "Blotched rat snake", we also suggest using the name "Sarmatian Rat Snake" for E. sauromates, instead of the older "Eastern Four-lined Rat Snake" derived as a subspecific name from the common name of E. quatuorlineata. The newly proposed name would decrease confusion and also better reflects the scientific name of E. sauromates.

Systematics of Elaphe sauromates and Elaphe urartica sp. nov
Since the publication of Boulenger's Catalogue of the snakes in the British Museum (1894) E. sauromates had been considered an eastern subspecies of E. quatuorlineata. Based solely on genetic distances between the two taxa,  and Lenk, Joger & Wink (2001) proposed their split into two sister species, which has been largely accepted. The distribution of both species is parapatric (Sindaco, Venchi & Grieco, 2013;Kornilios et al., 2014). While E. quatuorlineata and E. sauromates are clearly and easily distinguishable, mainly based on adult pattern and coloration (Böhme & Ščerbak, 1993;Schulz, 1996), E. urartica sp. nov. represents a cryptic species, only moderately phenotypically differentiated from E. sauromates. Based on the available evidence, both species currently have allopatric ranges, presumably separated by the Anatolian Diagonal and the Great Caucasus ( Fig. 1; see Schulz, 1996;Sindaco et al., 2000). Phylogenetically E. urartica sp. nov. and E. sauromates represent sister lineages, and together they form a lineage sister to E. quatuorlineata (Fig. 2). Relatively little is known about intraspecific variation of E. urartica sp. nov. but haplotype diversity indicates that it is much lower than in E. sauromates (Figs. 1 and 3; Table 3, Results). The species occurs in a relatively small range (Fig. 1), however, in the region known as a radiation and speciation center (Tuniyev, 1995). While it does not seem likely that E. urartica sp. nov. is further differentiated into the subspecies in eastern Anatolia or Transcaucasia, the Caucasian montane part of the range of E. sauromates s. l. has yet to be genetically analyzed, and thus it remains unclear what species inhabits this region. Using phenotypic identification (i.e., without genetic analysis), it seems that both species come into possible parapatry in Dagestan-with E. sauromates in northern plains and E. urartica sp. nov. in more southerly located montane habitats (B.S. Tuniyev, 2018, personal observation). E. sauromates has been found relatively close to the southern ranges of the Anatolian Diagonal (Kayseri Province, Turkey, sample 2,892, see Table 1; Fig. 1A). On the other hand, the Anatolian Diagonal has formed a natural barrier to this species' ability for dispersal southward from western and central Anatolia (Davis, 1971;Nilson, Andrén & Flärdh, 1990;Jandzik et al., 2018), so it is more likely that the Levantine populations actually belong to E. urartica sp. nov. Given the commonly observed biogeographic pattern in this region (Jandzik, Avcı & Gvoždík, 2013;Jandzik et al., 2018;Stümpel et al., 2016;Tamar et al., 2016;Kornilios, 2017), it is also possible that these, likely, isolated populations have diverged from one of these species (or their common ancestor) and represent a yet undiscovered taxon (Zinner, 1972).
Elaphe urartica sp. nov. represents a cryptic taxon within E. sauromates s. l. While we found some phenotypic differences, mainly being relatively smaller in size, having a moderate differences in the shape of the head, some head scalation dimensions, and rather different adult coloration, both species are still morphologically very similar (see Figs. 4 and 6-8 for E. urartica sp. nov. and 9, 10 for E. sauromates), especially when compared to their sister taxon E. quatuorlineata. However, this is also not unusual among closely related snake species, even within the genus Elaphe-E. dione, E. bimaculata, and E. zoigeensis Huang, Ding, Burbrink, Yang, Huang, Ling, Chen & Zhang, 2012 are virtually indistinguishable based on external morphology (Schulz, 1996;Huang et al., 2012).
This similarity could be the result of similar environmental conditions and, consequently, the absence of selective pressure toward other type of coloration, size, or proportions (Fišer, Robinson & Malard, 2018). Though none of the morphological characters seem to be species-specific, using their combination will allow for confident specific identification. More robust sampling for morphological analyses will shed additional light on the levels of differentiation between both species.

Historical biogeography
The fossil record indicates that the ancestors of modern E. sauromates s. l. were present in the region of the recent range since at least Late Pliocene. Snake remnants identified as either E. quatuorlineata or E. cf. quatuorlineata, presumably belonging to E. sauromates s. l., are known from the Early-Middle Pleistocene of the Balkans, Transcaucasia, and central Anatolia (Szyndlar, 1984(Szyndlar, , 1991Venczel & Şen, 1994;Venczel, 2000;Venczel & Várdai, 2000;Kornilios et al., 2014;Wallach, Williams & Boundy, 2014;Blain et al., 2014;Blain, 2016). In addition, Elaphe aff. quatuorlineata fossils were also described from the Middle Pliocene of Moldova (Redcozubov, 1991(Redcozubov, , 2005, indicating that parts of the range were probably inhabited even earlier. Molecular-phylogenetic analyses suggest that E. sauromates s. l. and E. quatuorlineata split from their last common ancestor sometime in the Late Miocene (7.3-8.3 Mya), presumably during the formation of the mid-Aegean Trench that separated the Balkans from Asia Minor (Lymberakis & Poulakakis, 2010;Kornilios et al., 2014). Further evolution of E. sauromates s. l. is mainly associated with Anatolia and adjacent areas of the southern Caucasus. The mtDNA genetic distances between E. sauromates s. s. and E. urartica sp. nov. (Table 3) place the split between these two lineages at the Miocene-Pliocene boundary (five to eight Mya; see Kornilios et al., 2014), so they presumably separated from each other not very long after their ancestors separated from E. quatuorlineata. Although our estimation is only approximate, it coincides with the splits among other snake and lizard taxa from the same region, most notably the montane vipers of Montivipera xanthina complex-M. xanthina clade from western Anatolia and M. bornmuelleri clade from the east and south of the Anatolian Diagonal separated around five to six Mya (Stümpel et al., 2016). However, these authors propose that the valley of Göksu River south-west of the Anatolian Diagonal formed the natural barrier for the vipers adapted to the life in high altitudes, which was probably not true for the eurytopic rat snakes. Nevertheless, similar patterns of the east-west Anatolian split have been described in numerous other reptilian taxa (see Fritz et al., 2008;Kyriazi et al., 2008;Kornilios et al., 2012;Ahmadzadeh et al., 2013;Jandzik, Avcı & Gvoždík, 2013;Jandzik et al., 2018;Kapli et al., 2013;Skourtanioti et al., 2016).
Elaphe sauromates s. s. then presumably survived the Pleistocene climatic oscillations in refugia located mainly in Anatolia, Crimea, and/or the southern Balkans, where the highest haplotype diversity was recorded (see Figs. 1B and 1C). Subsequently in the Late Pleistocene, E. sauromates s. s. reoccupied the territories of eastern Europe concomitant with dispersing north of the steppe and semi-desert habitats and the Black Sea regression (Ratnikov, 2014), similarly as proposed for the snake Dolichophis caspius (Nagy et al., 2010) and lizards Lacerta viridis (Marzahn et al., 2016), Pseudopus apodus (Jandzik et al., 2018), and Podarcis tauricus (Psonis et al., 2018). Before the terminal phase of the last Pleistocene glaciation, E. sauromates already lived in the Crimean Mountains as confirmed by the fossil record (Ratnikov, 2015), although remnants of E. quatuorlineata-sauromates from the Pleistocene and Holocene are not known from the East European Plain (Ratnikov, 2002). Interestingly, Mazanaeva & Tuniyev (2011) proposed that E. sauromates and other Mediterranean faunal elements might have colonized the Great Caucasus from the North through the Kuma-Manych Strait during the last phases of the Pleistocene. Our lack of samples from the south-eastern parts of E. sauromates s. s. range does not allow us to confirm or reject this scenario. Colonization of the Transcaspia and Central Asia is probably very recent as the region is inhabited by the common northern E. sauromates s. s. haplotypes (Figs. 1A-1C) and might have coincided with the Late Pleistocene or even the Early Holocene regression of the Caspian Sea (see Krijgsman et al., 2018).
The significantly lower genetic diversity of E. urartica sp. nov. (only two haplotypes in COI, three in ND4, one to three in the nDNA genes; Figs. 1 and 3) does not provide enough information for drawing viable biogeographic scenarios. The species probably survived the Pleistocene glaciations in one, or a few, Transcaucasian or east-Anatolian refugia (Kahle et al., 2011;Kornilios et al., 2011;Jandzik et al., 2018) and spread relatively fast to the north-eastern Caucasus and southern Anatolia. A very similar situation with low genetic variation in the same region was observed in an anguid lizard Pseudopus apodus (Jandzik et al., 2018). Additional samples from Iran, central Anatolia, and the Levant will be necessary to learn more about the history of this species. In particular, one might expect surprises from the Levantine region given that it hosts many endemic taxa and lineages, and several Anatolian species have their subspecies or sister species located in the region, for example, Natrix tessellata complex (Guicking, Joger & Wink, 2009), the L. trilineata complex (Ahmadzadeh et al., 2013), Z. hohenackeri (Jandzik, Avcı & Gvoždík, 2013), the Ablepharus kitaibelii complex (Skourtanioti et al., 2016), Xerotyphlops vermicularis (Kornilios, 2017), and the Mediodactylus kotschyi complex .

CONCLUSIONS
Here we studied the biogeography of the rat snake E. sauromates from the Balkans, Anatolia, Caucasus, and Ponto-Caspian region using both molecular and morphological data. We found that the taxon is, in fact, comprised of two distinct evolutionary lineages and the cryptic lineage represents a new species that we name E. urartica sp. nov. Both species split from their common ancestor around the Miocene-Pliocene boundary and their recent genetic structure was mainly influenced by Pleistocene climatic oscillations.
Çetin Ilgaz contributed reagents/materials/analysis tools, approved the final draft. Ekaterina Polyakova contributed reagents/materials/analysis tools, approved the final draft. Konstantin Shiryaev contributed reagents/materials/analysis tools, approved the final draft. Boris Tuniyev contributed reagents/materials/analysis tools, prepared figures and/or tables, approved the final draft. David Jandzik conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft.

Field Study Permissions
The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers): Field study permits were issued by the Republic of Turkey Ministry of Forest and Water Affairs (No. B.23.0.DMP.0.15.01-510-38417), by the Crimean authority under the title "Study of Biodiversity and Landscape Structure of the Southeastern Crimea, monitoring of Biotic and Abiotic Components of the Regional Ecosystems" (FASO: AAA-A16-116022510087-5) and by "Study of the structure and dynamics of land ecosystems in various climatic zones" No. AAAA-A19-119012490044-3).

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: The group of mitochondrial and nuclear DNA sequences described here is accessible via GenBank accession numbers MK640236 to MK640422.

Data Availability
The following information was supplied regarding data availability: All sequences are available in the GenBank database. Accession numbers for the DNA sequences included in this study are provided in Table 1.
The sequences are also available as a Supplemental File.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.6944#supplemental-information.