Using multilocus approach to uncover cryptic diversity within Pseudotrapelus lizards from Saudi Arabia.

The agamid Pseudotrapelus lizards inhabit the mountainous areas of the Arabian Peninsula and eastern North Africa. Currently six Pseudotrapelus species are recognised, though diagnostic morphological characters are still lacking, creating great difficulty in describing new species. Recently, two specimens of Pseudotrapelus were collected from the vicinity of Riyadh in central Saudi Arabia, an area that was not sampled in previous phylogenetic studies. In here we used both mitochondrial and nuclear data to investigate the phylogenetic position of the new samples, and assess their phylogenetic relationships with the other recognised species of Pseudotrapelus from across the distribution range of the genus. We used a multilocus approach of haplotype networks, concatenated datasets and species trees, performed mitochondrial and nuclear species delimitation analyses, and estimated divergence times. In general, our results support previous molecular studies and uncover the presence of cryptic diversity within Pseudotrapelus. The phylogenetic structure of the genus is of two major clades and within them seven distinct, delimited phylogenetic groups belonging to the six recognised species and the seventh to the individuals from Riyadh. The Riyadh specimens were distinct in all analyses performed. We suggest that the new specimens from the Riyadh area are a distinct lineage, forming a clade with their phylogenetic relatives, P. sinaitus and P. chlodnickii. The clade formed by these three species diverged during the Late Miocene around 6.4 Ma, with cladogenesis possibly facilitated by vicariance and isolation caused due to climatic fluctuations and the progression of sandy areas. Our results suggest further morphological research is necessary to revise the taxonomic status of this lineage and of the entire genus.

Introduction

Pseudotrapelus Fitzinger, 1843 are rock-dwelling agamid lizards occupying various rocky habitats in hilly and mountainous areas, including well vegetated wadis and slopes, barren rocky hillsides, and boulder-strewn plains (Arnold, 1980, Baha El Din, 2006, Disi et al., 2001, Gardner, 2013). These lizards are distributed throughout the mountainous areas of the Arabian Peninsula and the Red Sea (Fig. 1; Tamar et al., 2016a and see references therein): from the Hajar Mountains in Oman and the United Arab Emirates (UAE) through the montane ridges of southern Oman, Yemen, and western and central Saudi Arabia, across southern Israel northwards to southern Syria, and westwards through the Sinai Peninsula to northen Sudan and western Eritrea.

Fig. 1

Sampling localities of the recognised Pseudotrapelus specimens used in this study (circles), and the new samples from the vicinity of Riyadh, Saudi Arabia (star). Identification and localities of the known species correlate to specimens in Table S1 in Tamar et al. (2016a). Colours correspond to those in Fig. 3, Fig. 4.

Fig. 3

Phylogenetic relationships within Pseudotrapelus. Colours correspond to those in Fig. 1, Fig. 4. (A) Bayesian inference phylogenetic tree generated from the complete concatenated dataset. Support values are indicated near the nodes (Bayesian posterior probabilities/ML bootstrap). The results of the species delimitation analyses are represented by coloured bars (GMYC, I–XI and BPP, A–G). Sample codes correlate to specimens in Table S1 in Tamar et al. (2016a) and the samples from this study, the two individuals from Riyadh, Saudi Arabia (IBECN6252, IBECN13348) and Acanthocercus adramitanus from Yemen (JEM135). (B) Unrooted haplotype nuclear networks. Circle size is proportional to the number of alleles.

Fig. 4

Time-calibrated species tree of Pseudotrapelus. Mean age estimates are indicated near the nodes with bars representing the 95% highest posterior densities. White circles denote nodes with posterior probability values ≥0.95. Colours correspond to those in Fig. 1, Fig. 3.

To date, Pseudotrapelus includes six recognised species (Fig. 1; see Tamar et al., 2016b, Uetz et al., 2019): (i) P. aqabensis, from western Saudi Arabia, southern Israel, Jordan, and the Sinai Peninsula; (ii) P. chlodnickii, from Sudan and Egypt, including the Sinai Peninsula; (iii) P. dhofarensis, from the mountains of southern Oman; (iv) P. jensvindumi, endemic to the Hajar Mountains of Oman and the UAE; (v) the recently resurrected P. neumanni, from the mountains of southern Yemen and Saudi Arabia; and (vi) P. sinaitus, from Egypt to southern Syria and northern Saudi Arabia.

Throughout the long systematic history of Pseudotrapelus it was thought the genus included a single species, P. sinaitus, with several morphological forms (Anderson, 1896, Anderson, 1901, Arnold, 1980, Baha El Din, 2006, Fritz and Schütte, 1988, Schätti and Gasperetti, 1994). The conservative morphology of these lizards greatly complicated their systematics. Despite the flurry of new species descriptions, none included an adequate morphological revision, resulting in insufficient diagnostic characters to identify specimens or describe new species (see Tamar et al., 2016a and references therein).

In this study we aim to investigate the phylogenetic position and relationships of two important recently collected specimens of Pseudotrapelus from the surroundings of Riyadh in central Saudi Arabia. These individuals were collected from an area that was not sampled in previous studies and therefore of unknown phylogenetic position. We therefore used DNA sequences of four markers (two mitochondrial and two nuclear gene fragments) to investigate their phylogenetic relationships. We performed multilocus phylogenetic analyses of concatenated datasets, species trees, and haplotype networks, as well as species delimitation analyses, and re-estimated the divergence times and genetic diversity between and within all Pseudotrapelus species. We also provide several morphological characteristics of these specimens.

Materials and methods

The two Saudi Arabian Pseudotrapelus individuals were collected from the vicinity of Riyadh (see Fig. 1). The vouchers and their tissue samples were deposited in Salvador Carranza’s collection at the Institute of Evolutionary Biology (IBE), Barcelona, Spain. Specimen with IBE collection code IBECN6252 (field code 9642X; Fig. 2A), an adult male, was collected on March 25th 2016, from a hill in a wadi, north-west to Thumamah, 25.592N 46.401E, 668 m elevation (Fig. 2B and C). Specimen with IBE collection code IBECN13348 (field code 9689X; Fig. 2D), adult female, was collected on April 27th 2018, from the foothill of Jebel Baloum, 23.699N 46.173E, 803 m elevation (Fig. 2E and F). GenBank accession numbers of the two specimens are MK176908–MK176914.

Fig. 2

General appearance and habitats of Pseudotrapelus specimens from the vicinity of Riyadh, Saudi Arabia. (A) IBECN6252; (B, C) north-west to Thumamah. (D) IBECN13348; (E, F) the foothill of Jebel Baloum. Photos by Laurent Chirio.

The monophyly of Pseudotrapelus and the close phylogenetic relations with Acanthocercus have been shown in Tamar et al. (2016a). Therefore, our dataset included a total of 100 specimens: the two Riyadh individuals, 92 specimens of the six recognised species of Pseudotrapelus (ten of P. aqabensis, six of P. chlodnickii, 22 of P. dhofarensis, 31 of P. jensvindumi, three of P. neumanni, 20 of P. sinaitus), and six outgroup specimens of Arabian Acanthocercus (one specimen of Acanthocercus adramitanus from Yemen, code JEM135, was additionally sequenced in this study; GenBank accession numbers for this sample are MK182721–MK182723). The remaining sequences of Pseudotrapelus and Acanthocercus were retrieved from GenBank from the study of Tamar et al. (2016a).

DNA of alcohol-preserved tissue samples was extracted using the SpeedTools Tissue DNA Extraction kit (Biotools, Madrid, Spain). The final dataset comprised four gene fragments that were previously used in Tamar et al. (2016a) of a concatenated length of 2219 bp: two mitochondrial, the ribosomal 16S rRNA (16S; ∼503 bp) and the NADH dehydrogenase subunit 4 (ND4; 681 bp), and two nuclear, the oocyte maturation factor Mos (c-mos; 372 bp) and the melano-cortin 1 receptor (MC1R; 663 bp). Primers and PCR conditions were the same as detailed in Tamar et al. (2016a). Chromatographs were checked, assembled and edited using Geneious v.7.1.9 (Biomatter Ltd.). For the nuclear genes heterozygous positions were coded according to the IUPAC ambiguity codes and no recombination was detected using SplitsTree v.4.14.5 (Huson and Bryant, 2006) (P >0.3 in PhiTest; Bruen et al., 2006). No stop codons were detected in the protein-coding genes (ND4, c-mos, MC1R). We aligned the sequences for each marker using MAFFT v.7.3 (Katoh and Standley, 2013). We calculated inter- and intraspecific uncorrected p-distance of the Pseudotrapelus species for 16S and ND4, with pairwise deletion, in MEGA v.7.0.14 (Kumar et al., 2016).

We performed phylogenetic analyses using the complete concatenated dataset of the four markers. Partitions and substitution models were defined by PartitionFinder v.2 (Lanfear et al., 2016) with the following parameters: linked branch length; BEAST models; BIC model selection; greedy schemes search algorithm; single data block for 16S, and by codons for the protein-coding genes. We treated alignment gaps as missing data, and the nuclear gene sequences were not phased. We analysed the complete concatenated dataset with the following partitions and substitution models: 16S + ND4_1 (HKY + G), ND4_2 (HKY + G), ND4_3 (GTR + G), c-mos_1 + c-mos_2 + MC1R_1 (JC + I), c-mos_3 + MC1R_2 (HKY + I), and MC1R_3 (HKY + G). Phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian (BI) methods. We conducted the ML analysis in RAxML v.8.1.2 as implemented in raxmlGUI v.1.5 (Silvestro and Michalak, 2012), with the GTRGAMMA model, 100 random addition replicates, and 1000 bootstrap replicates. The BI analysis was conducted in BEAST v.1.8.4 (Drummond et al., 2012) with the following priors (otherwise by default): partitions and models as detailed above; unlinked parameter values for clock and substitution models, linked trees; Coalescent tree model; random starting tree; base substitution parameter (0–100); alpha prior uniform (0–10); uncorrelated relaxed clock for the mitochondrial partitions and strict clock for the nuclear partitions (uniform distribution; mean 1, 0–1). Three individual runs of 108 generations were carried out with sampling at intervals of every 104 generations. We assessed posterior trace plots and effective sample size values in Tracer v.1.6 (Rambaut et al., 2014), and LogCombiner and TreeAnnotator were used to infer the ultrametric tree after discarding 10% of the trees as burn-in.

We inferred genealogical relationships among the haplotypes of Pseudotrapelus species using the phased nuclear sequences. To resolve heterozygous sites we first used the on-line web tool SeqPHASE (Flot, 2010) to convert the input files, and used the software PHASE v.2.1.1 (Stephens et al., 2001, Stephens and Scheet, 2005) to resolve phased haplotypes using default settings, apart from the phase probabilities of 0.9 for c-mos and 0.5 for MC1R. We then used the TCS statistical parsimony network approach (Clement et al., 2000) implemented in the software PopART (Leigh and Bryant, 2015, Clement et al., 2000) to generate the networks.

To evaluate the number of distinct mitochondrial and nuclear lineages within Pseudotrapelus we applied two methods of species delimitation using the ingroup dataset only. To evaluate mitochondrial divergence within Pseudotrapelus we carried out the GMYC analysis (Pons et al., 2006) implemented in R (R development Core Team, 2018) using the ‘splits’ package (Ezard et al., 2009) and applying a single threshold algorithm. For this analysis, we used the concatenated mitochondrial dataset and reconstructed a tree using BEAST with partitions and models defined by PartitionFinder (with parameters as previously detailed). Partitions, models, and other priors were as described above. We carried out three individual runs of 5 × 107 generations with sampling at intervals of every 5 × 103 generations. To evaluate nuclear divergence within Pseudotrapelus we used the software BPP v.3.3 (Rannala and Yang, 2003, Yang and Rannala, 2010) with the full-length phased nuclear loci only. We carried out two analyses with “species” assignation based on the GMYC entities: (i) Conducting species delimitation analyses using a fixed guide tree (the topology of the BI tree using the complete concatenated dataset). (ii) Performing a joint analysis of species delimitation while estimating the species tree (Yang, 2015). For these analyses, algorithms 0 and 1 were used assigning each species delimitation model equal prior probability. As prior distributions on the ancestral population size (θ) and root age (τ) can affect the posterior probabilities for models (Yang and Rannala, 2010) we tested four different combinations of priors: θ = G (1, 10), τ = G (1, 10); θ = G (2, 2000), τ = G (2, 2000); θ = G (1, 10), τ = G (2, 2000); θ = G (2, 2000), τ = G (1, 10). The locus rate parameter that allows variable mutation rates among loci was estimated with a Dirichlet prior (α = 2). Since our dataset was autosomal only, the heredity parameter that allows θ to vary among loci was set as default. We ran each of the rjMCMC analysis twice to confirm consistency between runs, each run for 5x105 generations with 10% discarded as burn-in. We considered probability values ≥ 0.95 as strong evidence for speciation.

Divergence times for Pseudotrapelus were estimated using a multilocus coalescence-based Bayesian species tree using *BEAST (v.1.8.4; Heled and Drummond, 2010). We used jModelTest v.2.1.7 (Darriba et al., 2012, Guindon and Gascuel, 2003) to select the best model for each marker under the Bayesian information criterion (BIC): 16S (HKY + G), ND4 (TrN + I + G), c-mos (K80 + I), MC1R (TrN + I). In this analysis, using the Pseudotrapelus dataset only, we set the average sequence evolution rates estimated for the agamid genus Phrynocephalus according to Pang et al. (2003). Tamar et al. (2016a) also used these rates and compared them to an analysis conducted using different calibration points; both analyses resulted in almost identical dates (see Table 1 in Tamar et al., 2016a). The rates and priors were: 16S (0.0073–0.0132 substitutions/site/million years) and ND4 (0.0113–0.0204 substitutions/site/million years), uncorrelated relaxed clock with the ucld.mean prior of 16S and ND4 set to Normal (initial 0.01, stdv 0.0015 and initial 0.015, stdv 0.0025, respectively), and a strict clock for c-mos and MC1R with the clock.rate prior set to Uniform (initial 0.001, 0–0.0204). Other priors were as described above, apart from unlinked parameter values for clock, substitution models and trees (linked trees for the mtDNA partitions), the Yule process species tree prior, random starting trees, and ploidy type autosomal (mitochondrial for the mtDNA tree). Three individual runs of 108 generations were carried out with sampling at intervals of every 104 generations. Stationarity was assessed with Tracer, and LogCombiner and TreeAnnotator were used to generate the ultrametric tree.

Table 1

Uncorrected genetic distances (p-distances) of the mitochondrial markers between (16S lower-left; ND4 upper-right) and within (bold; 16S/ND4) all Pseudotrapelus species.

	1	2	3	4	5	6	7
1. P. aqabensis	0.9/0.4	19	11.8	15.2	11.2	18.6	18.9
2. P. chlodnickii	7.2	0.3/0.2	19.9	19.6	19	14.1	15.6
3. P. dhofarensis	2.2	6.6	0.1/0.4	14.9	12.9	17.8	18.2
4. P. jensvindumi	4.4	7.8	4.3	0/0.2	15.7	18.2	19.4
5. P. neumanni	2.3	7.1	2.2	4.1	1.1/0.8	16.5	18.6
6. P. sinaitus	7.5	5.2	6.9	7.5	7.6	0.1/0.1	15.2
7. Pseudotrapelus sp.	8.7	6.7	8.1	9	8.6	7.2	0.6/0.3

Morphological characteristics of the specimens include morphometric and mensural characters. Measurements were taken with a digital calliper to the nearest 0.1 mm and meristic traits from the right side. We measured the following morphometric characters: snout-vent length (SVL), distance from the tip of the snout to the cloaca; tail length (TL), measured from the posterior lip of the cloaca to the tip of the tail; head width (HW), measured at the point of greatest width; head height (HH), measured at the point of greatest height; head length (HL), measured from behind the tip of the retroarticular process to the tip of the snout. Mensural characters included: number of upper labial scales (ULS); number of lower labial scales (LLS); number of precloacal pores (PCP); number of subdigital lamellae under the 3rd and 4th toes (3L and 4L, respectively; not including claw scale). We also noted if the 3rd toe is equal or longer than the 4th toe.

Results

Our dataset comprised 94 Pseudotrapelus specimens sampled from localities across the distribution range of the genus in Arabia and North Africa, including type localities, and six specimens of Acanthocercus from Saudi Arabia, Yemen and Oman. The nuclear dataset included 17 haplotypes of c-mos and 43 haplotypes of MC1R.

The ML and BI phylogenetic analyses of the complete concatenated dataset resulted in identical topologies with high bootstrap support (ML) and posterior probabilities (BI) values (Fig. 3A). In both analyses, Pseudotrapelus, and each of the recognised species within it, were recovered as monophyletic, with a general phylogenetic structure of two major clades. The two Riyadh samples form a clade with P. chlodnickii and P. sinaitus, though the phylogenetic relationships within this clade are not supported in either the ML or BI analyses. The second major clade consists of the remaining species, with P. jensvindumi sister to a subclade comprising the other three species with high support. In this subclade P. dhofarensis is sister to the west Arabian group of P. aqabensis and P. neumanni, though this relationship is not supported in the ML analysis. The haplotype nuclear networks inferred for the phased c-mos and MC1R markers show a similar pattern of no allele sharing between the new Riyadh samples and the other species of the genus (Fig. 3B). In the MC1R network there was no allele sharing between any of the species, contrasting the c-mos network in which allele sharing was present among the phylogenetically close Arabian species P. aqabensis, P. dhofarensis, P. jensvindumi, and P. neumanni. This pattern of allele sharing in the c-mos marker is probably due to incomplete lineage sorting rather than gene flow among the species. The mitochondrial genetic distances among Pseudotrapelus species (Table 1), including the two samples from Riyadh, were relatively high, ranging in 16S between 2.2 and 9% (with most values higher than 5%) and in ND4 between 11.2 and 19.9%; intraspecific distances ranged between 0 and 1.1% and 0.3–7.4%, respectively. The genetic distances between the Riyadh samples and the other recognised species were among the highest in the genus, ranging in 16S between 6.7 and 9% and in ND4 between 15.2 and 19.4%; between the two Riyadh samples the genetic distance was 0.6% in 16S and 0.3% in ND4.

The delimited mitochondrial and nuclear lineages within Pseudotrapelus are shown in Fig. 3A. The GMYC species delimitation analysis resulted in 11 mitochondrial distinct entities within the genus (numbered I–XI). The BPP analyses, based on the nuclear dataset assigned to the 11 mitochondrial GMYC lineages, supported seven distinct lineages within the genus, identical to the six recognised species and the lineage of the Riyadh specimens (numbered A–G). In all BPP analyses carried out, the Riyadh specimens were distinct from the other species. The topologies of the BPP species trees, however, changed with the shifting phylogenetic positions within the clades as described above for the phylogenetic trees based on the concatenated dataset.

The time-calibrated species tree of Pseudotrapelus is presented in Fig. 4. The topology was generally similar to that of the complete concatenated dataset, apart from within the clade of the Riyadh individuals with P. chlodnickii and P. sinaitus, and within the subclade of P. dhofarensis, P. aqabensis, and P. neumanni. Cladogenesis within Pseudotrapelus began in the Late Miocene around 9 Ma (95% Highest Posterior Density [HPD]: 6.2–12.1) into the two major clades. During the Late Miocene ca. 6.4 Ma (95% HPD: 4.4–9) the clade comprising the two Riyadh samples with P. chlodnickii and P. sinaitus diverged. Divergence within the other clade began with the split of P. jensvindumi from the ancestral population approximately 5.5 Ma (95% HPD: 3.7–7.8). At 4.2 Ma (95% HPD: 2.8–5.9) divergence occurred within the subclade of P. dhofarensis, P. aqabensis, and P. neumanni.

The morphological characters of the two Riyadh specimens were similar. For specimen IBECN6252 (adult male): SVL 73.5 mm; TL 140 mm (original); HH 10.5 mm; HW 16.1 mm; HL 22.9 mm; ULS 17; LLS 16; PCP 6 in a straight undivided line; 3L 21; 4L 12; 3rd toe longer than the 4th toe. For specimen IBECN13348 (adult female): SVL 74.9 mm; TL 85 mm (tail cut); HH 12.3 mm; HW 16.9 mm; HL 22.5 mm; ULS 14; LLS 14; PCP 4 in a straight undivided line; 3L 21; 4L 12; 3rd toe longer than the 4th toe.

Discussion

In this study we aimed to investigate the phylogenetic affinity of two specimens of Pseudotrapelus from the vicinity of Riyadh in the interior of Saudi Arabia. We applied multilocus phylogenetic and coalescent based methods to all recognised Pseudotrapelus species, and the two unknown samples from Riyadh, to generate concatenated and species trees and mitochondrial- and nuclear-based species delimitation analyses. The inferred topologies in our study, based on both the concatenated dataset and the species trees, were mostly congruent across analyses and generally support the current taxonomy (Tamar et al., 2016a). The Riyadh individuals are phylogenetically close relatives of two recognised species, P. sinaitus and P. chlodnickii, forming together a well-supported clade. Although relationships within this clade are not supported and change between the different analyses (those of the concatenated dataset and the species trees), their phylogenetic affinity is strong. Additionally, the genetic distances of the two Riyadh specimens from the other known species is among the highest in the genus (16S, 6.7–9%; ND4, 15.2–19.4%; Table 1). Moreover, the lack of shared haplotypes between the two Riyadh individuals and the remaining six recognised species of the genus in the two nuclear markers analysed further supports the notion that the Riyadh lineage has been genetically isolated for a long time without gene flow. The species delimitation analyses based on the mitochondrial and nuclear data also recover these Riyadh specimens as distinct from the other species.

According to our molecular results, the two Pseudotrapelus specimens collected from Riyadh in Saudi Arabia, an area that was not sampled in previous phylogenetic studies and more than 700 km away from their closest phylogenetic relatives sampled, P. aqabensis, represent a separate, deep, and distinct genetic lineage. These results suggest that the two Riyadh samples probably belong to a cryptic species within Pseudotrapelus. In view of the new data presented here, we believe that in order to perform an integrative and comprehensive taxonomic assessment of this species and genus, it will be necessary to collect additional morphological and ecological data across the rocky areas surrounding Riyadh and central Saudi Arabia.

The estimated divergence time of the clade comprised of the Riyadh individuals, P. sinaitus, and P. chlodnickii is during the Late Miocene ca. 6 Ma. We hypothesize that the aridification and fluctuating climate that prevailed during the Middle Miocene onwards and the progression of sandy areas in the Arabian region during the Late Miocene (Edgell, 2006, Preusser, 2009) have created distributional restrictions for the ancestral Pseudotrapelus populations, limiting them to rocky habitats. These range constraints and habitat fragmentation processes may have promoted vicariance and isolation within montane or hard-substrate taxa such as Pseudotrapelus, especially in the interior of Saudi Arabia. Similar patterns were also suggested for Uromastyx agamids (Tamar et al., 2018), Ptyodactylus geckos (Metallinou et al., 2015), and snakes of the genus Echis (Pook et al., 2009).

The reptile fauna of Saudi Arabia currently comprises over 100 species (Roll et al., 2017, Uetz et al., 2019). These species occupy a diverse array of habitats, ranging from sandy areas to gravel and rocks zones, and from low valleys to high mountain ridges. Different regions in Saudi Arabia have been surveyed in past years, providing valuable information regarding the occurrence and distribution of species (e.g., Farag and Banaja, 1980, Al-Sadoon, 1988, Al-Sadoon, 2010, Al-Shammari, 2012, Al-Sadoon et al., 2016, Al-Sadoon et al., 2017 and see reference therein). However, genetic studies of Saudi Arabian reptiles, or those incorporating samples from Saudi Arabia, are relatively few, despite numerous studies of reptiles elsewhere in the Arabian Peninsula. Within the general region of the Arabia, broad molecular studies have contributed greatly to our understanding of the evolution, systematics, and distributions of the local reptilian fauna (e.g., Portik and Papenfuss, 2012, Metallinou et al., 2012, Metallinou et al., 2015, Šmíd et al., 2013, Tamar et al., 2016a, Carranza et al., 2018, Burriel-Carranza et al., 2019). Integrative studies carried out in other areas of Arabia have uncovered considerable levels of undescribed diversity (e.g., Carranza and Arnold, 2012, Metallinou and Carranza, 2013, Vasconcelos and Carranza, 2014, Carranza et al., 2016, Šmíd et al., 2015, Šmíd et al., 2017, Sindaco et al., 2018, Tamar et al., 2019a), including several remarkable examples of cryptic diversity (e.g., Badiane et al., 2014, Garcia-Porta et al., 2017, Simó-Riudalbas et al., 2017, Simó-Riudalbas et al., 2018, Tamar et al., 2019b, Tamar et al., 2019c). In the case of Oman, these studies have resulted in a dramatic increase of 17.8% in the total number of reptile species in the past 13 years (Carranza et al., 2018).

The availability of DNA data and the incorporation of molecular methods are a valuable tool to understand phylogenetic relationships, assess distributions, evaluate conservation efforts, identify species, and as presented in this study, to discover cryptic diversity. However, the lack of specimens from Saudi Arabia incorporated in phylogenetic or phylogeographic studies, as mentioned above, makes the identity and the phylogenetic relationships of the local reptile species uncertain. This scarcity of data demonstrates the importance and the potential for better integrated systematic and taxonomic studies of Saudi Arabian species, which still may include undescribed taxa, as that within Pseudotrapelus uncovered in this study. We therefore advocate for further herpetological surveys in Saudi Arabia and the availability of material for future systematic studies.