New records are based on both systematic and opportunistic field surveys conducted by the authors in northwestern Saudi Arabia between 2021 and 2023. An additional distribution record of Tropiocolotes wolfgangboehmei comes from fieldwork carried out by MA in central Saudi Arabia in 2010. Individuals of Tropiocolotes were observed or captured during nocturnal walking surveys. Captured individuals were photographed and released at the exact location of their capture, following the collection of a small tail tip sample for subsequent DNA extraction and laboratory processing. Tissue samples were stored in 96% ethanol. In total, we collected 11 Tropiocolotes tissue samples from Saudi Arabia for this study (Fig. 1).

Figure 1. 

Geographic distribution of the Tropiocolotes species in Saudi Arabia and the neighbouring territories. The coloured squares (blue = T. yomtovi, red = T. nattereri sensu stricto, green = the AlUla candidate species, and yellow = T. wolfgangboehmei) represent genetically analysed material, and the white circles represent observations and photographed records. The black arrow points to the location of the Saudi islands of Tiran and Sanafir mentioned in the text. The new record of T. wolfgangboehmei is highlighted by a dashed yellow buffer. Geographical sampling is overlaid by the majority rule consensus map of the two binarised distribution models of T. yomtovi. Dark blue indicates areas of suitability recovered by both models that differed in the thinning strategy; light blue indicates regions predicted as suitable by only one of the models (see Materials and Methods for more details). The dashed line indicates the spatial background used for developing the model. The specimen depicted is an unvouchered individual of T. yomtovi from approximately 30 km southeast of Haql, Tabuk Province, photographed by L. Pola.

Genomic DNA was extracted from the ethanol-preserved tissue samples using the DNeasy® Blood & Tissue Kit (Qiagen, Germany) following the manufacturer’s instructions. To assess the phylogenetic position of our new material and combine it with previously published data, we PCR-amplified and sequenced up to four genetic markers: two mitochondrial (mtDNA) gene fragments, the ribosomal 12S rRNA (12S) and the NADH dehydrogenase subunit 2 (ND2), and two nuclear (nDNA) gene fragments, the oocyte maturation factor MOS (c-mos) and the melanocortin 1 receptor (MC1R). The PCR products were visualised by electrophoresis, purified using EXOSAP-IT® PCR Product Cleanup Reagent (Thermo Fisher Scientific), and Sanger-sequenced in both directions at Macrogen Europe (Amsterdam, the Netherlands). Primers and PCR conditions used follow Machado et al. (2021).

We used Geneious R11 (Kearse et al. 2012) to inspect the raw sequence files, assemble contigs, generate consensus sequences, and concatenate alignments. Heterozygous positions in the nuclear markers were identified by the Heterozygote plugin, checked by eye, and coded according to the IUPAC ambiguity codes. In total, we generated 77 new sequences for this study. An additional 165 sequences, originating mainly from the studies by Metallinou et al. (2012), Machado et al. (2018, 2021), Ribeiro-Júnior et al. (2022b), Pinho et al. (2023), and Liz et al. (2025), were retrieved from GenBank. Trachydactylus hajarensis was used as an outgroup to root the tree.

Sequences of all genetic markers were aligned separately using the online version of MAFFT v7 (Katoh et al. 2019). For 12S, we applied the Q-INS-i strategy, which accounts for RNA secondary structure, while the sequences of the remaining genes were aligned using the default auto strategy. Sequences of protein-coding genes were translated into amino acids using the appropriate genetic codes (vertebrate mitochondrial for ND2 and standard for the nuclear markers), and no stop codons were detected, indicating that no pseudogenes were amplified.

All samples used in this study, including information on their country of origin, GPS coordinates (datum WGS84), GenBank accession numbers, and the original source, are listed in Suppl. material 1: table S1.

The final concatenated dataset, comprising two mtDNA and two nDNA markers, contained 75 terminals, including the outgroup species. The total length was 2,546 base pairs (bp) – 433 bp of 12S, 1,051 bp of ND2, 394 bp of c-mos, and 668 bp of MC1R.

We performed Maximum Likelihood (ML) and Bayesian inference (BI) analyses using the concatenated dataset of the four genetic markers. The ML analysis was carried out in IQ-TREE (Nguyen et al. 2015) using its online web service W-IQ-TREE (Trifinopoulos et al. 2016). The concatenated dataset was partitioned by gene, and the best substitution models were selected automatically for each gene by ModelFinder (Kalyaanamoorthy et al. 2017) as implemented in IQ-TREE. The best models of nucleotide substitution were TIM2+F+G4 for 12S, GTR+F+I+G4 for ND2, K2P+I for c-mos, and HKY+F+G4 for MC1R. Branch support was assessed by the ultrafast bootstrap approximation algorithm (UFBoot; Minh et al. 2013) and the Shimodaira–Hasegawa-like approximate likelihood ratio test (SH-aLRT; Guindon et al. 2010), both with 1,000 replicates.

The BI was carried out using MrBayes v3.2.1 (Ronquist et al. 2012) with the same partitioning strategy as for the ML analysis. The best-fit substitution models were selected using PartitionFinder v1.1.1 (Lanfear et al. 2012) with the following parameters: branch lengths linked; MrBayes models; AICc model selection; user schemes search; and each gene representing a separate partition. The best models of nucleotide substitution were GTR+G for 12S, GTR for ND2, HKY+I for c-mos, and HKY+G for MC1R. The BI analysis was set to run in three independent runs, each for 20 million generations with parameters and trees sampled every 10,000 generations. All parameters (statefreq, revmat, shape, and pinvar) were unlinked for the partitions. Stationarity was determined by examining the standard deviation of the split frequencies between the three runs, the Potential Scale Reduction Factor diagnostic, and by confirming that all of the parameters had sufficient effective sample sizes using Tracer v1.7.2 (Rambaut et al. 2018). After confirming convergence of the runs, 25% of the posterior trees from each run were discarded as burn-in. A 50% majority-rule consensus tree was then produced from all post-burn-in posterior trees. Nodes with UFBoot ≥ 95, SH-aLRT ≥ 80, and Bayesian posterior probability (pp) ≥ 0.95 were considered well supported.

To inspect the level of nuclear allele sharing within the Tropiocolotes nattereri group, we reconstructed haplotype networks of the two nuclear loci. To resolve the heterozygous single-nucleotide polymorphisms, each alignment was phased using the PHASE algorithm (Stephens et al. 2001) as implemented in Hapsolutely v0.2.3 (Vences et al. 2024). Prior to phasing, we excluded several shorter sequences to avoid potentially misleading results. Both allele and phase thresholds were set to 0.5. Haplotype networks were then inferred from the phased alignments using the TCS algorithm (Templeton et al. 1992).

Uncorrected inter- and intraspecific p-distances with pairwise deletion were calculated for the mitochondrial markers 12S and ND2 in MEGA v11 (Tamura et al. 2021).

We assembled a dataset of 224 records (see Suppl. material 1), of which 93 (37 were DNA barcoded) originated from previously published studies (Metallinou et al. 2012; Machado et al. 2018, 2021; Ribeiro-Júnior et al. 2022b; Liz et al. 2025) and 131 (11 were DNA barcoded) were represented by the new data from Saudi Arabia. To reduce possible model bias resulting from spatial clustering of distribution records, we applied spatial thinning using the ‘spThin’ R package (Aiello-Lammens et al. 2015). Two thinning strategies were implemented, using minimum separation distances of 25 km and 50 km, respectively, between any two records. Each thinning procedure was repeated ten times, yielding ten and nine distinct, randomly sampled datasets, consisting of 48 and 27 records, respectively. To define the spatial background, a 200 km buffer was created around each record of the original dataset, and then a convex hull was made enclosing all buffered records. We manually adjusted the spatial background not to extend further west beyond the Sinai Peninsula and further north of Israel.

We considered a set of 19 bioclimatic variables (CHELSA v2.1; Karger et al. 2017), elevation, and slope to predict the potential species range. The grain of the spatial analysis was 2.5 arc-minutes. We used ENMTools (Warren et al. 2010) to test for spatial correlation between all variables. Only the variables with Pearson’s r < 0.75 and biologically meaningful for the species were retained in the analysis. The final set of variables included elevation, slope, BIO1 (mean annual air temperature), BIO2 (mean diurnal air temperature range), BIO3 (isothermality), BIO13 (precipitation of the wettest month), BIO14 (precipitation of the driest month), and BIO15 (precipitation seasonality).

To develop the species distribution models and to assess the importance of each variable, we used the maximum entropy approach implemented in Maxent v3.4.1 (Phillips et al. 2006). For each of the ten and nine input datasets, ten model replicates were run using the subsample resampling method, with 10,000 background points and a maximum of 500 iterations. In each run, 25% of the records were randomly used as test samples. The final potential distribution models were averaged across the ten replicates and reclassified from continuous suitability values into binary presence–absence maps using the maximum training sensitivity plus specificity threshold. A majority-rule consensus map was then generated in ArcGIS Pro using the ‘Raster Calculator’ tool, based on the two binary models derived from different occurrence data thinning strategies. Model accuracies were determined by means of the area under the receiver-operating curve (AUC). Elevation for each record was extracted from the digital elevation model using the ‘Extract Multi Values to Points’ tool in ArcGIS Pro.

Both ML and BI analyses of the concatenated dataset yielded identical topologies with overall well-supported nodes (Fig. 2). The clade comprising T. nattereri sensu stricto, T. yomtovi, and the candidate species from AlUla County was recovered as monophyletic with high nodal support (SH-aLRT = 99.5/UFBoot = 100/pp = 1.0; support values are given in this order hereafter). Each of the species was also strongly supported (T. nattereri sensu stricto: 99.9/100/1.0; T. yomtovi: 99.9/100/1.0; the AlUla candidate species: 100/100/1.0). However, the relationships among the three lineages within the clade remained unresolved (58.7/73/0.93). Mean genetic divergences (uncorrected p-distances) between T. nattereri sensu stricto and the candidate species, and between T. yomtovi and the candidate species, are 7.1% and 7% in 12S and 12.2% and 13.4% in ND2, respectively. Genetic divergences between T. nattereri sensu stricto and T. yomtovi are 7.7% in 12S and 12.5% in ND2.

Figure 2. 

A. Maximum likelihood phylogenetic tree reconstructed from the concatenated dataset of 12S, ND2, c-mos, and MC1R (2,546 bp in total). The outgroup Trachydactylus hajarensis is not shown in the figure. Support values are indicated by the black circles at nodes. Samples for which new genetic data were generated are highlighted in bold; B. Haplotype networks of the two nDNA gene fragments (c-mos and MC1R). Circle size is proportional to the number of samples that share that given allele. Short transverse bars on the connecting lines indicate the number of mutational steps between alleles. Colourations correspond to Fig. 1.

All newly obtained samples from northwestern Tabuk Province (LP-prefixed) and the two southernmost ones from Medina Province (JIR1146 and JIR1147) are nested within T. yomtovi. Genetic divergences within T. yomtovi from Israel, Jordan, and Saudi Arabia range from 0.3–7% (mean = 3%) in 12S and 0.1–11.7% (mean = 5.7%) in ND2.

The Saudi Arabian endemic T. wolfgangboehmei was recovered as part of a well-supported clade (95.6/98/1.0) comprising T. bisharicus, T. somalicus, T. nubicus, T. steudneri, and T. naybandensis, the latter of which was recovered as its sister taxon, albeit with poor support (0/56/0.9).

Haplotype networks of the two nDNA markers reveal no allele sharing among T. nattereri sensu stricto, T. yomtovi, and the candidate species from AlUla County, each exhibiting private alleles (i.e., not shared with the other species).

The average test AUC values of the models ranged between 0.82 and 0.86 (mean = 0.84) for the dataset with the 25 km thinning strategy and 0.76 and 0.82 (mean = 0.78) for the dataset with the 50 km thinning strategy, indicating that the models could be classified as ‘good’ and ‘fair’ following standard criteria for distribution model evaluation (Araújo et al. 2005). For the dataset using the 25 km thinning strategy, the first four variables that best predicted the distribution of T. yomtovi were slope (contribution 28.3–40.7%; mean = 31.9%), BIO13 (contribution 16.5–23.1%; mean = 20.0%), elevation (contribution 10.9–18.6%; mean = 14.6%), and BIO14 (contribution 12.8–15.8%; mean = 13.8%). For the dataset with the 50 km thinning strategy, the top four contributing variables were slope (contribution 23.8–39.7%; mean = 32.3%), elevation (contribution 13.0–28.4%; mean = 19.1%), BIO2 (contribution 13.0–20.7%; mean = 16.0%), and BIO14 (contribution 11.5–17.1%; mean = 14.6%).

According to the consensus predictive model (Fig. 1), the habitat suitable for T. yomtovi extends from the Sinai Peninsula in eastern Egypt, the Negev Desert in Israel, the Dead Sea Fault and its associated wadis in Israel and the Palestinian West Bank, western Jordan, and eastern Israel, northward to the Sea of Galilee. Further southeast, the model identified suitable areas from the coastal zones to the Midian and Hejaz Mountains, across the highland volcanic fields (harrats) of Tabuk and northwestern Medina Provinces, with marginal extension into western Jouf and Ha’il Provinces. The predicted distribution is largely continuous across this range, with the notable exception of a gap in habitat suitability throughout Mecca and Bahah Provinces. South of this interruption, suitable conditions reappear and continue in the high elevations of Asir and Najran Provinces into northwestern Yemen. When the two southernmost, non-genotyped records were excluded from the analysis, the predicted southern extent of the suitable niche was minimal (Suppl. material 1: fig. S1).

According to Machado et al. (2021), who time-calibrated the Tropiocolotes diversification, the split between T. nattereri sensu stricto and T. yomtovi dates back to the Late Miocene (8.5 million years ago [Mya], 95% highest posterior density interval: 4.28–10.09 Mya). Although our analyses did not resolve the relationships within this group, the combined evidence of high mtDNA genetic divergences between T. nattereri sensu stricto, T. yomtovi, and the candidate species (sensu Liz et al. 2025), which are comparable to those observed between the two described species (Ribeiro-Júnior et al. 2022b), together with complete allele segregation at the two analysed nDNA markers, suggests that the third lineage may warrant species-level status. Due to the unresolved relationships within the T. nattereri group and the lack of voucher specimens for morphological examination, the status of this lineage is left unresolved. Whether it truly warrants recognition as a distinct species remains to be evaluated through an integrative taxonomic approach. Before any taxonomic conclusions are drawn, the genomic structure and the level of gene flow between T. nattereri sensu stricto, T. yomtovi, and the candidate species should be assessed to avoid unnecessary taxonomic inflation (Burriel-Carranza et al. 2023).

The overall inferred phylogenetic relationships are congruent with previously published findings (Machado et al. 2021; Ribeiro-Júnior et al. 2022b). This study is the first to assess the phylogenetic position of the Saudi Arabian endemic T. wolfgangboehmei within the genus using expanded taxon sampling, following the initial generation of genetic data for the species (Šmíd et al. 2021).

Numerous field observations, supported by DNA-barcoded individuals, indicate that T. yomtovi is widespread throughout the northwestern part of Tabuk Province. Our findings, along with those of Liz et al. (2025), indicate significant range extensions further into Saudi Arabia. The results of the ecological niche modelling conducted here suggest that suitable habitats can also be found in the western parts of Jouf and Ha’il Provinces, where Tropiocolotes geckos have not been recorded. The species may extend as far south as Asir and Najran Provinces of Saudi Arabia and likely to northwestern Yemen. The occurrence of the genus further southeast appears indisputable; however, the existing published records of T. tripolitanus and T. steudneri from southwestern Saudi Arabia (see Masood and Asiry 2012; Alqahtani 2018) almost certainly represent misidentifications. In contrast, a notable citizen science observation of T. aff. nattereri from Najran Province, reported on the iNaturalist platform by C. Emberson (https://www.inaturalist.org/observations/191034349), is worth highlighting. Currently, and pending verification of the species identification using genetic tools, this appears to be the only reliable record, potentially marking the southernmost occurrence of T. yomtovi.

At the western edge of the range of T. yomtovi, the predictive modelling suggests that nearly the entire Sinai Peninsula is suitable for its occurrence. However, its distribution extent, as well as that of T. nattereri sensu stricto (see Werner 1973; Baha El Din 2006), remains to be clarified through molecular data and additional sampling. While in some parts of Israel and the Palestinian West Bank the predictive modelling appears to accurately reflect the known distribution of the species, in the Jordan Valley the results seem to overpredict the presence in areas where no Tropiocolotes species has been recorded (Werner 2016; Bar et al. 2021). Similarly, numerous wadi systems in eastern Jordan are predicted as suitable, despite the species likely being restricted to the lower part of the Rift Valley, including the Dead Sea basin, Wadi Araba, and its tributaries (Disi et al. 2001; Handal et al. 2016).

Our field observations on the species’ natural history are consistent with previous findings (Baha El Din 2006; Bar et al. 2021; Ribeiro-Júnior et al. 2022b). Tropiocolotes yomtovi is a nocturnal and ground-dwelling species inhabiting a wide range of desert environments, with a particular preference for rocky and stony terrain (Fig. 3). In terms of elevation, the species has been recorded from –380 m below sea level at the Dead Sea basin to around 2,000 m above sea level. The highest recorded individuals, at 1,983 m asl, come from the foothills of the Jebel al Lawz range, the highest mountain in Tabuk Province.

Figure 3. 

Habitats and unvouchered individuals of Tropiocolotes yomtovi from Saudi Arabia. A. Near Haql, NEOM region, Tabuk Province at 29.094°N, 34.949°E (photo: NR); B. Subadult female (sample code LP105) from approx. 30 km south-southeast of Haql (photo: LP); C. Nearby Al Sourah at 27.965°N, 35.483°E (photo: NR); D. Adult female from nearby Haql (photo: LP); E. Foothills of Jebel al Lawz at 28.735°N, 35.394°E (photo: LP); F. Adult unsexed from nearby Al Bad at 28.47°N, 34.985°E (photo: EF); G. Volcanic fields of Harrat ar Raha at 27.911°N, 36.439°E (photo: LP); H. Adult unsexed from Bajdah at 28.188°N, 36.037°E (photo: EF).

As the global conservation status of T. yomtovi has not yet been evaluated, we recommend that any future assessment consider the species’ relatively broad distribution and apparent abundance of individuals in suitable habitats. Given these factors and the absence of any known imminent threats, it is likely that the species would not fall into any of the threatened categories as defined by the IUCN.

Until now, T. wolfgangboehmei has been known only from the type series and a few additional records from areas north of Riyadh, from the King Khalid Royal Reserve (formerly Ath-Thumama) and the King Abdulaziz Royal Reserve (Wilms et al. 2010; Šmíd et al. 2021; van Huyssteen et al. 2024). Herein, we present a new record (23.3281°N; 46.6836°E) originating from the Ibex Reserve, approximately 200 km south of the type locality (Fig. 1). The new locality reported here extends the known range of T. wolfgangboehmei further south along the Tuwaiq Escarpment. Likely, the species is more broadly distributed along this mountain ridge, which provides a suitable habitat for this rupicolous ground gecko, both to the north and to the south. With our current knowledge, the range of T. wolfgangboehmei resembles that of the recently described agamid Pseudotrapelus tuwaiqensis, another poorly known endemic species from central Saudi Arabia (Tamar et al. 2023), and follows the central Arabian part of the range of the lacertid lizard Mesalina cryptica, which also inhabits the rocky foothills of the escarpment (Šmíd et al. 2025). Its secretive lifestyle and small size may explain why it has not been encountered more frequently in field surveys around Riyadh.