Morphological and genetic differentiation in the anguid lizard Pseudopus apodus supports the existence of an endemic subspecies in the Levant

The Levant represents one of the most important reptile diversity hotspots and centers of endemism in the Western Palearctic. The region harbored numerous taxa in glacial refugia during the Pleistocene climatic oscillations. Due to the hostile arid conditions in the warmer periods they were not always able to spread or come into contact with populations from more distant regions. One large and conspicuous member of the Levantine herpetofauna is the legless anguid lizard Pseudopus apodus. This species is distributed from the Balkans to Central Asia with a portion of its range running along the eastern Mediterranean coast. Mitochondrial and nuclear DNA sequences, microsatellite genotypes, and morphology show that populations in this region differ from the two named subspecies and presumably had a long independent evolutionary history during the Quaternary. Here we describe the Levantine population as a new subspecies and present biogeographic scenarios for its origin and diversification. The new subspecies is genetically highly diverse, and it forms a sister lineage to Pseudopus from the remaining parts of the range according to mtDNA. It is the largest-bodied of the three subspecies, but occupies the smallest range. Vertebrate Zoology 71, 2021, 175–200  |  DOI 10.3897/vz.71.e60800 Copyright Daniel Jablonski. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Daniel Jablonski et al.: Morphological and genetic differentiation in the anguid lizard Pseudopus apodus 176

A very similar pattern of distribution is also typical for the Levantine populations of Pseudopus apodus (Pallas, 1775) from the southernmost part of the species range (Sindaco and Jeremcenko 2008). This anguid species is the only living member of the genus Pseudopus (previously considered congeneric with Ophisaurus; e.g. , a large and robust legless lizard that is among the relatively common inhabitants of the Mediterranean maquis leaf litter. Some of us ) recently studied its genetic structure using the analyses of the mitochondrial and nuclear DNA sequences, and showed clear genetic divergence partly corresponding to the traditionally recognized, morphologically defined, subspecies: P. a. apodus distributed in Crimea, Caucasus, south Caspian Region, and Central Asia, and P. a. thracius (Obst, 1978) from northern and western Anatolia and the Balkans ) Obst 1978(. Beside these two subspecific lineages, the phylogeny revealed a lineage, occurring exclusively in the Levant , with unclear morphological affiliations ) (. This Southern lineage represents an early divergence within recent populations of P. apodus. Despite very limited sampling from the Levant (n = 3),  proposed that the biogeographic history of the Southern lineage differed from histories of P. a. apodus and P. a. thracius. These two subspecies have large ranges, are genetically less diverse, and their present distribution pattern, and genetic structure, indicate they only recently expanded to inhabit most of their ranges. The Southern lineage, on the other hand, has a relatively small range and appears to be genetically more diverse than the other two. This may reflect historical fragmentation and disruptions in gene flow between populations inhabiting a relatively small range surrounded by hostile arid environment unsui table for this anguid species ).
Here we studied genetic diversity and morphology of P. apodus from across its range taking an advantage of an improved geographic sampling, especially in the Levant, which had previously not been analyzed in depth. This allows us to better understand the patterns of genetic differentiation and obtain deeper insights into the biogeographic history of the species. As a result, we propose to adjust the intraspecific taxonomy by describing a third P. apodus subspecies.

Material and methods
We used 156 tissue samples of P. apodus for genetic analyses and we analyzed the morphology of 479 specimens covering the entire distribution range of the species (Table S1, S2 and Figure 1).

Genetics Tissue sampling and laboratory procedures
As a source of DNA, we used buccal swabs, blood (taken from the caudal vein) or scale clips from living individuals, and liver or muscle biopsies from ethanol-preserved museum specimens and road-killed animals. We studied variation in two mitochondrial DNA (mtDNA) markers and eight microsatellite loci. In total, we newly obtained DNA data from 52 individuals for mtDNA and 155 individuals for microsatellites and combined them with previously published data (Table S1).
Total genomic DNA was extracted from the tissue samples using the NucleoSpin Tissue kit (Macherey-Nagel, Düren, Germany), following the manufacturer's instructions. For molecular-genetic analyses, we used fragments of two mtDNA markers (NADH dehydrogenase subunit 2, ND2, and cytochrome b, Cyt b) and eight microsatellite loci. The PCR strategy for mtDNA is described in detail in  and was followed without modifications. The complete ND2 dataset of P. apodus contained 33 sequences (21 newly obtained, 12 published) and the Cyt b dataset contained 85 sequences (51 newly obtained, 34 published; Table S1).
Eight microsatellite loci (Geiser et al. 2013;Miku líček et al. 2018( were amplified in two multiplex PCRs (multiplex 1: AF38-VIC, AF44-FAM, AF47-PET, An-Fr21-NED; multiplex 2: AF19-FAM, AF37-NED, AnFr12-VIC, AnFr35-NED) using the Type-it Microsatellite PCR Kit (Qiagen, Hilden, Germany). PCRs were performed in a total volume of 10 μL containing 1× Type-it Multiplex PCR Master Mix, 0.2 μM of each primer pair )forward primers were fluorescently labeled with FAM, VIC, NED, and PET(, 2 μL of DNA, and RNase-free water. PCR amplification involved an initial cycle of denaturation at 95°C for 5 min and 30 subsequent cycles of 95°C for 30 s, 58°C )multiplex 1( or 55°C )multiplex 2( for 90 s, and 72°C for 30 s, followed by a final extension step at 60°C for 30 min. PCR products were run on an ABI 3130 genetic analyzer (Applied Biosystems, Foster City, CA, USA) with a LIZ-500 size standard. Microsatellite fragments were visualized using the software GeneMapper 3.7 (Applied Biosystems) and were scored manually by a single observer.

Analyses of mitochondrial DNA
After manual aligning and trimming the low-quality sequence ends the resulting alignments contained sequences of 727-bp fragment of ND2 and 805-bp fragment of Cyt b. These were used in the tree analysis, along with a concatenated ND2+Cyt b dataset (1532-bp) which was constructed by combining the two former alignments. For the network analysis of all available sequences, we used the same length fragment of ND2, although only 423-bp fragment of Cyt b. Only the 575-bp long fragment of Cyt b was used for network analysis from the Southern lineage. In both cases this was due to variable lengths of the available sequences. No stop codons were detected when checked in the program DnaSP 6.00 (Rozas et al. 2017).
We constructed Bayesian Inference (BI; MrBayes 3.2.6; Ronquist et al. 2012) and Maximum Likelihood (ML; RAxML 8.0.0; Stamatakis 2014) phylo genetic trees, for all three datasets, using all available sequences. The best-fit model of sequence evolution was selected using partitionfinder 2 (Lanfear et al. 2017) with the following parameters (single locus trees): 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 as follows: the first position K80+G, second position HKY+I, third position HKY for Cyt b and HKY+I, HKY, and HKY+G for ND2. The ML analysis followed the same approach; the best model in this case was GTR+G with two subsets in Cyt b and three subsets in ND2. The concatenated dataset was divided to six subsets based on four user schemes (unpart, gene-part, 3 rd -pos-extra, codon-part( with following final models for BI analysis: HKY+I (subset 1), HKY (2, 3, 5, 6), and JC+I (4), and GTR+G in all subsets for ML analysis. The BI analysis was set as follows: two separate runs, with four chains for each run, 10 million generations with trees sampled every 100 th 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. Genealogical relationships between mtDNA haplotypes were separately assessed with haplotype networks (32 sequences in ND2, 85 in Cyt b) based on alignments of the network analyses as mentioned above, and from the concatenated dataset. Haplotype networks of both analyzed markers were constructed and drawn using PopArt (http://popart.otago.ac.nz) and the implemented parsimony network algorithm of TSC (Clement et al. 2000), with 95% connection limit. Well defined networks are considered distinct evolutionarily significant units, following Fraser and Bernatchez (2001), thus this analysis was also used to infer genetic structure within the studied taxa.
DnaSP 6.00 (Rozas et al. 2017) was used to estimate the number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (π) for particular detected lineages of the species based on single locus datasets. The same program was used for uncorrected p distance calculations.
Principal Component Analyses (PCAs) of mtDNA markers were carried out using 33 sequences of both the concatenated and ND2 datasets, and 85 sequences of Cyt b dataset in the R package Adegenet (Jombart 2008) implemented in the R statistical environment (R Core Team 2020).

Analyses of microsatellites
We used two approaches to assess genetic structure among the samples: a Bayesian-based method implemented in the Structure 2.3.4 software (Pritchard 2000) and a Principal Component Analysis (PCA) implemented in the R package Adegenet (Jombart 2008). Using Structure, direct posterior probabilities for the number of clusters K [LnP)K(] as well as the ad hoc statistic ΔK (Evanno et al. 2005) were estimated assuming uniform prior values of K between one and ten. ΔK statistic was calculated using online program Structure Harvester 0.6.94 (Earl and vonHoldt 2012). In Structure, both admixture and non-admixture models, always in combination with uncorrelated allele model, were applied. The analyses were based on runs of 10 6 iterations, following a burn-in period of 10 4 iterations. A series of five independent runs for each K was made with the same parameters to test the accuracy of the results. PCA were run using the same data set as Structure.
To assess genetic diversity and differentiation between clusters, the samples were pooled into three geographical regions based on the outcomes of the mtDNA phylogenetic, Bayesian and PCA analyses (see below). A number of different )N A ) and private (N P ) alleles, allele frequencies, percentage of polymorphic loci (P), and observed (H O ) and unbiased expected (H E ) heterozygosity were calculated using GenAlEx 6.5 (Peakall and Smouse 2006). Genetic differentiation among the three clusters )subspecies( was estimated using both unbiased Nei's genetic distance D S (Nei 1972) and F ST statistics calculated via an Analysis of Molecular Variance (AMOVA) (Meirmans 2006). The significance of genetic differentiation )F ST ) among clusters was based on 9999 permutations. To understand how genetic variation is partitioned within and among clusters, we applied AMOVA based on F ST statistics. Analysis of genetic diversity, AMOVA, and the calculation of genetic distances were carried out using GenAlEx 6.5.

Studied material and morphological data
We examined the external morphology of 479 P. apodus specimens (Fig. 1C): 358 specimens of P. a. apodus (189 males, 72 females, 35 juveniles, and 62 adult specimens of unknown sex originated in Armenia, Azerbaijan, Crimea, Georgia, Iran, Kazakhstan, Russia, Tajikistan, eastern Turkey, Turkmenistan, and Uzbekistan), 17 specimens of P. a. thracius (11 males and 6 females; from Bulgaria and western Turkey), and 104 specimens from the Southern lineage (39 males, 19 females, and 46 juveniles; all from Israel(. The subspecific status of the examined specimens was determined based on distribution data, morphological characters, and molecular analyses. Differences between the total number of examined specimens and the number of analyzed specimens for each character are due to poor condition of some specimens (e.g. damaged parts of the head or body), which only allowed some characters to be recorded. The complete list of analyzed specimens and their geographic origin (if georeferenced) are available in Table S3 and Fig. 1C. The material was obtained either directly in the field or from the collections of the following institutions: The Steinhardt Museum of Natural History, Tel Aviv University, Tel Aviv, Israel (TAU); Institute of Zoology of the Republic of Kazakhstan, Almaty, Kazakhstan (IZRK); Zoological Museum in Kiev, Ukraine (ZMK); Zoological Museum of the Zoological Institute of the Russian Academy of Sciences (ZISP); Museum of Nature of V. N. Karazin Kharkiv National University, Kharkiv, Ukraine (SMNK); Museum of Zoology of the Senckenberg Natural History Collection Dresden, Dresden, Germany (MTD). Manipulation of live individuals was limited to minimize animal discomfort. No specimens from Israel were collected specifically for this study. All originated from the TAU collections.
Coloration of the body was recorded via direct observation in the field or from photographs of live specimens (when available). Coloration in preservative is based on observation of the specimens the abovementioned collections.

Analysis of external morphology
For the initial analysis of the differences in metric and meristic characters among the three groups (P. a. apodus, P. a. thracius and the Southern lineage, below described as a new subspecies) we used analysis of co-variance (AN-COVA( with subspecies and sex as main effects and SVL as a covariate. Analyses were conducted in R statistical environment (R CoreTeam 2020). Next, we performed three principal component analyses (PCAs): 1) with all adult specimens )SVL ≥ 310 mm, sensu  and all characters, 2) with only adult males, and 3) with only adult females. Sexual dimorphism in the Southern lineage was analyzed using MANOVA for the meristic characters, t-test for SVL and Ltot, and MAN-COVA, with SVL as a co-variate, for the remaining measures. These analyses were carried out in SPSS Statistics 17.0.0. (SPSS Inc., Chicago, IL, USA).

Mitochondrial diversity and phylogeography
Our comprehensive sampling ( Fig. 1A and Table S1) split Pseudopus apodus into three lineages in all analyses and using all sequence alignments (Fig. 2, S1, S2). These lineages correspond to the subspecies P. a. apodus, P. a. thracius, and to the Southern lineage (below described as a new subspecies; Fig. 2A,C,E). The lineage corresponding to the nominotypic subspecies is distributed from Crimea, through the north-eastern Black Sea coast and southern Russia, Georgia, Azerbaijan, Armenia, northern Iran, Turkmenistan to Tajikistan, Uzbekistan and Kyrgyzstan. The lineage corresponding to the subspecies P. a. thracius inhabits the Mediterranean regions of the Balkan Peninsula, and western and northern Turkey. The Southern lineage inhabits southernmost Turkey, and northern and central Israel as well as western Syria and Lebanon, though we have no mtDNA data from these countries (but see data from Syria in Microsatellite genotyping).
The network visualization from the concatenated dataset corresponded to the phylogenetic division into three groups, P. a. apodus, P. a. thracius, and the Southern lineage. We found seven haplotypes of P. a. apodus forming a star-like pattern with shallow genetic substructure, four distinct haplotypes of P. a. thracius, and eight haplotypes of the Southern lineage, showing deep intra-lineage diversity. More than 40 mutation steps were recognized between haplotypes of P. a. apodus and the Southern lineage, and more than 33 between P. a. thracius and the Southern lineage (Fig. 2B).
In the 727 bp long fragment of ND2 (32 sequences) used to build the haplotype network, we found four haplotypes of P. a. apodus, two haplotypes of P. a. thracius, and six haplotypes of the Southern lineage. The most distant in-group haplotypes (P. a. apodus and the Southern lineage) spanned 31 mutation steps. Between the most distant haplotypes of P. a. apodus and P. a. thracius there were 19 mutation steps. A star-like pattern was partially found in P. a. apodus haplotypes (Fig. 2D). The highest number of mutations and the most divergent haplotypes were found in the Southern lineage -13 mutations steps between populations from southern Turkey and central Israel (and eight within Israel; see details below; Fig. 2D).
In the 423 bp long fragment of Cyt b (85 sequences) used for the haplotype network, we found seven P. a. apo dus haplotypes, ten haplotypes of P. a. thracius, and four haplotypes of the Southern lineage. Between the most distant in-group haplotypes (P. a. thracius and the Southern lineage) there were 29 mutation steps, nine mutation steps between P. a. apodus and P. a. thracius, and 13 mutation steps between P. a. apodus and the Southern lineage. Pseu dopus a. thracius had the highest number of haplotypes, and all lineages showed star-like pattern (Fig. 2F).
The dataset of the 727 bp long ND2 fragments of the Southern lineage revealed four, well-diverged haplogroups: two in southern Turkey and two in Israel. Three of these haplogroups comprised single haplotypes each. The most common haplogroup (19 sequences), detected in central and northern Israel, contains three haplotypes. All haplogroups are separated by at least six mutation steps. The most distant haplogroup was found on the southern border of the species range (at least eight mutation steps from the adjacent haplotypes; Fig. 3A).
In the 575 bp long Cyt b fragments dataset of the Southern lineage we detected lower haplotype separation with two main haplogroups separated by three mutation steps from each other. One haplogroup containing a single haplotype was detected in southern Turkey, another, containing four haplotypes and 21 sequences separated by one mutation step, was found in Israel (Fig. 3B).
The ranges of uncorrected p distances were as follows: ND2 + Cyt b )concatenated dataset, 1532 bp( 2.6% between P. a. apodus and P. a. thracius, 3.5% between apo dus and the Southern lineage, and 3.6% between thracius and the Southern lineage. In ND2 alone )727 bp(: 3.0% between apodus and thracius, 3.8% between apodus and the Southern lineage, 3.6% between thracius and the Southern lineage. With Cyt b )805 bp(: 2.6% between apodus and thracius, 3.3% between apodus and the Southern lineage, and 3.7% between thracius and the Southern lineage. In accordance with the phylogenetic and haplotype network analyses, the PCAs revealed three highly distinct clusters corresponding to P. a. apodus, P. a. thracius and the Southern lineage with non-overlapping 95% confidence intervals ) Fig. 4A(. Due to the similarity of sequences and haplotypes, the individual values appear only as few dots despite the large numbers of samples. First (PC1) and second (PC2) principal components explain 45.9% and 21.3% of the observed variance in the concatenated dataset, 44.9% and 17.4% in the ND2 dataset, and 26.3% and 20.7% in the Cyt b dataset, respectively.

Microsatellite genotyping
For detection of genetically homogeneous groups of individuals in Structure, we used two approaches: direct posterior probabilities [LnP)K(] and ad hoc statistic ΔK. Direct posterior probabilities increased sharply from K=1 to K=3 and then, with larger Ks, increased just slightly and reached a plateau. The statistic ΔK estimated the most likely number of clusters as K=2, followed by K=3. For higher Ks the analysis obtained much lower probabilities. Assuming K=2, all individuals from the Balkans and western Turkey were assigned to cluster 1, corresponding to P. a. thracius. Specimens from the Crimea, through the Caucasus region, Iran up to Central Asia were assigned to cluster 2, corresponding to P. a. apodus. Lizards from the Levantine region (Israel, Syria, and an adjacent part of southern Turkey) were either assigned to cluster 1 or showed admixture. Assuming K=3, Levantine specimens formed a separate cluster 3, corresponding to the Southern lineage (Figs. 1B and 4, Table S2). The only exceptions were two individuals from Biranit (Israel) and Misyaf (Syria), which showed admixed genomes in the admixture model but were assigned to the Southern lineage in the non-admixture model. PCA analysis of microsatellite data showed similar results ( Fig. 4B) to those obtained using mtDNA and the Structure clustering (K=3; Fig. 4C), presented a clear separation among the three sampling regions, i.e. 1) the Balkans and western Turkey (P. a. thracius), 2) Crimea eastward to Central Asia (P. a. apodus), and 3) the Levant  Table S1) of the genetically investigated populations of Pseudopus apodus levantinus ssp. nov. (ND2 & Cyt b) in the Levant with the haplotype networks for each DNA marker. Circle sizes in the haplotype networks are proportional to the relative frequency of haplotypes and small black circles represent missing haplotypes. The expected distribution range (orange) of the Levantine subspecies (see also Figure S4) follows Sindaco and Jeremcenko (2008) and .
To assess genetic diversity and differentiation, all samples were assigned to three groups corresponding to the Structure and PCA clusters. The percentage of polymorphic loci (P) was lower in P. a. apodus and P. a. thra cius )87.5% polymorphic loci( than in the Southern lineage, in which all microsatellite loci were polymorphic. Individuals from the Levant showed substantially higher genetic diversity [numbers of different )N A ) and private (N P ) alleles and observed (H O ) and expected (H E ) heterozygosity; Table 1] followed by lizards assigned to P. a. apodus and P. a. thracius, respectively.
Genetic differentiation, as measured by unbiased Nei's genetic distance (D S ) and F ST statistics, reached similar values between the Southern lineage and both of the subspecies ) Table 2(. The highest differentiation was found between P. a. thracius and P. a. apodus (D S = 0.796, F ST = 0.395). Average F ST summarized over all samples reached the value of 0.333. Genetic differentiation among all three groups was highly significant ) Table 2

Morphology Quantitative analyses of morphological data
Snout-vent length )SVL( differed significantly between subspecies and sexes. Males are longer than females, specimens of the Southern lineage are longer than P. a. apodus, and P. a. thracius specimens are the shortest (Tables 3 and 4). All other characters were therefore compared between the three groups while controlling for SVL and sex (no interactions were tested). All metric and meristic characters differed significantly among all three groups, except for the frontal-snout length (FSL ;  Table 3(. Morphological differences among the three groups revealed variable degrees of divergence between the Southern lineage and P. a. apodus and P. a. thracius. The Southern lineage has the longest pileus (PL) and the highest intermaxillar shield (IMH). It also has wider head (HW1 and HW2) and intermaxillar shield (IMW), and longer rudiment (RL) than P. a. apodus but does not differ in these characters from P. a. thracius. The head at the level of the suture between the third and fourth supraorbital scales (HW3) is wider in the Southern lineage than in P. a. apodus, and narrower than in P. a. thracius. The Southern lineage does not differ from P. a. apodus but has significantly longer heads than P. a. thracius (HL1 and HL2) and frontal shields (FL), and wider snouts (SW). The frontal shield (FW) is narrower in the Southern lineage than in both P. a. apodus and P. a. thracius. In meristic characters, the Southern lineage has the highest number of preanal scales (PAN), and the number of these scales increases with SVL. Members of the Southern lineage have more supralabials (SPL) than P. a. apodus and fewer than P. a. thracius for a given SVL. The Southern lineage has fewer dorsals (DSL), ventrals (VSL), and subcaudals (SCR) than P. a. apodus, and more ventrals (VSL), and subcaudals (SCR) than P. a. thracius (Table  3). A summary of the variation of the characters in P. a. apodus, P. a. thracius, and the Southern lineage is presented in Table 4.
The three PCA analyses with (1) both males and females; 2) only males; 3) only females, presented similar results. The first component explained 98.7, 98.6, and 98.9%, respectively, of total variance of the 17 analyzed characters, with the most important character being SVL (Table S5). Most of the specimens from Israel (the Southern lineage) had much higher PC1 scores than any specimen of the other subspecies, while the majority of P. a. apodus and P. a. thracius specimens had a combination of high scores on PC2 and low scores on PC1, not seen in the Southern lineage. However, there were areas of morphological overlap among all three lineages, especially between P. a. apodus, and P. a. thracius, which showed near-complete overlap (Fig. 5, Table S5).
We found no evidence for sexual dimorphism in meristic characters in the Southern lineage. None of the eight measures significantly differ between males and females and three characters (DST, VST, PAN) do not vary at all. The sexes do not differ in their overall size, but males have generally wider heads than females )significant in HW1-3, SW, and IMW; Table S6). In the PCA analysis of the sexual dimorphism of the Southern lineage the first axis explained 99.2% of the total variance of the four analyzed characters, and each analyzed character was presented as the most important in each component (HL1 in the first component; HW1 in the second; PL in the third one; and HW2 in the fourth component; Fig. S3).

Systematic account
Pseudopus apodus comprises three clearly genetically and slightly morphologically differentiated allopatric populations (Obst 1978, this study).  Table 3. Intercepts )in mm( of the analysis of co-variance )ANCOVA( models with subspecies and sexes main effects and SVL as a covariate (for details see Table S4).  Counts and measurements (in mm) are presented as minimum-maximum (mean±standard deviation); n = total number of specimens studied. Only intact and non-regenerated tails were considered in SCR, TL, and TotL (for details see Table S4). Two of them are taxonomically recognized as subspecies: P. a. apodus (Pallas, 1775) and P. a. thracius (Obst, 1978). The type locality of P. apodus (Pallas, 1775) (Pallas, 1814: 33 p.). This Pallas' auto-correction provided ground for the assumption that the information about the occurrence of Lacerta apoda in the Volga-Ural interfluve )Naryn steppe( was erroneous. Obst (1978) probably overlooked the above mentioned original range description provided by Pallas and because he knew that P. apodus did not occur in the Naryn Steppe (northern coast of the Caspian Sea; today's Kazakhstan laying outside of the species range according to the current knowledge), restricted the type locality to "Terek-Gebiet" (Obst 1978: 136 p.   and by us, and to microsatellite data from the Bulgarian Black Sea coast (Fig. 1, Table S1). Our data also confirm that P. a. thracius has a wide distribution (Fig. 1). This range also includes the type locality of the name Pseudopus Durvilii Cuvier, 1829, which is "L'Archipel", Greece. Therefore, we designate the name Oph isaurus apodus thracius Obst, 1978  The proposed common name in English, Hebrew, and Arabic: Levantine Glass Lizard, ‫לבנטיני‬ ‫חורש‬ ‫,קמטן‬ ‫,الرشق‬ respectively.
Comparisons. Pseudopus apodus levantinus ssp. nov. can be distinguished from P. a. apodus and P. a. thracius by having 10 preanal scales (PAN) in all morphologically examined specimens (vs. 5-11, with the mean of 8 in P. a. apodus, and 5-8 and mean of 6 in P. a. thracius). Other differences can be used in combination to support differences between P. a. levantinus ssp. nov., P.  Fig. 6, Table S7), specimen in a good state of preser vation in 70% ethanol, midbody oval, and robust. All measurements of the holotype were tak- The hind limb rudiments are well visible with a length 7.42 mm. The complete tail is clearly differentiated from the body by the cloacal opening and 10 preanal scales on the ventral side. There are 154 subcaudal scale rows (SCR) overall, which is very similar to the specimens with regenerated tails. It has 106 longitudinal dorsal scale rows (DSL), and 121 ventral longitudinal scale rows (VSL).
The coloration of the holotype in life was not recorded. The coloration of the holotype in preservative is brownish or slightly orange with some of the scales on the body that have darker coloration creates an impression of tiny dark spots on the body. The head is lighter than the body, which is especially apparent on its dorsal side. The dorsal side of the body is light brown to gray, again with the impression of darker spots on the tail.
Variation. Details on variation among the type specimens of P. a. levantinus ssp. nov. are presented in Table  S7. The overall morphology and coloration of the paratypes (TAU 17076, 17235, 17895) are very similar to that of the holotype; TAU 19403 was recently collected and has an evident dark-brown coloration pattern (TAU 19403, Fig. 8).
Etymology. No name is available for the glass lizards from the Levantine region. We hence suggest a new name, Pseudopus apodus levantinus, as a reference to the isolated and allopatric distribution of this subspecies exclusively in the Levant region. This region covers present-day )western( Syria, Lebanon, )north-western( Jordan, Israel, West Bank, Cyprus, and most of Turkey south-east of the middle Euphrates, which is almost identical to the known distribution range of the new subspecies. The term "Levant" is derived from the Italian "Levante", meaning "rising" and implying the rising of the Sun in the east as an equivalent to the Arabic "al-Mashriq" ‫)المشـرق(‬ and the Hebrew "Mizrāḥ" ‫ח)‬ ‫ְרָ‬ ‫ז‬ ‫)מִ‬ both meaning "east".
Distribution. According to the genetic data, P. a. levan tinus ssp. nov. occurs in southern Turkey, western Syria, northern and central Israel this study). although we did not examine specimens from these regions. We therefore hypothesize that it ranges from southern Turkey (Hatay province) to central Israel (throughout the Mediterranean zone but excluding the deserts to the south and east; Figs. 3 and S4). Its northern distribution edge is presumably south of the Nur Mountains, the biogeographic barrier for the split between Anatolian and Levantine populations of other species Tamar et al. 2015;Jablonski and Sadek 2019;Šmíd et al. 2021). This, however, needs further investigation as there is a possibility that the new subspecies could be found further north in Mersin or Adana provinces (Figs. 1, S4 and Sindaco and Jeremcenko 2008; see also Hofmann et al. 2018 Fig. S4; Disi and Amr 1998;Roll et al. 2017, the Steinhardt Museum records).
Habitat and ecology. The new subspecies is known in Mediterranean habitats of the Levant (see the type locality Fig. 7), preferring relatively well-shaded light woodlands (maquis and garrigue), dry and warm hillsides, stream banks, and agricultural fields. It is found from ~50 m below sea level )in Hefzibah, Israel; 32. 52°N nov. is a leaf-litter diurnal carnivore (Bar and Haimovich 2012). It is oviparous and females typically lay a single clutch comprised of 3-12 eggs (data from the records of the Zoological Research Garden, Tel Aviv University). Its habitat characteristics, ecology, and diet were studied in Jordan and presented by , however data from other countries (see Distribution) are scarce (see Arbel 1984;In den Bosch et al. 1998;Bouskila and Amitai 2001;Bar and Haimovich 2012;Werner 2016).
Diversity. The specimens analyzed in this study showed certain variation in both genetic and morphological markers. Some intra-lineage variation in P. a. levantinus ssp. nov. can be found in ND2 and Cyt b sequences as well as in microsatellites (see details in the Results). The haplotype pattern of nuclear genes PRLR and RAG1 indicates incomplete lineage sorting among all three subspecies, with one recorded heterozygote of P. a. levantinus ssp. nov. . The calculated genetic distances among the three main mtDNA lineages )~2.6-4.0%( and the estimated time of divergence in Pseudopus mostly based on the nuclear loci (~8 to 3 Mya) support the subspecies status of the Levantine population (Jandzik et al. 2018, Lavin andGirman 2019). This is in accordance with subspecific taxonomy of the closely related species Anguis colchica (approx. split between 2.8-2.5 Mya, genetic distances ~3.1-4.7% among defined subspecies/ lineages based on ND2 mtDNA fragment; , while the genetic distance among the Anguis species are significantly higher )7.0% between A. fragilis and A. colchica and even 9.2% between A. fragilis and A. veronensis; . We also recorded some variation in the measurements and coloration both among the three subspecies (Figs. 5,9) and within P. a. levantinus ssp. nov. Males of P. a. levantinus ssp. nov. differ from females by having relatively wider heads (HW1-3, SW, IWM; Table  S6). The coloration of this subspecies is similar to the nominotypic subspecies ) but more data are needed to allow for drawing any more robust conclusions. A melanistic individual has been recorded in the Levantine population (Jablonski and Avraham 2018).
Conservation. Based on the data presented here, the distribution range of P. a. levantinus ssp. nov. covers approxi mately 30,000 km 2 . Together with human overpopulation and accelerated development in the Mediterranean parts of the Levant, high traffic density )most individuals nowadays brought to the Steinhardt Museum of Natural History in Tel Aviv are road-kills), development of mass tourism, extensive use of pesticides in the agricultural areas, proliferation of human commensals such as domestic cats, cattle egrets, dogs, rats, and golden jackals, and challenging political situation have potential to worsen the conservation status of this endemic subspecies. Pending a formal assessment, we preliminary recommend the

IUCN category of Vulnerable (VU) based on the criteria
A2c,e and strongly encourage further surveys benefiting from international collaboration allowing to open a dialogue across the conflict zones )e.g. EcoPeace Middle East, http://ecopeaceme.org; Roulin et al. 2017). Interestingly, recent data suggest that populations of P. apodus from the Levant were likely a part of modern human diet for millennia )Natufian culture, around 13,050 to 7,550 BC), which provides a rare piece of evidence of a longterm civilization pressure on the local biota (Lev et al. 2020).

Pseudopus apodus: a species of a monotypic genus with three subspecies
Pseudopus apodus is the single extant species of the genus Pseudopus that has been known as a relatively common member of the fossil assemblages from the Early Miocene or possibly even the Oligocene (Klembara et al. 2010;Klembara 2015;Čerňanský et al. 2016;Villa and Delfino 2019(. It has been suggested that the species gradually replaced a larger, ecologically more plastic species P. pannonicus during the Pleistocene (Klembara et al. 2010). The possibly lower ecological plasticity of P. apodus has been used to explain its more limited distribution in comparison to its closest relatives from the genus Anguis .
The simplified morphology and lack of conspicuous markings may be among the reasons why the species had been considered monotypic until F. J. Obst described the subspecies P. a. thracius based on a slightly distinctive morphology (Obst 1978). Recently, an additional support for recognition of P. a. thracius has been provided by a molecular-phylogeographic study . These authors also revealed a genetically distinct lineage in the Levant, which could potentially represent morphologically cryptic endemic subspecies. After adding further evidence from both molecular-genetic markers and from morphology, we described it here as P. a. levantinus ssp. nov. The newly described taxon meets the criteria of the subspecific category according to which the subspecies should occupy a distinct geographic range, and should have distinct hereditary characters including mtDNA ha plo types, nuclear alleles and ideally, but not necessarily, morphological traits (Patten 2015, Kindler and Fritz 2018, Hillis 2020. In contrast to species, subspecies are not reproductively isolated and are capable of extensive gene flow. Although we do not have evidence for present or historical gene flow between the three P. apodus lineages (see below), we opted for subspecific and not specific rank for it taking into account the expected time of divergence and genetic distances in relation to the divergence of closely related Western Palearctic species within the family Anguidae (see also Diversity).
This new subspecies is currently known from the area south of the Nur Mountains in Turkey, to the Mediterranean regions of Syria and Israel. While we lack direct genetic evidence, it is highly probable that also the populations from the Mediterranean coast of Lebanon and from Jordan belong to P. a. levantinus ssp. nov. Not so clear are the taxonomic relationships of the populations from the south-central Turkey around the Gulfs of Iskenderun and Mersin (particularly Mersin and Adana Provinces), the only area where two subspecies likely come into contact and, potentially, hybridize (Fig. 1). While the distribution shows an almost 400 km distribution gap from the closest P. a. thracius populations in the west (see  and Fig. 1), this area also lies across the Nur Mountains which are known to form a biogeographic boundary between lineages of several other reptile taxa Jablonski and Sadek 2019;Šmíd et al. 2021). Although we expect the new subspecies to occur in these parts of southern Turkey (i.e. similar phylogeographic pattern of Zamenis hohenackeri; Hofmann et al. 2018), we cannot exclude that the populations west of this mountain range are distinct (P. a. thracius; ( from those in the east and southeast, reflecting a possibly longer separation and range fragmentation. Alternatively, the occurrence of the nominotypic subspecies in the southern-central Turkey, is also possible, although less likely (see the biogeographical pattern of Bufo ver rucosissimus; Özdemir et al. 2020). The most enigmatic populations of P. apodus are, however, in north-eastern Iraq and mid-western Iran (Schmidt 1955;Reed and Marx 1959). They are isolated from the remaining species range (Levant, Turkey, Caspian region) by vast areas of unsuitable arid habitats of the Syrian Desert and Mesopotamia in the west and likely by Zagros Mountains montane habitats in the north (Fig. 1). While the range of the nominotypic subspecies is the closest to this Iraqi and Iranian populations, they could also potentially represent old relict populations or even a distinct evolutionary lineage. It would be thus extremely interesting to learn more about their phylogenetic affiliations. With our present knowledge we cannot exclude the possibility of the historical contact with subsequent extinction but the current species range seems to be fragmented (Fig. 1). This is congruent with our genetic data and, therefore, studying possible contact zones and inter-lineage hybridization in P. apodus offers a further challenge.
All three subspecies are genetically differentiated in both mtDNA and nDNA. Strongly supported subspecific lineages in the phylogenetic analysis and relatively long distances among mtDNA haplotypes (Figs. 2, 3) are mirrored in non-overlapping, well-supported and distant clusters of the highly polymorphic microsatellite genotypes (Fig. 4B). All unambiguously correspond to the three subspecies. Microsatellites, as fast evolving neutral markers, show better discrimination among closely related evolutionary lineages than most nuclear genes (e.g. Spilani et al. 2019), including RAG1 and PRLR used in our previous study ). In the case of the sequences of these two genes, two and one out of the three identified haplotypes, respectively, were shared be-tween two or even among all three subspecies (see Fig. 3 in ), which we interpreted as indicative of incomplete lineage sorting. This is not unusual among closely related evolutionary lineages of reptiles, which show monophyly in mtDNA but paraphyly in nuclear genes (e.g. Abreu et al. 2020). We discuss alternative biogeographic hypotheses explaining this incongruence below.
In addition to the genetic level differentiation, we provide morphological evidence supporting independent evolutionary history of the three subspecies. The new subspecies differs from P. a. apodus and P. a. thracius in several traits, most notably in the overall size. It attains significantly longer average and maximum lengths )Tables 3, 4). We note, however, that the maximum recorded size for the species (1,440 mm; Obst 1981) was reported from a specimen observed in Bulgaria (P. a. thracius).
In several other characters P. a. levantinus ssp. nov. differs from both remaining subspecies, while in some it is more similar to one or to the other (Tables 3, 4). These overlaps are visualized in the PCA plots (Fig. 5), in which P. a. levantinus ssp. nov. shows greater variation. Its morphospace, as manifested in the first two principle components, encompasses most of the variation existing within both P. a. apodus and P. a. thracius, and its range along both axes (especially PC1) is greater than the combined range of the other two subspecies. So far, the only available morphological analysis of the specimens that most likely represent this Levantine subspecies comes from Jordan )not confirmed by molecular-genetic data yet), where  found that the lizards are characterized by morphology intermediate between P. a. apodus and P. a. thracius. The results of our morphological analyses (and partly also of microsatellites) show a similar though not exactly the same picture. We posit that the overlap in morphology could be explained by higher variation within P. a. levantinus ssp. nov. However, overall size is a very plastic trait, often characterized by high variation even within genetically uniform groups. It can be influenced by a number of environmental variables and their interactions (Horvathova et al. 2013;Mikulíček et al. 2013;Pincheira-Donoso and Meiri 2013;Ortega et al. 2019;Slavenko et al. 2019). Interestingly, Glavaš et al. )2020( recently found striking size differences between two geographically relatively close populations of P. a. thracius in Croatia. The population attaining larger size was an island population from a more extreme environment than the mainland population, generally resembling the situation among P. a. levantinus ssp. nov. and the other two subspecies. The question whether and to what extent the distinct morphology of P. a. levantinus ssp. nov. is a result of environmentally driven plasticity or independent evolutionary history of the lineage remains to be answered and would require at least a correlative study between morphological and environmental variation in all three subspecies. Studying environmentally driven phenotypic plasticity of recent P. apodus populations could be very informative for re-evaluation of the taxonomic identity of very similar fossil species. These are mainly characterized by different sizes and osteolog-ical characters associated with size such as ridges, crests, and tubercles of the cranial bones or cranial and vertebral processes (e.g. P. pannonicus; Klembara et al. 2010Klembara et al. , 2019.

Historical biogeography of Pseudopus in the Levant
The Levant is a major biodiversity hotspot (Myers et al. 2000) and has the highest reptile richness (Sindaco and Jeremcenko 2008;Roll et al. 2017) and endemism (Ficetola et al. 2018;Gumbs et al. 2020) in the Western Palearctic. The actual Levantine reptile diversity still appears to be underestimated Šmíd et al. 2019, 2021Kornilios et al. 2020). Our description of a new taxon in a very common, large and conspicuous lizard species is a further evidence for this. Endemism in the Levant is often associated with complex paleogeographic events that occurred in the Levant mainly in the middle Miocene and separated Eurasia from Africa-Arabia (e.g. Tamar et al. 2016) with subsequent geological instability, volcanic activity, and environmental changes, mainly aridification )Edgell 2006(. All this probably facilitated the emergence of different evolutionary lineages in the Levant. The origin of endemic species in this region can be dated back to the Late Eocene (split of Latonia spp. and Discoglossus spp., or Parvilacerta spp. from Anato lolacerta spp.; Biton et al. 2013;Garcia-Porta et al. 2019) or the Miocene (split of Mediodactylus orientalis and M. kotschyi, Acanthodactylus beershebensis and A. ma cu latus, Phoenicolacerta kulzeri and P. laevis, or Xero typhlops syriacus and X. vermicularis; Tamar et al. 2015Tamar et al. , 2016Kornilios 2017;Kotsakiozi et al. 2018;Kornilios et al. 2020). Similarly, dating of some of the currently recognized lineages representing endemic subspecies places their splits to the Miocene (Zamenis h. hohenackeri and Z. h. tauricus; Hofmann et al. 2018), though the Pleistocene origin is more typical (Pelobates s. syriacus and P. s. boettgeri, Bufotes cf. sitibundus, and Platyceps c. collar is and P. c. rubriceps; Dufresnes et al. 2019a,c;Šmíd et al. 2021). Dating of the basal split among the P. apodus subspecies remains controversial.  placed the basal split to about 2-3 Mya (Plio-Pleistocene) and discussed strong incongruence between the genetic and fossil data in the split of the sister genera Anguis and Pseudopus, while Lavin and Girman (2019) using more robust genetic data (25 nuclear loci and three mtDNA genes) but limited taxon sampling (missing P. a. thracius) dated the split to the Late Miocene.  presented a biogeographic scenario in which P. a. apodus and P. a. thracius recently expanded to their ranges after they survived the Pleistocene climatic oscillations in one or only few refugia. Deeper divergences within P. a. levantinus ssp. nov. lead us to hypothesize that this lineage did not experience significant glacial retraction and postglacial expansion and that recent diversity is of an older origin than in the other two subspecies. Our new data confirm the highest within-group genetic variation of P. a. levantinus ssp. nov., i.e. the highest number of ND2 haplotypes and the highest microsatellite polymorphism and diversity (Figs. 3, 4, Tables 1, 2(. This distinct structure is also reflected in morphology, with higher variance in the new subspecies than in P. a. apodus and P. a. thracius (Fig. 5, Table 4). Several newly revealed, and distantly related, mtDNA haplotypes of P. a. levantinus ssp. nov. are separated by numerous missing haplotypes suggesting longer independent evolution (and/or sampling artefacts) in more stable environment, rather than repeated extinctions and expansions.
The molecular-genetic data we obtained from different markers, and by different methods, show some incongruence in detailed relationships among the three subspecies. The genetic distances, network, tree and PCA analyses analyses using mtDNA sequences (Figs. 1-4A) show that P. a. levanti nus ssp. nov. is more distant from both P. a. apodus and P. a. thracius than they are from each other (Figs. 2, 3). However, closer relationships of P. a. levantinus ssp. nov. to P. a. thracius than to P. a. apo dus based on microsatellite analyses (Fig. 4B, Table 1, 2) corroborate findings from the sequences of the nuclear genes , in which P. a. levantinus ssp. nov. and P. a. thracius share more nDNA haplotypes than any other combination of the remaining subspecies. This pattern suggests three alternative scenarios: i) it could be an artefact of incomplete lineage sorting, ii) sharing the nuclear haplotypes and higher variability in microsatellites could have resulted from the past hybridizations of P. a. levantinus ssp. nov. with the remaining two subspecies over large portions of the range, resulting in massive nuclear introgressions. Alternatively, iii) P. a. levantinus ssp. nov. could have split from P. a. apodus (since they have the lowest Nei's distance), after the latter diverged from P. a. thracius. Then the mtDNA diversity would be largely influenced by bottlenecks following the range restrictions. Neither of these scenarios is substantially supported by the existing data and a better picture of the biogeographic history of P. apodus awaits further study.