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 Jandzik et al. (2018) 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; Mikulíček et al. 2018) were amplified in two multiplex PCRs (multiplex 1: AF38-VIC, AF44-FAM, AF47-PET, AnFr21-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.

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) phylogenetic 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).

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⁶ iterations, following a burn-in period of 10⁴ 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.

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).

Figure 2. 

Phylogenetic relationships of Pseudopus apodus reconstructed from concatenated and single locus trees of ND2 & Cyt b sequences and presented as collapsed Bayesian trees (A, C, E) and haplotype networks (B, D, F). The numbers at the tree nodes represent Bayesian Posterior Probabilities/Maximum Likelihood Bootstraps showing branch support (for details see Figures S1 and S2). The numbers at the end of main lineages or with haplotypes show the sampling code (DJab acronym) or GenBank numbers (see also Table S1). Circle sizes in the haplotype networks are proportional to the relative frequency of the haplotypes and small black circles represent the missing haplotypes.

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. apodus 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. Pseudopus 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).

Figure 3. 

Maps showing the localities and sampling numbers related to the collection with DJab acronym or GenBank number (see 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 Jandzik et al. (2018).

We found clear differences in mtDNA polymorphism between the lineages with values of Hd = 0.800 and π = 0.18% for P. a. apodus (n = 6 sequences), Hd = 0.667, π = 0.09% for P. a. thracius (n = 3), and Hd = 0.458, π = 0.34% for the Southern lineage (n = 23) on the ND2 dataset and Hd = 0.722 and π = 0.28% for P. a. apodus (n = 28), Hd = 0.494, π = 0.17% for P. a. thracius (n = 35), and Hd = 0.271, π = 0.11% in the Southern lineage (n = 22) in the Cyt b dataset.

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 apodus 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.

Figure 4. 

Principal Component Analysis (PCA) of three lineages of Pseudopus apodus for mitochondrial DNA (A) and microsatellite data (B) with the Structure plot of posterior probabilities for K=3 (C). The oval outlines in PCAs represent 95% confidence intervals. 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 ND2 dataset, 26.3% and 20.7% in Cyt b dataset, and 10.8% and 6.7% in microsatellites, respectively. The plot of two principal components shows a clear separation between three groups of samples. Color scheme corresponds with the one used in Fig. 2: P. a. apodus (red), P. a. thracius (blue) and P. a. levantinus ssp. nov. (yellow).

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 (Southern lineage). First (PC1) and second (PC2) principal components explain 10.8% and 6.7% of the observed variance (Fig. 4B).

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. thracius (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). As shown by AMOVA, genetic variance within individuals was the highest (51%, DF = 155, SS = 244.00, MS = 1.574), followed by genetic variance among the three geographical regions (33%, DF = 2, SS = 199.92, MS = 99.96) and among individuals (15%, DF = 152, SS = 383.28, MS = 2.52).

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).

Figure 5. 

Principal Component Analysis (PCA) of the examined phenotypic characters among the adult specimens of Pseudopus apodus apodus (red dots), P. a. thracius (blue dots), and P. a. levantinus ssp. nov. (yellow dots): A - adult females and males together, B – females only, C – males only. The oval outlines correspond to 95% confidence intervals. For A, PC1 explains 98.7% of variance, PC2 0.59%; for B, PC1 explains 98.9% of variance, PC2 0.45%; for C, PC1 explains 98.6%, PC2 0.59% (see details in Table S5), respectively.

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).

Pseudopus apodus comprises three clearly genetically and slightly morphologically differentiated allopatric populations (Obst 1978, Jandzik et al. 2018, this study). 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), described under the name Lacerta apoda Pallas, 1775, was extensively but not precisely reported by Pallas (1775: 448–449 pp.) as a wide area between “… sabulosi Naryn inter Rhymnum et Volgam … Rarius etiam ad Sarpam rivum occurrit; sed frequentior in eodem, ex quo Sarpa fluit, deserto Kumano versus ipsum Kumam fluvium … circaque Terekum lecta fuit.” In 1776 (703 p.), Pallas published an abbreviated version of the report on the occurrence of the species: “Habitat in conuallibis herbidis deserti Naryn et ad Sarpam, Kumam, Terekum fluuios”. Additionally, for unclear reasons, Pallas subsequently restricted the range of the species to the Caucasus (with the Ciscaucasia, i.e. around Terek region) and the mountainous coast of the Crimean peninsula: “Ad omnem Caucasum, in hortis ad Terekum et in montosae orae Tauricae Chersonesi …” (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.). The valid type locality of P. a. apodus however, encompasses all the localities given by Pallas (1775, 1776) and mentioned above, i.e. from the Terek region to the Naryn Steppe the Naryn Steppe, included. Pallas (1776) clearly mentioned P. apodus from the more humid parts of Naryn Steppe though we wish to stress here that the current occurrence in this region has not been confirmed (Bakiev et al. 2020). Therefore, we used a sequence from a geographically closely related specimen from the locality of Voskresenkoe (approx. 25 km NNW of Gudermes, Tersko-Kumskaya, Terek River region, lowland Russia, 43.55°N, 46.38°E) of the ND2 marker (GenBank accession number AF085623, Macey et al. 1999; the specimen from the herpetological collection of the California Academy of Sciences [CAS]: 182911) and microsatellite data from the locality Bilbil’-Kazmalyar, Samur Forest, Russia-Dagestan, (41.90°N, 48.48°E, distance 230 km from the Terek River; see Fig. 1, Table S1) as proxies of the genetic lineage of the nominotypic subspecies from the type locality of the species. The type locality of P. a. thracius, described as Ophisaurus apodus thracius Obst, 1978 (holotype MTD 5464, J. Fritzsche leg. May 1972; see Fritz 2002 for the complete details on the type series), is on the Black Sea coast of Bulgaria, Primorsko, District Burgas (~42.27°N, 27.74°E). It is related to the sequences from Primorsko i.e. MF547732 (ND2) and MF547717 (Cyt b) analyzed by Jandzik et al. (2018) 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 Ophisaurus apodus thracius Obst, 1978 as nomen protectum and P. Durvilii as nomen oblitum according to the Article 23.9.1.2 of the ICZN (1999; see Appendix to this article). The animals from the Levant (analyzed here from the territory of southern Turkey, Syria, and Israel), represent a separate, allopatric lineage (Jandzik et al. 2018; Lavin and Girman 2019), genetically and morphologically differentiated from the two described subspecies (see details below). The Southern lineage of P. apodus, according to our results and under the current practices in reptile taxonomy and nomenclature (see Kindler and Fritz 2018; Fritz and Schmidtler 2020 based on Moritz 1994), can be taxonomically defined and we described it here as a new subspecies:

Pseudopus apodus levantinus ssp. nov. Jablonski, Ribeiro-Júnior, Meiri, Maza, Mikulíček, Jandzik

http://zoobank.org/D1621AC7-53A4-4CD5-A958-659510199169

The proposed common name in English, Hebrew, and Arabic: Levantine Glass Lizard, الشرق ,קמטן חורש לבנטיני, respectively.

Holotype

TAU-R 16895 (Fig. 6, Table S7), an adult male, collected on the 16th of May, 2014 by Erez Maza, at Giv’at Ada (גבעת עדה), Israel (32.52°N, 35.01°E, Fig. 7). The molecular-genetic data of the holotype are available for ND2 (MW400924), Cyt b (MW400903), and microsatellites (Fig. 2, Table S1, S2).

Figure 6. 

Holotype (TAU 16895) of Pseudopus apodus levantinus ssp. nov. from Giv’at Ada (הדע תעבג), Israel. A – dorsal view, B – ventral view, C – dorsal view of the head, D, E – lateral views, F – ventral view. For details and measurements see Tables S1, S2 and S7.

Figure 7. 

Habitat of Pseudopus apodus levantinus ssp. nov. at the type locality Giv’at Ada (גבעת עדה) in Israel (photo by Ilana Rosenstein).

Paratypes

TAU-R 17076 (MW400927, MW400906), an adult male, collected on 19^th May 2014 by Talia Oron, at Biranit (בירנית), Israel (33.05°N, 35.34°E); TAU-R 17235 (MW400928, MW400907), an adult male, collected on 6^th May 2015 by Ron Elazari, at Hadera (חֲדֵרָה), Israel (32.47°N, 34.89°E); TAU-R 17895 (MW400930, MW400909), an adult male, collected on 18^th April 2016 by Ofer Shimoni, on Mt. Gilboa’ (הר הגלבוע), Israel (32.45°N, 35.43°E); TAU-R 19403, an adult male, collected on 16^th May 2019 by Amir Arnon, at Ramat HaNadiv (רמת הנדיב), Israel (32.55°N, 34.95°E) (Figs 2, 8, Table S1, S2, S7).

Figure 8. 

Paratypes of Pseudopus apodus levantinus ssp. nov. from Israel showing A – dorsal part of the body, B – ventral part of the body, C – dorsal part of the head, D – ventral part of the head. For the locality details and measurements see Tables S1, S2, and S7.

Diagnosis

A large Pseudopus (up to 610 mm snout-vent length, 1,367 mm total length and a mass of 1,100 g) that can be distinguished from the other two subspecies by a combination of the following characters (means followed by standard deviations; Table 4): (1) preanal scales (PAN) 10; (2) long body (SVL; 490.22 mm ± 76.15); (3) long tail (TL; 600.70 mm ± 107.04); (4) long head (HL1; from the ear aperture to the tip of the snout, 47.6 mm ± 7.99); pileus length (PL; mean 41.98 mm ± 6.78); and distance from the frontal shield to the tip of the snout (FSL; 17.27 mm ± 2.96); (5) wide head (HW1; maximum width, 32.04 mm ± 6.34); and width at the level of the posterior edge of the eyes orbits (HW2; 28.21 mm ± 4.95); (6) relatively long limb rudiments (RL; 5.55 mm ± 1.53) (see and Tables 3 and 4 for comparisons to P. a. apodus and P. a. thracius).

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. a. apodus, and P. a. thracius: longer body (SVL; 610 mm maximum length and mean of 490.22 mm in P. a. levantinus ssp. nov., vs. 487 mm maximum and mean of 387.61 mm in P. a. apodus, and 410 mm maximum and mean of 360.06 mm in P. a. thracius); longer distance between the ear aperture and the tip of the snout (HL1; 63.75 mm maximum length and mean of 47.60 in P. a. levantinus ssp. nov., vs. 58 mm maximum and mean of 38.61 mm in P. a. apodus, and 49.50 maximum and mean of 39.55 mm in P. a. thracius); longer pileus (PL; 54.79 mm maximum length and mean of 41.98 in P. a. levantinus ssp. nov., vs. 46.30 mm maximum and mean of 35.82 mm in P. a. apodus, and 28.80 mm maximum and mean of 36.93 mm in P. a. thracius); longer distance between the anterior frontal scale to the tip of snout (FSL; 23.39 mm maximum length and mean of 17.27 mm in P. a. levantinus ssp. nov., vs. 18 mm maximum and mean of 13.90 mm in P. a. apodus, and 15.98 mm maximum and mean of 12.58 mm in P. a. thracius); wider head (HW1; 44.90 mm maximum width and mean of 32.04 mm in P. a. levantinus ssp. nov., vs. 36.5 mm maximum and mean of 24.47 mm in P. a. apodus, and 29.50 mm maximum and mean of 24.35 mm in P. a. thracius); wider distance between the posterior edge of the orbits (HW2; 39.33 mm maximum width and mean of 28.21 mm in P. a. levantinus ssp. nov., vs. 28.76 mm maximum and mean of 19.44 mm in P. a. apodus, and 26.30 mm maximum and mean of 22.33 mm in P. a. thracius); and longer limb rudiments (RL; 8.56 mm maximum length and mean of 5.55 mm in P. a. levantinus ssp. nov., vs. 6.55 mm maximum and mean of 3.5 mm in P. a. apodus, and 4.8 mm maximum and mean of 3.31 mm in P. a. thracius). The remaining differences are presented in Tables 3–4.

Description of the holotype

An adult male (TAU 16895; Fig. 6, Table S7), specimen in a good state of preservation in 70% ethanol, midbody oval, and robust. All measurements of the holotype were taken post-preservation. Snout-vent length 532 mm, tail length 602 mm, weight 900 g. The number of transversal dorsal scale rows at midbody 12 (DST), the number of transversal ventral scale rows 10 (VST). Head large, length (HL1) 59.01 mm, clearly distinct from the neck. The pileus length is 50.52 mm, maximum head width 42.25 mm. The right part of the mouth is slightly open. The number of supralabials on both parts of the head is 11. Head and body scales smooth. Slightly keeled scales are visible on the ventral part of the tail. Rostral curved toward the top of the head. Frontal scale well visible and big with the length 17.73 mm and width 13.15 mm. Eyes oval, both closed. The ear and nose openings are visible. 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. levantinus ssp. nov. occurs in southern Turkey, western Syria, northern and central Israel (Jandzik et al. 2018; this study). Based on the published distribution data, the new subspecies could also be expected in the Mediterranean regions of Lebanon and north-western Jordan (In den Bosch et al. 1998; Disi et al. 2001; Hraoui-Bloquet et al. 2002; Rifai et al. 2005; Sindaco and Jeremcenko 2008), 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 S4b). 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 (Jandzik et al. 2013; 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 Baran et al. 1988; Sindaco and Jeremcenko 2008; see also Hofmann et al. 2018). The southern limit of the subspecies range is near the border of the Negev desert, around Lahav (31.38°N, 34.87°E), and Kisufim (31.37°N, 34.40°E) in Israel, and Dhibam (31.47°N, 35.80°E) in Jordan (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, 35.42°E) to the hilly sub-mountain areas in Israel, Lebanon, Syria, and Turkey (elevation maximum ~1000 m; In den Bosch et al. 1998; Disi et al. 2001; Rifai et al. 2005; Werner 2016; the Steinhardt Museum records). Pseudopus a. levantinus ssp. 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 Rifai et al. (2005), 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. (Jandzik et al. 2018). 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 and Girman 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; Gvoždík et al. 2010; Jablonski et al. 2016), 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; Gvoždík et al. 2013). 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 (Rifai et al. 2005) 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).

Figure 9. 

Color and pattern variation in the Pseudopus apodus subspecies: A – adult male of P. a. apodus in a typical steppe habitat from Samarkand, Uzbekistan (photo by Daniel Jablonski). B – juvenile individual of P. a. apodus from Kyz-Kermen, Bakhchisarayi, Crimea (photo by Oleg V. Kukushkin). C – adult female of P. a. thracius from Dadia, Greece (photo by D. Jablonski). D – juvenile individual of P. a. thracius from National Park Paklenica, Croatia (photo by D. Jablonski). E – adult female of P. a. levantinus ssp. nov. from Antakya, Turkey (photo by David Jandzik). F – juvenile individual of P. a. levantinus ssp. nov. from Ness Ziyona, Israel (photo by David David).

Conservation

Based on the data presented here, the distribution range of P. a. levantinus ssp. nov. covers approximately 30,000 km². 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 long-term civilization pressure on the local biota (Lev et al. 2020).