Some like it hot: Past and present phylogeography of a desert dwelling gecko across the Arabian Peninsula

Deserts represent dynamic ecosystems that support communities of endemic and specialised species. We analysed the role of present and past climatic conditions in shaping the distribution of the widespread Bunopus geckos in the Arabian and south‐west Asian deserts. We studied their phylogeographic and demographic history to test whether the Bunopus geckos colonised Arabia from Asia or, vice versa, Asia from Arabia and to identify migration corridors that have historically enabled the dispersal of Bunopus geckos.


| INTRODUC TI ON
Past climatic oscillations have played a crucial role in shaping the current distribution of biodiversity worldwide.In contrast to the temperate biomes that were characterised by recurrent warm-to-cold climate shifts during the Quaternary, subtropical deserts have experienced alternations of humid and arid conditions (Glennie, 2020;Hesse et al., 2004).Arid deserts have received considerably less attention from biodiversity researchers compared to other ecosystems (Durant et al., 2012), and our understanding of the effects of past climatic fluctuations on their biota remains rather incomplete (Douglas et al., 2006;Pepper & Keogh, 2021).
The Arabian Peninsula is an isolated subcontinent that was historically part of Africa, from which it drifted away to the northeast after their split in the mid-Oligocene to Early Miocene (Bosworth et al., 2005).The peninsula is rimmed by mountains that run along the seas which flank Arabia from the west, south and east.The Arabian interior is dominated by basalt flows and salty plains, but most notably by sand and gravel deserts with the Rub' Al Khali sand sea (also called the Empty Quarter) being the dominant feature (Edgell, 2006).
The extent of these deserts was, however, not stable throughout the history of Arabia as it responded to fluctuating climatic conditions.The climate of Arabia is believed to have been humid with well-developed systems of seasonal river valleys (termed wadis) that drained the peninsula during the Late Pleistocene and Early Pliocene (Anton, 1984;Dabbagh et al., 2020).At that time, these deserts were covered by open savanna woodlands and the giant sand dunes in southern Arabia were interspaced with lakes and swamps (Edgell, 2006;McClure, 1976;Vincent, 2008).During the Quaternary, the climate of Arabia fluctuated regularly between hyper-arid (similar to those of today) and humid that was characterised by increased precipitation and the reactivation of river and lake systems in the interior (Breeze et al., 2015;Dinies et al., 2015).This generated a complex spatial and temporal mosaic of habitats that likely impacted the population dynamics of the desert dwelling biota and provided windows of opportunity for dispersal for animals and hominins (Parker, 2010;Stimpson et al., 2016).
The mountains of Arabia have been shown to support unique diversity of squamates with exceptional levels of endemism and as such have received considerable scientific attention (Carranza et al., 2016;Carranza & Arnold, 2012;Garcia-Porta et al., 2017;Metallinou et al., 2015;Šmíd et al., 2013, 2017).By contrast, the fauna of the inland deserts has been overlooked until relatively recently (Metallinou et al., 2012;Pola et al., 2021;Šmíd et al., 2021).The phylogeographic histories of widespread fauna may reveal how past climatic oscillations have affected the whole peninsular faunal assemblage.
Geckos of the genus Bunopus Blanford, 1874 are habitat generalists that inhabit a broad range of habitats throughout Arabia and the Iranian Plateau, ranging from southern Israel in the west to central Pakistan in the east (Sindaco & Jeremčenko, 2008;Šmíd et al., 2014, 2021).Their ground-dwelling habits in combination with their widespread distribution and high local population densities make them a suitable model for studying present and past dispersal dynamics in these hyper-arid environments.Two to three species of Bunopus are currently recognised, B. crassicauda, B. tuberculatus and presumably also B. blanfordii.While the first is endemic to Iran, the second occupies the rest of the genus' range including the entire Arabian Peninsula.The status of the third species, B. blanfordii, remains questionable, and the species is often considered a synonym of B. tuberculatus (see Bauer et al., 2013).The apparent uniformity of B. tuberculatus across its range was first doubted by Červenka et al. (2008) and later on by Khosravani et al. (2017), who found cryptic diversity within the species indicating that it likely represents a species complex.Nonetheless, the status of the populations occurring in the Arabian Peninsula has not been thoroughly investigated.
In this study, we analyse the role of present and past climatic conditions in shaping the distribution of a widespread generalist species in the Arabian deserts.We analyse sequence data of two mitochondrial and two nuclear markers to untangle the phylogenetic and phylogeographic history of the Bunopus species and populations across the entire distribution of the genus range.We use the genetic data and a dense sampling from throughout the Arabian Peninsula to infer the demographic history of the Arabian populations since the Pliocene to the present.Finally, we apply species distribution modelling to identify the extent of suitable habitats for Bunopus in Arabia in the present and in the past.The integration of the genetic and spatial results allows us to analyse the connectivity of landscapes across the peninsula and its role in the migration of this broad-ranging genus.This ultimately leads to the identification of dispersal corridors that enable migration and promote gene flow among Bunopus populations in Arabia.(DJ).We assembled a total of 88 samples covering densely the Arabian part of the genus range.We retrieved sequences for 93 additional samples from the GenBank and BOLD (www.bolds ystems.com) databases.The final dataset included 174 ingroup samples from across the entire range of the genus (Figure 1).We adopted the code system proposed by Khosravani et al. (2017) for the undescribed candidate species (Bunopus sp.1-5).

| Taxon sampling and outgroup selection
As for the outgroup taxa, there are disputes with regards to what gecko genus is the closest relative to Bunopus.Numerous molecular phylogenetic studies showed a sister relationship between Crossobamon and Bunopus (Bauer et al., 2013;de Pous et al., 2016;Gamble et al., 2012;Machado et al., 2019Machado et al., , 2021;;Metallinou et al., 2012).However, other studies recovered Crossobamon to be nested within Bunopus (Agarwal et al., 2014;Pyron et al., 2013;Zheng & Wiens, 2016), making the latter paraphyletic.The phylogenetic position of the two genera with respect to each other remains disputed and using only Crossobamon as the outgroup for the phylogenetic analyses might affect the results.We therefore used samples of both known Crossobamon species, C. eversmanni and C. orientalis, that cover broadly their ranges, but we also included more distant taxa Agamura persica and Trachydactylus spatalurus to root the tree.

| DNA extraction, amplification and sequence analysis
Genomic DNA was extracted from ethanol-preserved tissue samples using Tissue Genomic DNA Mini Kit (Geneaid) following the manufacturer's instructions.We PCR-amplified up to four genetic markers: two mitochondrial (mtDNA): the 12S rRNA (12S) and the cytochrome c oxidase subunit 1 (COI), and two nuclear (nDNA): the recombination activating gene 2 (RAG2) and the oocyte maturation factor MOS (c-mos).The PCR products were purified using EXOSAP-IT® PCR Product Cleanup Reagent (Thermo Fisher Scientific) and were Sanger-sequenced in both directions in Macrogen Europe (Amsterdam, the Netherlands).Primers, their sequences and PCR conditions are provided in Table S1.
Raw sequence data were inspected and contigs assembled using Geneious R11 (Kearse et al., 2012).Heterozygous positions in the nuclear markers were identified by the Heterozygote Plugin and were coded according to the IUPAC ambiguity codes.Sequences of each genetic marker were aligned independently by MAFFT (Katoh et al., 2019) using the default auto strategy for all genes except the 12S, where the Q-INS-i strategy that considers the secondary structure of RNA was applied.For the 12S alignment, we used Gblocks (Castresana, 2000) to trim poorly aligned regions with gaps.Sequences of protein-coding genes were translated into amino acids and no stop codons were detected.Samples used in this study are listed in Table S2.

| Phylogenetic and nuclear network analyses
We performed maximum likelihood (ML) and Bayesian inference (BI) analyses using the concatenated dataset of the four markers.The ML was carried out in IQ-TREE (Nguyen et al., 2015) using its online web interface W-IQ-TREE (Trifinopoulos et al., 2016).The dataset was partitioned by gene with models selected automatically by ModelFinder (Kalyaanamoorthy et al., 2017) as implemented in IQ-TREE.Branch support was assessed by the Shimodaira-Hasegawalike approximate likelihood ratio test (SH-aLRT; Guindon et al., 2010) and the Ultrafast bootstrap approximation algorithm (UFBoot; Minh et al., 2013), both with 1000 replicates.
The BI was performed using BEAST v.2.5.2 (Bouckaert et al., 2019).The dataset was again partitioned by gene, site and clock models were unlinked across partitions.We applied the reversible-jump based method for best model selection with four gamma-distributed rate categories (Bouckaert et al., 2013).The relaxed lognormal clock model was applied to each partition.We used the coalescent constant population tree prior with a 1/X population size prior.Lognormal prior distributions were selected for the clock parameter priors (ucldMean), with the mean = 0.1 and standard deviation = 1.25.Rate variation across lineages (ucldStdev) of each partition was estimated using an exponential prior distribution (mean = 0.5).The analysis ran three times for 5 × 10 7 generations through the CIPRES Science Gateway (Miller et al., 2010) with trees and parameters sampled every 2 × 10 4 generations.Tracer v.1.7.1 (Rambaut et al., 2018) was used to check the effective sample size of all parameters and to ensure that stationarity and convergence had been reached.Tree files were then combined using LogCombiner after discarding 10% of the posterior trees as burn-in.The maximum clade credibility tree was identified using TreeAnnotator.Nodes that received SH-aLRT ≥80, UFBoot ≥95 in the ML analysis, and Bayesian posterior probability (pp) ≥ 0.95 were considered strongly supported.
Inter-and intraspecific relationships were inspected by reconstructing haplotype networks of the nuclear loci.To resolve the heterozygous single nucleotide polymorphisms, the alignments of RAG2 and c-mos were phased using the PHASE algorithm (Stephens et al., 2001) as implemented in DnaSP v.6 (Rozas et al., 2017) with probability threshold set to 0.7.Prior to phasing, we excluded several shorter sequences and the outgroups to avoid misleading results.Haplotype networks were constructed from the phased alignments using the TCS algorithm (Clement et al., 2000;Templeton et al., 1992) implemented in PopART (Leigh & Bryant, 2015).

| Estimation of divergence times
We calibrated the phylogeny with the substitution rate of the 12S gene estimated by Carranza and Arnold (2012), with the mean clock rate of 0.00755 and standard deviation of 0.00247.Similar approach has proven useful when calibrating trees of other gekkonid taxa in the region (Carranza & Arnold, 2012;de Pous et al., 2016;Machado et al., 2021).The analysis was run in BEAST through CIPRES with parameters and priors as described above.The only difference was that we applied the Yule tree prior that assumes a constant lineage birth rate with sampling limited to one sample per lineage.The described and candidate species of Bunopus and all the outgroup species were thus represented by one sample each.Only the Arabian candidate species, Bunopus sp. 4, was represented by four samples to be able to estimate divergence times between the major geographic lineages as recovered by the ML and BI analyses (see Results below).The analysis ran three times for 3 × 10 7 generations and was sampled every 3000 generations.

| Ancestral area reconstruction
To infer the biogeographical history and ancestral ranges of Bunopus, we used the R package 'BioGeoBEARS' (Matzke, 2013).
We used the calibrated tree as input and pruned the outgroup species prior to the analysis.We also retained only one tip for the Bunopus sp. 4 candidate species.We defined three biogeographic areas based on the geological history of the region (Popov et al., 2004): (i) Arabia, for the Arabian Peninsula including the desert in southern Jordan; (ii) Mesopotamia, for the lowlands along the Euphrates and Tigris Rivers; (iii) mainland Asia east of the Zagros Mountains in Iran.We assigned each tip to one or more of these areas based on the current distribution of that lineage.
We restricted the maximum number of areas in which ancestral nodes could occur to two and performed ancestral reconstructions using the three models available in BioGeoBEARS: Dispersal-Extinction-Cladogenesis (DEC; Ree & Smith, 2008), DIVALIKE (Ronquist, 1997) and BAYAREA (Landis et al., 2013).We examined the plausibility of the results of each model empirically and we also assessed the fit of the models by the Akaike information criterion corrected for sample size (AICc; Akaike, 1973).The + J parameter that is implemented in BioGeoBEARS and allows including founder-event speciation (jump-dispersal; Matzke, 2014;Ree & Sanmartín, 2018) was not included in the models.

| Inferring the demographic history
To estimate population size changes through time, we used Extended Bayesian Skyline Plots (EBSP) using BEAST v.2.5.2 (Heled & Drummond, 2008).Since the focus of the study lies on the Arabian populations of Bunopus, we pruned the dataset for this analysis to only include samples of the candidate species from Arabia, Bunopus sp. 4, of which there were 83.The BEAST settings followed those described above for the BI analysis.The average number of population changes was modelled with a Poisson prior distribution.The analysis ran three times for 2 × 10 8 generations, and 10% of the posterior parameter values were discarded as burn-in.

| Modelling potential distribution in the present and in the past
We compiled a database of available distribution records by searching published literature, museum catalogues, public databases (e.g., GBIF) and gathering field observations.In total, we assembled 1314 records of Bunopus representing 920 unique localities.We thinned the dataset using the 'spThin R' package (Aiello- Nineteen bioclimatic variables were downloaded from CHELSA (Karger et al., 2017) at the resolution of 2.5 arc-minutes and cropped to the study area.BIO8, BIO9 and BIO18 were excluded because they showed spatial artefacts and BIO14 because it showed no variation across the study area.In addition to the bioclimatic layers, we used layers for elevation and slope.To be able to project habitat suitability in the past when sea level was different from today, we created a layer of elevation that also contained negative values for areas below the sea level (bathymetry data downloaded from GEBCO; https:// www.gebco.net).We tested for collinearity between the layers using ENMTools (Warren et al., 2010) and of those with correlation over 0.75 we retained only the more biologically meaningful ones (Elith & Leathwick, 2009).The final set contained these variables: elevation, slope, mean diurnal air temperature range (BIO2), temperature seasonality (BIO4), mean daily mean air temperatures of the coldest quarter (BIO11), precipitation seasonality (BIO15), mean monthly precipitation amount of the wettest quarter (BIO16), and mean monthly precipitation amount of the driest quarter (BIO17).
We used Maxent 3.3 (Phillips et al., 2006) to develop the species distribution model and to assess the importance of each variable.Ten model replicates with the cross-validate resampling method were run for each of the ten input datasets, using 10,000 background sample points and with 5000 maximum iterations.The area under the curve (AUC) was assumed as a measure of individual model fit.
The final model of potential distribution was averaged over the ten replicates.To test whether the models performed better than random, we generated 100 null models, each for a set of 161 records randomly generated within the study area and with settings similar to the models based on real data.
To assess the dynamics and stability of the Bunopus distri- The elevation and slope layers were also included in the paleo projections.The elevation was manually adjusted for each of the past time periods to reflect the sea level difference at that time compared to the present.For projections to the mid-Pleistocene and mid-Pliocene, mean diurnal air temperature range was excluded as it was not available for those time periods.The paleo projections were run ten times each with the final model averaged over the ten runs.Input layer and parameter details of the distribution modelling are reported in an ODMAP protocol file (Zurell et al., 2020) in the Supplementary Material.

| Identifying contemporary dispersal corridors
We analysed contemporary spatial connectivity of the Bunopus pop- The percentage of least-cost path value was used to select the LCC with the high, mid and low cut-off values being respectively 5, 2 and 1.

| Spatial analysis of population structure
We assessed the genetic structure of the Arabian populations of Bunopus and identified spatial genetic neighbourhoods using the 'MEMGENE' R package (Galpern et al., 2014).MEMGENE regresses Moran's Eigenvector Maps (MEM), it is variables describing patterns of positive and negative spatial autocorrelation, against genetic distances to detect genetic structure and visualises spatial components of genetic dissimilarity among individuals (Galpern et al., 2014).
Based on the results of the phylogenetic analyses (see below) we included in this analysis only samples and localities of the candidate species Bunopus sp. 4. We calculated pairwise genetic distances between all samples on the ML tree using the Geneious software.We also visualised the correlation between the geographic and genetic distances.We used functions from the 'MASS' (Venables & Ripley, 2002) and 'adegenet' (Jombart & Ahmed, 2011) packages to create a kernel density plot of Bunopus records in Arabia to highlight regions of increased point density in the plot.We calculated linear geographical distances between sampled sites using the 'raster' R package (Hijmans et al., 2014) and correlated them with the genetic distances calculated above.In addition to the Euclidean (straight) distances, we also calculated least-cost path distances between all pairs of points using the corridor layer identified in section 2.8 as a cost layer and correlated this distance matrix with the genetic distances.We are aware that the correlation between genetic and geographic distances does not account for spatial autocorrelation, we however find it useful for visualising the relationships between the variables.

| RE SULTS
For this study, we generated 351 new sequences for 110 samples of the total of 202 samples used in the analyses.Sampling completeness (i.e., all four gene sequences available per sample) was 79.8% for the samples newly sequenced in this study and 55.6% with the GenBank and BOLD sequences included.

| Phylogenetic analyses
Both ML and BI analyses resulted in almost identical topologies in most nodes.According to the results (Figure 1 trees.In the ML tree, it was supported as sister to the whole Bunopus clade (98.9/100), in the BI tree it was sister to Agamura persica, although with low support.
In light of the paraphyly of the genus Crossobamon recovered in both ML and BI analyses, we run additional analysis to test whether the genus is significantly non-monophyletic.We constrained the topology of the tree and forced the two Crossobamon species to form a clade.We used the approximately unbiased (AU), the

| Estimation of divergence times
The initial split within the genus that separated the Iranian (including C. orientalis), and the Arabian clades was estimated to take place

| Ancestral area reconstruction
The results of the biogeographic reconstruction were largely congruent between the tested biogeographic models.DIVALIKE was the most plausible of the models (Table S3) and we therefore present only the results based on this model (Figure 3).

| Inferring the demographic history
The reconstruction of the demographic history of the Arabian populations shows a stable population trend since the split between Bunopus sp. 3 and sp. 4 at 3.4 Ma until about 200 ka (Figure 3).
At that time the population size started decreasing considerably, which continued until after the LIG (ca.130 ka).At about 80 ka, the trend turned, and the population increased almost to the predrop level.

| Present and past potential distribution
Mean AUC for the present ranged between 0.744 and 0.768, with the mean being 0.759.The consistency of the AUC values across the models along with extremely low standard deviation values of all runs (0.055-0.06; mean = 0.057) implies model stability regardless of the input data.The models performed significantly better than the null models (AUC: 0.586-0.688;mean = 0.637).The AUC values of the models would be categorised as 'fair' according to standard criteria for distribution model evaluation (Araújo et al., 2005).It should however be noted that it has been shown that predictive models of generalist species with broad environmental niches, such as Bunopus, achieve lower AUC values compared to habitat specialists (Connor et al., 2018).The most important environmental predictors were the elevation (contribution 48.6%-54.6%;mean = 50.8%),precipitation seasonality (contribution 13.9%-19.4%;mean = 17.6%), temperature seasonality (contribution 7.2%-10.0%;mean = 8.5%) and mean diurnal air temperature range (contribution 6.3%-9.8%;mean = 8.0%).
The predictive model based on the present environmental conditions showed that large parts of eastern Arabia, coastal western Arabia and coastal Iran support habitat that is suitable for Bunopus (Figure 4).LGM (ca.21 ka), the western Arabian coast was also suitable, but its extent has since been retreating.In the late Pleistocene (14.7-12.9ka) and early Holocene (12.9-11.7 ka), the range expanded and covered most of eastern Arabia, including the inland deserts which, however, became unsuitable again in the late Holocene (4.2-0.3 ka).

| Dispersal corridors
The dispersal corridors inferred for the three schemes showed congruent spatial patterns (Figure S3).Al Khali sands (Figure 5).
F I G U R E 3 (a) Time-calibrated tree and ancestral area reconstruction of Bunopus.All nodes were supported with posterior probabilities higher than ≥0.95.Mean age estimates for the branching events are provided below each node with the 95% HPD interval in brackets and also likely that at least some of the candidate species will warrant formal taxonomic recognition.
The three lineages of the Arabian clade show clear differentiation at the mitochondrial level but a certain overlap in the nuclear markers (Figures 1, 2).It may be a result of their relatively recent split that simply did not provide enough time for the slowly evolving nuclear genes to differentiate.Alternatively, it could be caused by events of introgression.This remains to be tested with a broader sampling of genomic loci (e.g., SNPs; work in progress).The broad distribution of Bunopus in Arabia provides suitable grounds for comparison with other widespread Arabian genera.Most previous studies of other pan-Arabian squamates have uncovered cryptic diversity present across the peninsula and concluded that the diversity of species is in fact much higher than had been previously thought.These findings often resulted in descriptions of new microendemic species, with geckos being the most taxonomically dynamic group of reptiles in this respect (e.g., Carranza et al., 2016;Machado et al., 2019;Simó-Riudalbas et al., 2017, 2018;Šmíd et al., 2013, 2015, 2017, 2023;Tamar et al., 2019;Vasconcelos & Carranza, 2014).Bunopus, however, shows a completely different pattern.The results of our phylogenetic analyses imply that, despite the broad distribution of Bunopus sp. 4 lineage, it harbours only low genetic diversity in Arabia.

| Biogeographic and demographic history
The biogeographic reconstructions together with the timecalibrated analysis indicate that Bunopus originated in mainland Asia in the mid-Miocene (Figure 3).All species of the Iranian clade are confined to the Iranian Plateau and adjoining parts of mainland Asia and do not seem to ever have expanded anywhere else.On the other hand, the Arabian clade was estimated to have dispersed from the Iranian Plateau to the Arabian Peninsula during the Pliocene/ Pleistocene.At that time, the two landmasses were connected by a continental land bridge (Popov et al., 2004), which likely facilitated biotic interchange between Arabia and mainland Asia (Badiane et al., 2014;Simó-Riudalbas et al., 2019;Tamar et al., 2018;Tamar et al., 2021).
The colonisation of the Arabian Peninsula was a very successful one indeed.The historical range reconstruction shows that Bunopus managed to disperse from Mesopotamia in the north throughout the Arabian Peninsula to its eastern-and southernmost margins (Figure 3).The facts that the distribution of Bunopus spans across the entire Arabia and that its genetic structure throughout the peninsula is rather shallow point to the present distribution of Bunopus sp. 4 being a result of rapid dispersal with ongoing gene flow.This is further supported by the spatial patterns of genetic data that shows genetic homogeneity across most of the Arabian Peninsula (Figure 6).
Interestingly, the process of range expansion does not seem to have been associated with expanding population size which was inferred here to have been stable since its split from the Mesopotamian lineage until about 250 ka when it dropped substantially (Figure 3).
After the decline, the population however returned rapidly to its original size.This drop may also be discernible in the projections of the predictive distribution model to past climatic conditions that indicate range contraction at the time of the Last Interglacial (130 ka).
At that time, Arabia underwent a predominantly moist climatic phase that was interwoven with short windows of semi-arid to arid conditions (Edgell, 2006;Vincent, 2008) which might have resulted in habitat fragmentation and subsequent population isolation.Since the Last Interglacial, however, the conditions became generally more arid again and Bunopus started repopulating Arabia.The considerable range expansion in the latest Pleistocene to early Holocene may be linked with the hyper-aridification of Arabia and the expansion of sand dunes at that time (Vincent, 2008).It is worth noting that the seabed of what is today the shallow Arabian Gulf presented suitable habitat for Bunopus during the glacial sea-level drops and likely formed a corridor for migration between the Iranian and Arabian populations (Lambeck, 1996).Taken together, the dynamic system of pulsating habitats that oscillated in response to the changing climatic conditions between humid and hyper-arid seems to have played a crucial role in shaping the present and past distribution of the desert adapted Bunopus geckos in Arabia.

| Dispersal across Arabia
The past distribution models imply that the range of Bunopus sp.et al., 2019;Sindaco et al., 2018).
It should be stressed that although the distribution models performed well for such a broadly distributed and generalist taxon, there were still regions where Bunopus was not predicted to occur despite the presence of records in these areas.For example, several distribution records are available from northern Saudi Arabia, but the region was not found suitable for the geckos.If these places are truly suboptimal for Bunopus and the existing distribution points represent sinking populations or if the predicted absence is caused by the scarcity of data is at the moment uncertain.The latter possibility seems very plausible.However, our field experience has taught us that Bunopus population densities vary considerably across Arabia and that while in some regions it is the most abundant reptile species (e.g., in central Saudi Arabia around the city of Riyadh), in other seemingly suitable desert habitats it is extremely rare (e.g., southwestern Arabia).Hence, until more field work is conducted in northern Arabia, we prefer not to draw conclusions on the predicted absence of Bunopus in these places.

| Taxonomic account
The results of all phylogenetic analyses conducted for this study support the paraphyly of Bunopus with Crossobamon orientalis being nested within Bunopus.Similar results have been confirmed in some previous studies, which were however always based on much sparser taxon sampling (Agarwal et al., 2014;Pyron et al., 2013;Zheng & Wiens, 2016).The genus Crossobamon currently contains two species -C.eversmanni and C. orientalis.
The phylogenetic position of the former species in the tree of (Blanford, 1876) that should be used from now on.A detailed list of chresonyms is available in the Supplementary Information.The cryptic diversity within the genus Bunopus (Červenka et al., 2008;Khosravani et al., 2017) as well as the status of the enigmatic B. blanfordii (Bauer et al., 2013) remain a task for future taxonomic investigation.Besides the polyphyly of Crossobamon found in our study we also noted deep genetic divergences within C. eversmanni throughout its range, suggesting possible cryptic diversity.

ACK N OWLED G EM ENTS
We are indebted to the National Center for Wildlife (NCW), Saudi Arabia, and its staff who participated on the field trips, their sup- The final concatenated alignment of the four markers was 1842 base pairs (bp) long -378 bp of 12S (after Gblocks trimming), 663 bp of COI, 408 bp of RAG2 and 393 bp of c-mos.F I G U R E 1 (a) Maximum likelihood phylogenetic tree reconstructed from the concatenated dataset of 12S, COI, RAG2 and c-mos genes (1842 bp).The tree was rooted using Trachydactylus spatalurus and Agamura persica (not shown in the figure).Support values (SH-aLRT/ UFBoot/pp) are indicated by the circles at nodes with colours explained in the legend under the tree.Colours of tree branches match those of the sampled sites in (b).(b) Map showing the geographical sampling across the Middle East.Complete trees with original ML and BI support values are provided as Figures S1, S2, respectively.Taxon names correspond to changes proposed in this study.Specimen depicted is an individual photographed in south Jordan (Photo: Lukáš Pola).
Lammens et al., 2015) to reduce possible model bias resulting from high concentrations of distribution records from thoroughly explored areas (e.g., the UAE, Oman;Carranza et al., 2018Carranza et al., , 2021;;Burriel-Carranza et al., 2019).We used a radius of a minimum of 50 km to separate any two records and run the thinning ten times, which produced ten different and randomly sampled datasets of 161 records.Since we focus in this study on the Arabian populations of Bunopus, we only used records of the clade containing Bunopus tuberculatus sensu stricto and the candidate species Bunopus sp. 3 and sp. 4. We pooled records of these three lineages together for the modelling purposes.The reasoning was that although they show a certain degree of genetic differentiation, it is mostly in the mitochondrial DNA and only in a limited way in the nuclear DNA, and with our current knowledge, it cannot be ruled out that the three lineages represent a single species.The species and candidate species from the Iranian Plateau were not included in the modelling since their environmental niches may differ from those of the Arabian clade and because they were not the primary aim of this study.
schemes as follows: (i) all sites of the three lineages -B.tuberculatus, Bunopus sp. 3, and Bunopus sp. 4 -were pooled together; (ii) samples were assigned to the three lineages, which were treated as distinct evolutionary entities; and (iii) the three lineages were treated as separate groups, and samples of Bunopus sp. 4 were further divided to five groups based on the intraspecific structure of the phylogeny.

Forward
selection of positive and negative MEM eigenvectors against genetic distance added eigenvectors to a regression model until they ceased to improve model fit.Principal component scores of the predicted values are defined as Memgene variables and we used those with the highest R 2 values to produce maps of the spatial patterns of genetic relationships.
Shimodaira-Hasegawa (SH)  and the Kishino-Hasegawa (KH) tests to compare this enforced topology with the unconstrained tree.Persite log likelihoods were calculated in raxmlGUI v.1.5(Silvestro & Michalak, 2012) and p-values were calculated using CONSEL v.0.1(Shimodaira & Hasegawa, 2001).The results indicate that the monophyly of Crossobamon can be significantly rejected (AU: 0.001; SH: 0.002; KH: 0.002).The haplotype networks (Figure2) show a certain degree of allele sharing between the Bunopus species, including Crossobamon orientalis, in both nuclear markers.The only species that have all alleles private (i.e., not shared with other species) are Crossobamon eversmanni, Bunopus crassicauda and Bunopus sp. 1.All the other species share alleles of one or both nuclear markers with some other species.Within the Arabian clade, B. tuberculatus sensu stricto possesses unique haplotypes in RAG2, and all the three species of that clade share one common allele in c-mos.
The suitable habitat covers most of Oman and the UAE except the Hajar and Dhofar Mountains and regions adjoining the Rub' al Khali Desert.It extends along the Arabian Gulf through Qatar and Kuwait to southwestern Iran and then further along the Gulf to south-eastern coastal Iran.There is a narrow band of suitable habitat along the Red Sea coast in northwestern Arabia.It connects to the eastern part of the suitable habitat through a longitudinal belt that crosses central Arabia.Interestingly, most of southern Arabia (Yemen and southwestern Saudi Arabia including the Rub' Al Khali Desert) and northern Arabia (the An-Nafud Desert) were not found to be suitable for Bunopus.Projections to past climatic conditions showed that eastern Arabia and most of the Arabian Gulf coast have harboured suitable habitat F I G U R E 2 Haplotype networks of the RAG2 and c-mos nuclear markers.Circle size is proportional to the number of samples that share that allele.Transverse bars on the connecting lines indicate the number of mutational steps between alleles.Colours correspond to those in Figure 1.Taxon names correspond to changes proposed in this study.(Figure 4).The extent of suitable habitat was very similar in mid-Pliocene (3.3 Ma) and mid-Pleistocene (ca.787 ka).It retreated during the LIG (ca.130 ka) and covered only central Oman and southern coasts of the Arabian Gulf.This habitat reduction was followed by a subsequent north-westerly expansion along and into the desiccated Arabian Gulf during the LGM (ca.21 ka).During the Most of lowland Oman and the coastal UAE are suitable for the dispersal of Bunopus, and the main migration corridor stretches from there along the Arabian Gulf coast through Qatar and Saudi Arabia to Kuwait, from where it continues across central Arabia in a broad, longitudinal belt all the way to the Red Sea coast.The northern part of the Saudi Red Sea coast from around the city of Jeddah to the border with Jordan also promotes Bunopus population connectivity.The isolated populations in Yemen and southern Saudi Arabia are connected by dispersal routes to the southern Oman and Arabian Gulf populations, some of which run across the Rub' indicated with the blue horizontal bars.The biogeographic areas defined for the analyses are in the inset map.Pie charts at the nodes show the probability of each ancestral area.(b) Extended Bayesian Skyline Plot showing the temporal dynamics of the effective population size of the Arabian populations (Bunopus sp.4) since its split from its sister lineage, Bunopus sp. 3. Taxon names correspond to changes proposed in this study.The first axis (MEMGENE1) showed a significant genetic structure that separates the northwestern Arabian Bunopus populations from the rest of the peninsula (Figure6).Curiously, a sample from Qatar showed genetic affinity to the northwestern populations.The second axis (MEMGENE2) supported the division of eastern Arabian populations from the rest.The plot of the relationship of the genetic and geographic distances showed a non-linear pattern of several structured populations, some of which were close both geographically and genetically (Figure6).Interestingly, the plot of the correlation between the genetic and geographic distances did not change regardless of whether we used the Euclidean distances or distances calculated as least-cost paths through the dispersal corridors.We therefore show only the latter plot.Most between-sample comparisons were separated by a geographic distance between 1000 and 2000 km and a genetic distance of about 0.05 estimated substitutions per site, which is consistent with the results of the phylogenetic analyses that indicated the presence of several clades within Bunopus sp. 4. 4 | DISCUSS ION 4.1 | Diversification within Bunopus According to our estimates of the evolutionary history of Bunopus, the crown diversification took place in the mid-Miocene, about 14 Ma (confidence interval: 10.8-17.6Ma) and resulted in the split between the Iranian and Arabian clades.The Iranian clade subsequently and gradually diversified into up to five lineages that may correspond to five distinct species (Khosravani et al., 2017).The Arabian clade radiated considerably later at 5.8 Ma (4.2-7.8Ma) and gave rise to the lineages occurring around the Arabian Gulf and in the Arabian Peninsula.This clade contains B. tuberculatus sensu stricto from southern Iran and two candidate species, a Mesopotamian one that is referred to as Bunopus sp. 3, and one that is widespread in Arabia (Bunopus sp.4).Our results that are based on a broad geographic and genetic sampling support the findings of Khosravani et al. (2017), who found that Iran supports several genetically distinct lineages within Bunopus presumably representing cryptic species.Our results show that the differentiation of the Iranian clade is older and deeper than that in the Arabian clade, and that it is quite F I G U R E 4 Contemporary habitat suitability model of Bunopus in the Arabian Peninsula (upper left panel; based on pooled records of B. tuberculatus sensu stricto, Bunopus sp. 3 and Bunopus sp.4), and habitat suitability models projected to different past time periods as indicated on top of each panel.Warmer colours denote higher probability of presence.Note that the Arabian Gulf dried out during the Quaternary sea-level low stands and the seabed provided a suitable habitat for Bunopus.F I G U R E 5 Contemporary dispersal corridors for the Bunopus geckos across the Arabian Peninsula.Large and coloured dots show sampled localities with colours corresponding to different lineages within Bunopus and matching the colours used in Figures 1, 2. Only the Arabian clade that consists of B. tuberculatus sensu stricto, Bunopus sp. 3 and sp. 4 was included in this analysis.Black dots indicate distribution records that were used for developing the potential distribution model.The population connectivity visualises landscape corridors that enable dispersal and promote gene flow between Bunopus populations.Specimen depicted in an individual from south Israel (Photo: Doubravka Velenská).

F
I G U R E 6 (a) Results of the spatial analysis of population structure conducted in MEMGENE.A total of 82 genotyped localities were used for the analysis.Circles of similar size and colour indicate individuals with similar scores along the first (left) and second (right) MEMGENE axes (large black and large white circles describe opposite extremes).The heatmap in the background shows the dispersal corridors for Bunopus across Arabia.(b) Plot of genetic distances against distances calculated as least-cost paths connecting all pairs of sampled localities among the Arabian populations of the Bunopus geckos.Warmer colours indicate higher densities of points.Note the large cluster of points at the genetic distance of about 0.05, which indicates the presence of several shallow phylogenetic lineages within the species.

4
oscillated substantially according to the prevailing climatic conditions in Arabia.For example, most of eastern Arabia seemed to have supported suitable conditions for Bunopus continuously since the mid-Pliocene, while central Arabia became inhabitable only very recently in the late Holocene (Figure4).The uninterrupted presence of Bunopus sp. 4 in eastern Arabia might have been allowed by the absence of dispersal barriers in the region.Most of the region is and has been suitable for the geckos since the mid-Pliocene, with only the massif of the Hajar Mountains always presenting an insurmountable barrier for this lowland-dwelling species.The continuous presence of Bunopus in eastern Arabia also likely explains the genetic homogeneity of local populations along the second MEMGENE axis (Figure6).Based on the suitable habitat models, the strongest environmental predictor of the genus' distribution in Arabia is the elevation.Bunopus rarely occurs above 500 m in Arabia(Šmíd   et al., 2021)  and it is thus not present in the mountain ranges that rim the peninsula: the Hajar Mountains in the east, Dhofar in the south, and the Asir and Hejaz Mountains in the west.This may be paralleled in the lineages of the Iranian clade of Bunopus that inhabit the uplifted Iranian Plateau.Although we did not include them in the distribution modelling, it is obvious from the available distribution data that they also avoid high-elevation regions such as the Zagros Mountains in the southwest of Iran(Šmíd   et al., 2014).Whether the Iranian and Arabian Bunopus lineages show some differences in the environmental niches they occupy should be addressed in a separate study.The wide belt of suitable habitat that stretches longitudinally across Arabia from the Arabian Gulf to the Red Sea constitutes the dominant contemporary dispersal corridor for Bunopus (Figure5).By connecting the eastern and western margins of the peninsula it enables longitudinal migration between geographically disparate regions with subsequent population connectivity and genetic homogenisation over this vast territory (Figure6).This corridor turns northwards at the Red Sea coast and runs to Jordan and thus provides connection between the southern Jordanian and Israeli populations with the central Arabian ones.Of note is the origin of the isolated Yemeni populations.Although they are geographically closer to those from southern Saudi Arabia, they more likely originated from southern Oman to which they are also genetically most similar.Such a biogeographic route also conforms to the general distributional patterns in the area(de Pous et al., 2016; Machado Khosravani et al. (2017) remained unresolved, and the monophyly of Bunopus was not supported.Our sampling of C. eversmanni covered broadly the distribution of the species and our phylogenetic results enable us to infer its position with more confidence.In summary, the results of all our analyses indicate that the genus Bunopus is paraphyletic with respect to C. orientalis and the genus Crossobamon is polyphyletic.Crossobamon eversmanni(Wiegmann, 1834)  is the type species of the genus Crossobamon Boettger, 1888 and as such retains its generic name.To resolve the above issue of para-and polyphyly, we propose to reassign Crossobamon orientalis to the genus Bunopus, the new combination being Bunopus orientalis comb.nov.
fragmented during moist climatic phases, and it expanded in times of increased aridity.The genus requires taxonomic revision to formally assess its diversity.Based on the results obtained in this study, Crossobamon orientalis is reassigned to Bunopus.
Tissue samples included in this study originated from targeted field trips of the authors and colleagues.They were supplemented by samples obtained from museum voucher specimens from the following collections: Zoological Research Museum Alexander Koenig, Bonn, Germany (ZFMK); Museum of Vertebrate Zoology, Berkeley, was Bunopus crassicauda was inferred to be sister to the remaining species, but the topology was only partially supported (78.3/92/1.00).The Arabian clade is formed by Bunopus tuberculatus sensu stricto from southern Iran, and the candidate species Bunopus sp. 3 from Mesopotamia and Bunopus sp. 4 from the after).It was formed by two strongly supported sister clades: Iranian (96.4/100/1.00)and Arabian (100/100/1.00).The Iranian clade consists of species occurring on the Iranian plateau and further East and North: Bunopus crassicauda, the candidate species Bunopus sp. 1, Bunopus sp. 2 and Bunopus sp. 5, and also Crossobamon orientalis from Pakistan and India.The relationships within this clade were only partially resolved.