The demise of a wonder: Evolutionary history and conservation assessments of the Wonder Gecko Teratoscincus keyserlingii (Gekkota, Sphaerodactylidae) in Arabia

About the Authors:

Karin Tamar

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (KT); [email protected] (SC)

Affiliation: Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Marítim de la Barceloneta, Barcelona, Spain

ORCID logo https://orcid.org/0000-0003-2375-4529

Johannes Els

Roles Conceptualization, Data curation, Investigation, Project administration, Resources, Validation, Visualization, Writing – review & editing

Affiliation: Environment and Protected Areas Authority, Sharjah, United Arab Emirates

Panagiotis Kornilios

Roles Formal analysis, Investigation, Writing – review & editing

Affiliations Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Marítim de la Barceloneta, Barcelona, Spain, The Molecular Ecology Backshop, G. Lekka, Loutraki, Greece

Pritpal Soorae

Roles Conceptualization, Data curation, Investigation, Resources, Validation, Writing – review & editing

Affiliation: Environment Agency, Abu Dhabi, United Arab Emirates

Pedro Tarroso

Roles Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – review & editing

Affiliations Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Marítim de la Barceloneta, Barcelona, Spain, CIBIO/InBIO, Research Centre in Biodiversity and Genetic Resources, Universidade do Porto, Campus Agrário de Vairão, Rua Padre Armando Quintas, Vairão, Vila do Conde, Portugal

Evanthia Thanou

Roles Formal analysis, Investigation, Writing – review & editing

Affiliations The Molecular Ecology Backshop, G. Lekka, Loutraki, Greece, Section of Animal Biology, Department of Biology, School of Natural Sciences, University of Patras, Patras, Greece

John Pereira

Roles Resources, Validation, Writing – review & editing

Affiliation: Environment and Protected Areas Authority, Sharjah, United Arab Emirates

Junid Nazeer Shah

Roles Conceptualization, Data curation, Resources, Validation, Writing – review & editing

Affiliation: Natural Resource Conservation Section, Environment Department, Dubai Municipality, Dubai, United Arab Emirates

Esmat Elfaki Mohammed Elhassan

Roles Data curation, Resources, Writing – review & editing

Affiliation: Natural Resource Conservation Section, Environment Department, Dubai Municipality, Dubai, United Arab Emirates

ORCID logo https://orcid.org/0000-0001-6018-0099

Jeruel Cabadonga Aguhob

Roles Data curation, Resources, Writing – review & editing

Affiliation: Natural Resource Conservation Section, Environment Department, Dubai Municipality, Dubai, United Arab Emirates

Saoud Faisal Badaam

Roles Data curation, Resources, Writing – review & editing

Affiliation: Natural Resource Conservation Section, Environment Department, Dubai Municipality, Dubai, United Arab Emirates

Mohamed Mustafa Eltayeb

Roles Conceptualization, Data curation, Resources, Validation, Writing – review & editing

Affiliation: Natural Resource Conservation Section, Environment Department, Dubai Municipality, Dubai, United Arab Emirates

Ricardo Pusey

Roles Resources, Writing – review & editing

Affiliation: Environment Agency, Abu Dhabi, United Arab Emirates

Theodore J. Papenfuss

Roles Data curation, Resources, Writing – review & editing

Affiliation: Museum of Vertebrate Zoology, Valley Life Sciences Building, University of California, Berkeley, CA, United States of America

J. Robert Macey

Roles Data curation, Resources, Writing – review & editing

Affiliation: Peralta Genomics Institute, Chancellor’s Office, Peralta Community College District, Oakland, CA, United States of America

Salvador Carranza

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (KT); [email protected] (SC)

Affiliation: Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Marítim de la Barceloneta, Barcelona, Spain

ORCID logo https://orcid.org/0000-0002-5378-3008

Introduction

Arid and semi-arid regions cover almost 20% of the world’s land surface and are important for understanding global biodiversity patterns [1–4]. Biodiversity in arid regions is often perceived as poor and as less valuable than in other biogeographic regions (but see [5]), thus it is often neglected and attracts less scientific attention and conservation concern [3]. Ecosystems of arid areas, despite their wide extension and limited accessibility [6], are severely affected by anthropogenic disturbances. Threats to arid environments continue to erode local biodiversity, mostly due to urban expansion, overgrazing, desertification, agricultural expansion, drought, and armed conflicts [7–9].

The United Arab Emirates (UAE) in the Arabian Peninsula, similar to other oil producing countries of the Middle East, is developing rapidly. As a result of this extensive development, human populations are increasing, threatening wildlife populations and natural habitats [9,10]. Rapid construction and industrialization in Arabia have expanded widely, especially in recent decades, destroying and degrading unique arid habitats, particularly coastal sand dunes and salt flats (sabkhahs) [11]. Destruction or degradation of natural habitats and loss of original cover have continuously driven declines in biodiversity and ecosystem services worldwide [12,13]. Fragmentation and destruction of natural habitats through anthropogenic activity are known threats to the viability of endangered species and are recognised as important drivers to current species extinctions [14,15]. Loss of natural habitats and their fragmentation result in a landscape composed of isolated patches of suitable habitats, separated by an inhospitable matrix, and thus to the isolation of previously connected populations. This segregation leads to population size reduction, which in turn, increases the probability of extinction by demographic and/or environmental changes [16,17].

One of the most widely used indicators for assessing the health of ecosystems and their biodiversity is the conservation status of flora and fauna, which is also an important component of priority-setting exercises for species conservation. Reptiles (i.e., non-avian sauropsids) comprise a major element of the global biodiversity [18] and constitute the world’s most diverse group of terrestrial vertebrates, especially in arid regions [19]. Their population trends thus represent an indicator of a general decay of environmental health, as well as the declines of other species [20]. However, despite their relevance for many ecosystems and role as excellent models for conservation studies, knowledge or assessments for the conservation of reptiles are generally lacking, falling behind those of birds, mammals, and amphibians [19,21,22].

The drastic population decline of the Wonder Gecko, Teratoscincus keyserlingii Strauch, 1863, in the UAE provides particular evidence for strong anthropogenic pressures [9] and thus should be used as a model for conservation in arid regions. Teratoscincus keyserlingii is a member of the family Sphaerodactylidae and is one of nine recognised species of the genus Teratoscincus, distributed across desert regions from South-West Asia to Central Asia (Fig 1) [18]. Teratoscincus keyserlingii is mainly distributed in Iran, South Afghanistan and West Pakistan, with a small population in the Arabian Peninsula [18,23]. It is a psammophilous species and is often found between sand dunes, but occasionally also recorded on hard soil [24–26]. In Arabia, this species is the sole representative of the genus, being isolated from con-specific populations and con-geners by the Arabian Gulf and the Strait of Hormuz. It is locally restricted to the UAE, where it occurs over a distribution range (extent of occurrence) of less than 5,000 km2, with an area of occupancy of less than 100 km2 [27]. It occurs mostly within 30–40 km off the Gulf Coast, inhabiting sand sheets, low sand dunes, and sandy plains with vegetation below 200 m of elevation (Fig 2) [26,28]. The UAE T. keyserlingii population has been recorded since the early 1970’s [29], and currently is seldom seen being restricted to just a few localities. This population exhibits severe decline and is highly localised, divided into sub-populations, which at present are severely fragmented from each other.

[Figure omitted. See PDF.]

Fig 1. Distribution ranges of Teratoscincus taxa.

Colours correspond to species in Figs 3 and 4. Localities of specimens sampled in this study are detailed in S1 Table. Teratoscincus sistanensis was not included in this study. Credits: Natural Earth contributors, SRTM.

https://doi.org/10.1371/journal.pone.0244150.g001

[Figure omitted. See PDF.]

Fig 2. Localities of T. keyserlingii in the UAE used in this study.

Images presenting the UAE coastline in the years 1984 and 2020 indicating the massive development in recent decades, which destroyed and fragmented the natural habitat of T. keyserlingii. Localities of specimens used in this study are indicated in red, with numbers corresponding to localities detailed in S1 Table. Localities 6, 7, 9, 12 are extinct due to urban development. Left map credits: OpenStreetMap contributors, SRTM. Middle and right images credits: Satellite images of environmental change, Earth Resources Observatory and Science (EROS) Center, USGS.

https://doi.org/10.1371/journal.pone.0244150.g002

Conservation status of reptiles of the Arabian Peninsula assessed by the IUCN stated that Teratoscincus keyserlingii was Regionally Endangered due to its restricted range in the UAE and habitat loss, owing to massive estate development [9] (see Fig 2). Recently, in the regional IUCN Red List specific assessment of T. keyserlingii in the UAE (the still unpublished conclusions of the 2018 UAE Red List of Reptiles and Amphibians, IUCN unpublished data), the species was evaluated as Critically Endangered. Its Arabian population is declining rapidly and predicted to decline by over 80–90% in the future. Major threats to the species in the UAE include habitat loss due to extensive and intensive development and urbanization, leading to severe fragmentation and isolated small populations. Minor threats include off-road driving to target animals due to local stigma, collection for international trade, and feral cats. Thus, the Arabian population experiences declines in the extent of occurrence, area of occupancy, and quality of habitat. Ongoing and planned estate development across the distribution of T. keyserlingii in the UAE may well drive the species to the brink of extinction in the Arabian Peninsula.

For effective conservation planning, genetic characterization is playing a considerable role in understanding the likely impact of threats for ensuring the long-term viability of a species [12,30]. Thus, in this study we describe the Arabian population dynamics within the framework of the genus–we generate a geographically broad phylogeny of Teratoscincus, with the broadest sampling to date of T. keyserlingii in the UAE. We evaluate phylogenetic structure and biogeographic history of the genus using multi-locus data and a time-calibrated Bayesian coalescent species tree. As conservation planning requires detailed genomic structure, we use genomic data and population genetics analyses to infer the T. keyserlingii Arabian population assemblage. We aim to help promote proper conservation assessments and implications that will be applied for the long-term survival of the Arabian population of T. keyserlingii.

Materials and methods

Ethics statement

No in vivo experiments were performed. Research in the UAE was conducted and specimens were sampled and manipulated, with the authorization and under strict control, and permission, of the Environmental Agencies of the three Emirates that have participated in the study: Environment and Protected Areas Authority, Government of Sharjah; Natural Resources Conservation Section, Environment Department, Dubai Municipality, Government of Dubai; Environment Agency, Government of Abu Dhabi. No animals were killed as a result of the study and all efforts were made to minimise animal suffering when obtaining a tissue sample for genetic analyses. Tissue samples from outside the UAE used for phylogenetic analyses were legally loaned from museum collections at CAS San Francisco, USA; MVZ Berkeley, USA; and KU Lawrence, USA.

This research is not institutional. As a result of its characteristics and the total control and compliance with laws, regulations, and procedures of this kind of conservation study in the UAE, this study did not need an approval by an Institutional Animal Care and Use Committee (IACUC) or ethics committee.

Spatial data

Maps were generated using QGIS v.3.14.15 Pi [31]. The global map was extracted from the public domain Natural Earth (http://www.naturalearthdata.com) with elevation data at 1 km resolution. Spatial layers with the UAE political borders, its seven emirates and coastline were extracted from the open data initiative OpenStreetMap (https://www.openstreetmap.org). Digital elevation model (DEM), used for the elevation data, was downloaded from the reverb tool from NASA (http://reverb.echo.nasa.gov). Elevation data originated from the Shuttle Radar Topography Mission at a spatial resolution of 1 arc-second (~30 m). UAE Environmental change images were extracted from USGS Earth Resources Observatory and Science (EROS) Center [32]. Terrestrial protected areas of the UAE shapefile (S1 Appendix) was provided by the authors from the Environmental Agency from Abu Dhabi, the Environment and Protected Areas Authority, Sharjah, and the Natural Resource Conservation Section, Environment Department, Dubai Municipality.

Taxon sampling, DNA sequencing, and alignment preparation

To evaluate phylogenetic structure within our focal species, T. keyserlingii, we gathered a dataset of the entire genus. We included a total of 123 specimens of Teratoscincus in our analyses, 92 newly sequenced specimens and 31 sequences retrieved from GenBank: 12 of T. bedriagai, 13 of T. microlepis, six of T. roborowskii, 12 of T. przewalskii, one of T. rustamowi, 10 of T. scincus, three of T. mesriensis, and 66 of T. keyserlingii. Monophyly of Teratoscincus has been supported in several studies (e.g., [33–35]), thus we used another member of the family Sphaerodactylidae, Pristurus rupestris, as an outgroup. Detailed information on specimens, GenBank accession numbers, and localities is presented in S1 Table.

Genomic DNA of alcohol-preserved tissue samples was extracted following standard high-salt protocols [36]. We amplified via Polymerase Chain Reaction (PCR) and bi-directionally sequenced two mitochondrial gene fragments, Cytochrome c oxidase subunit I (COI) and NADH Dehydrogenase subunit 2 (ND2), and two nuclear gene fragments, Melanocortin 1 receptor (MC1R) and Recombination activating gene 1 (RAG1). Primers used, their sequences, lengths of amplified regions, PCR conditions and sources are given in S2 Table. Chromatograms were checked and assembled, and the final alignments were generated using MAFFT v.7 [37] in Geneious v.9.0.5 [38]. No stop codons were detected after translation and heterozygous positions in the nuclear genes were coded according to the IUPAC ambiguity codes.

Phylogenetic analyses

Phylogenetic relationships of Teratoscincus were inferred using the concatenated mtDNA-nucDNA and the COI datasets, performing Maximum Likelihood (ML) and Bayesian Inference (BI) analyses. The concatenated dataset was used for the general phylogenetic structure of Teratoscincus, and the COI dataset allowed broader geographic sampling and allocation of the UAE specimens to a particular T. keyserlingii lineage. Partitions and substitution models were selected using PartitionFinder v.2 [39] with the Bayesian Information Criterion (BIC) and the dataset partitioned by genes: COI (TrN+G), ND2, RAG1 (HKY+G), and MC1R (HKY+I). ML analyses were performed in raxmlGUI v.1.5 [40], with the GTRGAMMA model and 100 random addition replicates. Nodal support was assessed with 1,000 bootstrap replicates. BI analyses were performed in MrBayes v.3.2.6 [41] with all parameters unlinked across partitions. Analyses were performed with four chains per run for 2 million generations, with sampling of every 200 generations and burn-in of the first 25% of posterior trees of each run. We calculated inter- and intraspecific uncorrected p-distance between and within Teratoscincus taxa for COI and ND2 mitochondrial fragments, with pairwise deletion, in MEGA v.7 [42].

Estimating a temporal framework of divergence

We estimated divergence times for Teratoscincus with a species-tree approach using *BEAST v.1.8.4 [43]. We assigned “species” to lineages from the phylogenetic trees described above. We calibrated the species-tree using biogeographical data from Macey et al. (1999, 2005) [44,45], supported in other studies using different calibration approaches [e.g., 34,35,46,47]: the split between T. microlepis and the remaining Teratoscincus species due to the rise of the Hindu Kush ~20 Mya (normal distribution; mean 21, stdev 3) and the T. scincus-T. roborowskii divergence caused by the Tien Shan–Pamir uplift ~10 Mya (normal distribution; mean 11, stdev 3; including also T. przewalskii, T. keyserlingii, and T. mesriensis). We analysed the data with the following priors (otherwise by default): COI+ND2 and RAG1 (HKY+G), MC1R (K80+G), alpha prior uniform (0–10), uncorrelated relaxed clock for the mitochondrial partition and strict clock for the nuclear loci (uniform distribution; mean 0.01 and 0.0001, respectively, 0–1). We run *BEAST three times for 200 million generations and a sampling frequency of every 20,000 generations, with 10% discarded as burn-in. Convergence, posterior trace plots, and effective sample sizes (>200) were evaluated with Tracer v.1.6 [48]. LogCombiner and TreeAnnotator were used to combine the runs and to generate the ultrametric tree.

Ancestral range reconstruction

To infer biogeographic history and estimate ancestral ranges of Teratoscincus, we used a reduced dataset. Specimens selected represent lineages delimited after applying the multi-rate Poisson Tree Processes (mPTP) [49] model, and one additional specimen representing a lineage occurring in more than one discrete biogeographic area (i.e., T. keyserlingii in South-West Asia and Arabia). As the mPTP analysis relies on single locus data, we reconstructed and used a ML mitochondrial phylogenetic tree as specified above (i.e., same mitochondrial partitions and analysis parameters). Ancestral area reconstruction analysis included a reduced concatenated mtDNA-nucDNA dataset of representatives and performed using the Bayesian Stochastic Search Variable Selection (BSSVS) [50] model implemented in BEAST [51]. Partitions, models of nucleotide substitution, and priors for divergence times were as detailed above. Prior settings, MCMC chain length and sampling strategy were as in the divergence-time estimations analysis, with additional specification of symmetric discrete trait substitution model and an exponential prior for the discrete location state rate. Teratoscincus representatives were assigned to four discrete biogeographic areas based on their modern day distributions: (i) Central Asia (China, Mongolia); (ii) West Asia (Turkmenistan, Uzbekistan); (iii) South-West Asia (Iran, Afghanistan, Pakistan); (iv) Arabia (UAE).

ddRAD library preparation and sequencing

We collected ddRAD data for 26 samples of Teratoscincus keyserlingii from the UAE (see S1 Table), following the protocol, adapters, and indices of Peterson et al. (2012) [52]. We double-digested 500 ng of genomic DNA for each sample with 20 units each of two restriction enzymes (SbfI and MspI, New England Biolabs) for 6 h at 37°C. Fragments were purified with Agencourt AMPure beads before ligation of barcoded Illumina adapters. Samples with unique adapters were pooled, and each pool of eight samples was size-selected for a distribution of fragments with a peak of 480 bp, using E-Gel SizeSelect 2% agarose gels (Thermo Fisher Scientific). Illumina multiplexing indices were ligated to individual samples using a Phusion polymerase kit (high-fidelity Taq polymerase, New England Biolabs). Final pools were sequenced on an Illumina NexSeq 500, under a 50 bp single-end read protocol at the UPF Genomics Core Facility, Barcelona, Spain.

ddRAD bioinformatics and genome-wide SNPs

We processed raw Illumina reads using the program iPyRAD v.0.7.8 [53]. We demultiplexed samples using their unique barcode and adapter sequences. Sites with Phred quality scores under 99% (Phred score = 33) were changed into ‘N’ characters, and reads with ≥3 N’s were discarded. Within the iPyRAD pipeline, filtered reads for each sample were clustered using VSEARCH v.2.4.3 [54] and aligned with MUSCLE v.3.8.31 [55]. We assembled ddRAD data using a clustering threshold of 92%. As an additional filtering step, consensus sequences that had low coverage (<6 reads), excessive undetermined or heterozygous sites (>3) or too many haplotypes (>2) were discarded. Consensus sequences were clustered across samples using the within-sample clustering threshold (92%). Again, alignment was done with MUSCLE, applying a paralog filter that removes loci with excessive shared heterozygosity among samples (paralog filter = 200). Maximum number of SNPs per locus was set to 6. We generated final datasets that included no missing data, i.e., with all loci present for all samples. Concatenated ddRAD data generated in this study are available in S2 Appendix.

Two data-matrices were used in the analyses: the first included the complete concatenated sequence of all invariant and variable sites (“ddRAD” dataset), while the second included one random SNP from each putatively unlinked locus (“uSNPs” dataset).

Population structure analyses and measures of genetic variation

We estimated genetic variation parameters, including percentage of polymorphic loci and percentage of polymorphic sites in the total alignment. The concatenated ddRAD dataset was used in MEGA v.7 [41] to estimate nucleotide diversity π [56] and to perform the neutrality test of Tajimaʼs D [57] for the UAE T. keyserlingii population. The R package hierfstat [58] was used for the estimation of basic F-statistics, specifically per locus expected and observed heterozygosity (He and Ho, respectively), overall fixation index (FST) and inbreeding index (FIS), as well as the hybridization coefficient (FIC) for each studied individual. We tested whether loci departed significantly from Hardy-Weinberg (HW) equilibrium with a Fisher’s exact test and 1,000 permutations. Finally, for clusters of individuals, we estimated pairwise FST values [59] and pairwise Nei’s D genetic distances [56], which corrects for sample size.

To investigate population structure, we used the uSNPs dataset and ran a discriminant analysis of principal components (DAPC) implemented in the R package adegenet [60,61]. We inferred the number of clusters (max K = 10) with the find.cluster function using the Bayesian Information Criterion (BIC) and used the optim.a.score function to optimise the number of PCAs.

We used the same uSNPs dataset in a Bayesian clustering analysis, as implemented in STRUCTURE v.2.3.4 [62], in order to detect population structure under specified models of evolution and to define admixed individuals. Given close geographic proximity of sampling locations, we avoided the assignment of individuals into pre-defined populations, and expected some degree of gene-flow throughout the studied region. Thus, we applied a model that assumed correlated allele frequencies and admixture among clusters and tested the probability of different numbers of independent genetic clusters (K = 1–10) with 15 iterations (250.000 burn-in length and 500.000 post burn-in replicates). Distribution of allele frequencies lambda (λ) was estimated for our dataset, setting K = 1. As for prior alpha (a), specifying the degree of admixture between each sub-population, the default setting presumes clusters that are expected to contribute equally to the sample and share the same estimated value of alpha. However, unbalanced samples (different sample size per K), extensive admixture or low FST values (<0.05) can seriously affect the number of inferred Ks and the use of an alternative setting (different alpha per K, initial a<1) can better detect the true number of clusters in most cases [63]. Thus, we opted for this alternative setting with initial a = 0.3. Inferred number of population clusters was estimated by examining log-likelihood per K [lnP(X/K)] [62] and Evanno’s ΔK [64], in Clumpak online web server [65], with the ‘greedy’ option and 2,000 random input orders.

Results

Phylogenetic and biogeographic inferences

Our phylogenetic analyses included 123 individuals of Teratoscincus. The dataset totalling 2,868 bp comprised mitochondrial fragments of COI (n = 93, 651 bp), ND2 (n = 98, 546 bp), and nuclear fragments of MC1R (n = 52, 663 bp), and RAG1 (n = 56, 1,008 bp). Uncorrected p-distances of COI and ND2 mitochondrial fragments between and within each species/lineage are summarised in S3 Table.

The ML and BI phylogenetic trees of Teratoscincus using the concatenated mtDNA-nucDNA and COI datasets (Fig 3), the time-calibrated species tree (Fig 4), and the ancestral range reconstructions (Fig 5; mPTP tree is in S1 Fig) show a similar phylogenetic structure. Teratoscincus began diverging approximately 19.3 Mya (95% highest posterior density [HPD]: 13.9–25.1 Mya), most likely in South-West Asia (63% probability) with subsequent splits to north and north-west (i.e., West and Central Asia). A sister-taxa relationship between T. bedriagai and T. microlepis is recovered in the species tree and biogeographic analysis, but not supported in the concatenated mtDNA-nucDNA tree or at all recovered in the COI tree. Teratoscincus microlepis shows high intraspecific diversity and is divided into two distinct lineages, which split around 5.6 Mya (95% HPD: 2.6–9.6 Mya) and have a genetic distance of 8.3% in COI and 8.5% in ND2 (S3 Table). A sister taxon relationship between T. roborowskii and T. przewalskii was supported in all trees, diverging ca. 3.5 Mya (95% HPD: 1.3–6.1 Mya). The remaining species comprise a well-supported clade diverging approximately 4 Mya (95% HPD: 2.4–6.2 Mya), yet relationships among them are not supported in all trees.

[Figure omitted. See PDF.]

Fig 3. Maximum likelihood phylogenetic trees of Teratoscincus.

(left) Phylogenetic tree using the concatenated mtDNA–nucDNA data (COI, ND2, MC1R, RAG1). (right) Phylogenetic tree using the COI data. Support values indicated near the nodes (ML bootstrap/Bayesian posterior probabilities). Sample codes correlate to specimens in S1 Table. Specimen depicted is T. keyserlingii from the UAE. Photo credit: Salvador Carranza.

https://doi.org/10.1371/journal.pone.0244150.g003

[Figure omitted. See PDF.]

Fig 4. Time-calibrated species tree of Teratoscincus.

Mean age estimates are provided above the nodes with horizontal bars representing the 95% highest posterior densities. Black circles represent nodes with posterior probability values ≥0.95. Specimen depicted is T. keyserlingii from the UAE. Photo credit: Salvador Carranza.

https://doi.org/10.1371/journal.pone.0244150.g004

[Figure omitted. See PDF.]

Fig 5. Ancestral area reconstructions of Teratoscincus.

A Bayesian phylogenetic tree estimated using the concatenated mtDNA–nucDNA data with BSSVS. A pie chart describing the probability of each inferred area is presented near the major nodes. Branch colours indicate inferred ancestral range and the mean age estimates are provided above the nodes with horizontal bars representing the 95% highest posterior densities. Black circles represent nodes with posterior probability values ≥0.95.

https://doi.org/10.1371/journal.pone.0244150.g005

Teratoscincus keyserlingii is paraphyletic with respect to T. mesriensis, and shows great intraspecific diversity, with four (I–IV) and five (A–E) distinct lineages in the concatenated mtDNA-nucDNA and COI datasets, respectively. These lineages comprise specimens from Pakistan (lineage I/A), Iran (lineages II/B, C, III/D), and Iran and the UAE (lineage IV/E). Genetic distance among these lineages of T. keyserlingii ranged between 3.4–9.3% in COI and 5.1–9.4% in ND2. Lineage IV/E of T. keyserlingii comprised all the specimens sampled from across the UAE. These specimens clustered together with one specimen from southern Iran collected near the city of Bandar Abbas on the shores of the Strait of Hormuz (Fig 1 and S1 Table). Divergence within this lineage, between southern Iran and the UAE was recovered in our analyses at 0.5 Mya.

Genome-wide data

For the 26 specimens of T. keyserlingii from the UAE, we obtained more than 70 million reads, with the number of reads per individual ranging from 1.17 to 4.03 million (mean 2.71). The total number of pre-filtered loci was 34,487, which dropped to 33,123 after applying all quality control filters, and to 3,912 in the final dataset, after applying a conservative filter of “all loci present in all samples”. Most loci were found to be monomorphic (3,542), while 282 included one SNP, 56 had two SNPs, 28 had three SNPs and only four loci were found to include the maximum observed number of four polymorphisms. Consequently, the final concatenated ddRAD dataset was 136,177 bp long and included a total of 494 SNPs and 370 unlinked SNPs (uSNPs).

Population structure and genetic variation within T. keyserlingii from the UAE

Percentages of polymorphic loci and polymorphic sites in our concatenated ddRAD dataset were approximately 0.095 and 0.004, respectively, and the estimated π value was 1.28x10-4. Tajima’s D test of neutrality resulted in a value of -2.233. Basic F-statistics estimated from the uSNPs dataset showed a significant departure from the HW equilibrium in less than 5% of the loci (p<0.05), but an excess of homozygosity in most of them. Average FST and FIS values over all loci were 0.19 and -0.18, respectively. The individual inbreeding coefficient (FIC) ranged between 0.44–0.52.

The optimal number of clusters obtained from DAPC was K = 1 (S2 Fig). STRUCTURE analysis recovered a slightly higher possibility for K = 3 (mean LnProb = -6,275 and mean similarity score = 0.926 for all 15 iterations) than K = 1 (mean LnProb = -6,403) (S3 Fig). This result persisted under both alternative settings of the alpha prior and even when different evolutionary models were tested (uncorrelated frequencies or no-admixture model). Most individuals were grouped in one cluster widely distributed in all sampling locations. Each of the other two presumed clusters grouped two individuals; in each case these were specimens sampled in the same exact location or in locations less than 300 m apart. Individuals belonging to the widespread cluster were also captured in the same locations, while some admixture was found between the widespread and the two local clusters (S3 Fig). Genetic differentiation among clusters was high, with Nei’s D ranging from 0.15 between the widespread cluster and each of the two local ones to 0.47 between the two local clusters. Pairwise FST values were similarly high (0.11–0.53).

Discussion

The coupled phylogenetic-genomic approach carried out in this study for Teratoscincus accomplishes the task of combining broad taxonomic and geographic sampling for time-calibrated species tree inference to investigate the origin and evolution of the genus, together with including valuable information gained from genomic data for the population structure of T. keyserlingii in the UAE.

Evolutionary history of Teratoscincus

Our emphasis on geographic range-wide diversification patterns enabled us to estimate relationships among almost all described Teratoscincus species and identify natural groupings within taxa. Congruent phylogenetic structure and monophyletic groups of species were recovered in previous studies of Teratoscincus [44,45,66], supported by studies incorporating other taxa (e.g., [33–35,46,67]), and corroborated by our study. We recovered monophyly of most recognised species, along with high genetic and geographic diversity within T. microlepis and our non-monophyletic focal species, T. keyserlingii. Sister-species topology is mainly supported in the species tree and the biogeographic reconstruction, though with only weak or no support in the concatenated mtDNA-nucDNA and COI trees. This is especially evident in geographically proximate species, and at least in one sympatric location, of T. microlepis and T. bedriagai. Divergence time estimations using the concatenated mtDNA-nucDNA dataset were generally similar though older, as expected, from the species tree.

Geographic limits of the major phylogenetic groupings of Teratoscincus correspond to mountain ranges uplifting in Asia caused by the Indo-Eurasian and Arabian-Eurasian tectonic collisions (Fig 1). These continuous ridges likely created geographic barriers restricting past and current species’ distributions, promoting vicariance, separating the genus into biogeographical regions. Macey et al. (1999, 2005) [44,45], having representatives of major phylogenetic groupings within Teratoscincus, previously hypothesised tectonic scenarios responsible for key divergences within the genus. Topology and timing of the proposed orogenic events have been supported in later studies (e.g., [34,35,46,67]), and the broader geographic and taxonomic sampling in this study further strengthen their assessments. Divergence within Teratoscincus began in the early Miocene and likely originated in South-West Asia. An ancestral population probably diverged due to tectonic uplifts of the Hindu Kush and Karakorum mountain ranges in northern Afghanistan and Kashmir [68–70]. This event divided Teratoscincus into two–a southern group of T. microlepis and T. bedriagai in the general region of Afghanistan-Iran, and a northern one, the ancestral population of the remaining species. Rise of the Tien Shan and Pamir mountain ranges in western China and eastern Tajikistan during the middle Miocene [71–73] likely promoted an east-west divergence of remaining species–an eastern population dispersed into western China and southern Mongolia, and a western population spread to the south and east of the Caspian Basin. The eastern population diverged into two recognised species, T. roborowskii and T. przewalskii, after the rise of the Kuruk Tagh mountain range (an eastern extension of the Tien Shan range) in the middle Pliocene [74–76], restricting T. roborowskii to the Turpan Basin in China. The western population diverged into northern T. scincus and southern T. keyserlingii, likely as a result of the Kopet-Dagh mountain uplifting along the Iran-Turkmenistan border in the early Pliocene [77,78].

Taxonomic concerns are raised regarding T. microlepis and T. keyserlingii due to considerable intraspecific diversity within them. Teratoscincus microlepis presents a late Miocene diversification, with one lineage (lineage II/B) restricted to the Hamun-e Jaz Murian Basin in southeast Iran, where it is sympatric with T. keyserlingii, spanning across the southern provinces of Kerman and Sistan and Baluchistan. A second lineage (lineage I/A) is widespread in eastern and southern Iran, southern Afghanistan and western Pakistan. Mitochondrial genetic distance values between these lineages are greater than between the two recognised species T. roborowskii and T. przewalskii (S3 Table). These lineages thus should be taxonomically revised including specimens from the type locality of T. microlepis in Iran. A taxonomic revision should also include T. sistanensis, which is morphologically closest to T. microlepis [79]. These two species are distinguished by the number of dorsal scales around mid-body (85–110 in T. microlepis and 145–165 in T. sistanensis) [79]. Each lineage of T. microlepis from our study exhibits characters belonging to the two species (examined specimens from the MVZ collection by TJP; lineage I/A, 98–188 dorsal scales around mid-body and lineage II/B, 87–165). Since T. sistanensis was described based on morphological data alone and the lack of additional distinguishing characters known, further studies are needed to investigate these taxa, considering the type localities of the two species.

A second case is the non-monophyly of T. keyserlingii, with respect to T. mesriensis, and the presence of five distinct lineages within it. The lineages recovered within T. keyserlingii, including that of T. mesriensis, diverged in the Pliocene and Pleistocene, and are geographically structured within Iran and neighbouring countries. For both T. microlepis and T. keyserlingii we hypothesise that the tectonic consequences of the northward motion of the Arabian plate and collision with the Eurasian landmass, resulted in intense geological instability and uplift of mountain ranges, as well as associated climatic changes prevailing in the Miocene to Pleistocene, most likely prompted these diversifications [80]. A positive association between intraspecific diversity, exhibited as genetic structure and/or structured lineage distribution, and habitat discontinuity is expected [81]. We therefore call for thorough systematic revisions of these two species, with broad geographic sampling of specimens to account for the large amount of spatial and genetic diversity, with the plausible finding of additional lineages.

Origin of the Arabian population of T. keyserlingii

The UAE specimens of T. keyserlingii cluster together in lineage IV/E with one specimen from southern Iran, near the coastal city of Bandar Abbas. Low levels of genetic divergence between the coastal southern Iran sample and the Arabian population leads us to hypothesise that Arabia was colonised from southern Iran relatively recently, during the Pleistocene. Connections between Arabia and Iran have been recorded for several reptile taxa, including Acanthodactylus [82], Mesalina [83,84], Hemidactylus [85,86], and Pristurus [87]. As has been hypothesised in the case of Pseudocerastes persicus [88], colonization could have occurred by transmarine dispersal across the Strait of Hormuz, or directly by land dispersal during the last glaciations when sea levels were low [89]; another plausible scenario is human-mediated dispersal, although being a secretive, sand-dwelling, nocturnal reptile makes this hypothesis less likely. Additionally, recent divergence is indicated by extremely low genetic variability within the Arabian population of T. keyserlingii. Several UAE specimens differentiate by a single mutation in the mitochondrial data indicating rapid recent expansion across the UAE.

Population structure of T. keyserlingii in the UAE

Genetic structure analysis with DAPC shows the UAE T. keyserlingii population lacking any internal geographic structure, hence representing a single population. In fact, this pattern is in agreement with the genetic variation measures (percentage of polymorphic loci and polymorphic sites) that highlight the extremely low genetic diversity within the studied system. Recently, a study that also used a ddRAD approach estimated within-population π values for more than 80 lizard populations ranging from 2.14x10-4 to 3.99x10-3 [90]. Our estimated π value within T. keyserlingii in the UAE is much lower than the lowest estimated value in the aforementioned study. Accordingly, Tajima D test of neutrality returned a negative value of -2.233. Negative values are indicative of recent selective sweeps or, more possibly in our case, population expansions after a recent bottleneck or a founder event, which is congruent with our divergence time estimation of Arabian colonization during the late Pleistocene.

Results from DAPC and STRUCTURE seem contradictory at first glance: the former identified no genetic structure within T. keyserlingii in the UAE, while the latter detected three sub-populations. The DAPC analysis does not require a model of admixture and describes the relationship between the population clusters, optimizing the variance between groups while minimizing variance within groups. In this sense, it is expected to perform better when organisms exhibit continuous spatial population structure [91]. On the other hand, the algorithm implemented in STRUCTURE assumes a specific evolutionary model, under which a more or less equal number of individuals are randomly sampled from K populations, each of which is characterised by a set of allele frequencies at each locus. It is further assumed that within populations, loci are at HW and linkage equilibrium, while the effect of genetic drift and lack of gene-flow will lead to linkage disequilibrium and departures from Hardy-Weinberg proportions, when there is population structure. The STRUCTURE algorithm uses this signal to detect the “true” number of populations [62].

Any major violation of the former assumptions can lead to severe under- or overestimation of the number of K. For example, representing some populations with only very few individuals is one possible violation and simulation studies have shown that this may affect the estimated K [63]. Following Wang (2017) [63], our analyses were run under a scheme aiming to correct for unbalanced sampling violations, without altering the outcome of K = 3. However, there is strong evidence that this clustering does not represent “true” population structure of T. keyserlingii in the UAE. Firstly, where several values of K give similar estimates of LnP(X/K), as here for K = 1 and K = 3, the smallest value is often “correct” [92]. Secondly, the split to three clusters most probably reflects departures from the STRUCTURE model that are detected as HW or linkage disequilibrium within the assumed populations. Inbreeding, specifically, can lead to a weak statistical signal for K>1, even in the absence of population structure [93]. The mean inbreeding coefficient (FIC) per individual in our sample showed high values of FIC>0.45 (Fig 6). This indicates an almost 50% probability that individuals with common alleles at any locus show this genetic similarity by descent from genetically related parents. This would be a serious violation of the random mating assumed by the STRUCTURE model, forcing the algorithm to group together the most closely related individuals and in a distinct cluster all non-related (or less related) ones. Upon inspection of the aligned consensus sequences, we found that each of the two local clusters defined by STRUCTURE included individuals that shared an extremely high number of common alleles in approximately 94% of the total polymorphic sites, or 348–350 of 370 total uSNPs. The effect of inbreeding also explains the unreasonably high divergence estimated among the three clusters; classical F statistics consider independent populations and are not suitable for populations connected with extensive gene-flow [94].

[Figure omitted. See PDF.]

Fig 6. Inbreeding frequency distribution.

Calculated for the 26 individuals of T. keyserlingii from the UAE.

https://doi.org/10.1371/journal.pone.0244150.g006

Conclusively the bulk of evidence suggests that there is no true structure within our samples and that T. keyserlingii specimens from the UAE belong to a single population. This population shows a low level of genetic diversity, an excess of homozygotes, and high levels of inbreeding, all consistent with a recent founder event. Non-random mating leading to high values of FIC can be induced by bottlenecks or founder events, small population size and habitat fragmentation in natural populations (e.g., [95,96]), but it can also be the result of breeding programs on captive animals [97]. In any case, increased inbreeding may lead to the fixation of deleterious alleles and is associated with higher mortality, lower fertility and even extinction [98,99].

Conservation implications for T. keyserlingii in the UAE

Effects of anthropogenic activities on natural environments and biodiversity by habitat destruction and fragmentation are considered major drivers for the decline and categorization of Teratoscincus keyserlingii as Critically Endangered in the UAE. Its decline is related to the environmental characteristics of small patches of natural ecosystems: they comprise lower environmental heterogeneity, they are occupied by small populations that are more vulnerable to local extinction, they have limited immigration rate, and their edges represent degraded habitat for many species [15].

Biodiversity conservation starts with knowledge of genetic diversity within and among populations. A high level of genetic diversity in a population is thought to increase the probability of surviving a disastrous event or extinction [12]. Thus, in order to implement effective management policies, i.e., to prioritise protection status or to determine re-introduction or translocation options, it is imperative to know the genetic structure of focal populations. Our genomic conclusions of a single population of T. keyserlingii in the UAE, with no internal structure, rather simplify possible immediate conservation strategies (i.e., not accounting for the loss of genetic diversity). This genetic inference therefore has significant implications for conservation management.

The UAE currently includes 42 protected areas, covering 16.6% of the terrestrial area of the country [28], albeit most of the protected range is located at the southern edge of the Abu Dhabi Emirate (Fig 7). Burriel-Carranza et al. (2019) [28] calculated the percentage of T. keyserlingii distribution included within these protected areas, showing that it did not reach the strict Aichi biodiversity conservation target of 17% and only slightly surpassed the less restrictive 12% target. Their study indicated an immediate problem regarding habitat loss of T. keyserlingii and its low presence in protected areas. Currently, three protected areas (Jabal Ali, Al Marmoun Desert, and Misanad) incorporate areas of T. keyserlingii distribution in the UAE, three include past localities of the species (Ed-dhelaimah, Meleiha and Elfaya, and Ras Ghanada), and the rest are either too small and scattered, or located beyond the natural habitat of the species (Fig 7). Unfortunately, most protected areas covering currently observed specimens of T. keyserlingii are small, containing minimal overlap with the species’ overall distribution in the UAE, are located near developing areas, or are too scattered to maintain natural environments with habitat connectivity. Species distribution modelling presented in Burriel-Carranza et al. (2019) [28] indicated that the distribution of T. keyserlingii may be extended to southern Dubai and northern Abu Dhabi Emirates. As protection of the environment is being recognised as an essential component of the natural development policies in the UAE, one achievable proposition is expansion of current protected areas (e.g., northern expansion of Al Marmoun Desert reserve in Dubai), or establishment of new areas with suitable habitats for T. keyserlingii (e.g., in northern Abu Dhabi Emirate, west of Al Marmoun Desert reserve). These spatial zones are of critical importance to the natural preservation of the species in the UAE. In this regard, improved land management and natural habitat restoration in disturbed areas are also plausible efforts to improve the quality, connectivity, and permeability for future dispersal of different species, as well as preserving the biological diversity and stability of the ecosystem [100].

[Figure omitted. See PDF.]

Fig 7. Distribution of T. keyserlingii in the UAE within the current terrestrial protected areas.

Terrestrial protected areas are indicated in green (see details in Burriel-Carranza et al., 2019) [28]. Circles indicate T. keyserlingii specimens used in this study (red) and past known localities (blank). Credits: OpenStreetMap contributors, SRTM. Shapefile of the terrestrial protected areas of the UAE is provided in S1 Appendix.

https://doi.org/10.1371/journal.pone.0244150.g007

In the UAE there are currently two breeding centres for endangered Arabian wildlife, including T. keyserlingii–in Abu Dhabi and Sharjah Emirates. The goal of many ex-situ conservation agendas, such as captive breeding or re-introduction programs, is retention of a high percentage of genetic diversity from wild populations, until its potential re-introduction to natural habitat to ensure the long-term viability of the founder populations [101–103] (see also [104] for further references and discussion). However, increased inbreeding originating from a small number of stock parents may lead to fixation of deleterious alleles and is associated with higher mortality, lower fertility and even extinction [98,99]. Additionally, the strategy of releasing inbred animals locally may result in extreme local homogeneity. Both cases will result in fixation of alleles throughout the UAE population, or at least locally. Coupled with small home ranges and limited dispersal distances per generation, re-introduction should thus be taken into consideration for future release of inbred animals and a more efficient re-location strategy. Thus, ex-situ conservation programs have to be carefully integrated within in-situ programs, so that captive populations can constitute an insurance against the loss of natural populations and a source for re-introduction projects [15]. Yet, when the decline of species has reached catastrophic dimensions, as is predicted for T. keyserlingii, captive breeding and other ex-situ conservation programs should not be considered alone as efficient conservation measures without additional in-situ measures [105], and must not become an excuse to avoid dealing with the preservation of natural habitats.

Several decades of habitat destruction and fragmentation have led to the gradual disappearance of the Arabian population of T. keyserlingii. As shown in Fig 2, between 1984 and 2020, massive development in the coastal areas of the Emirates of Abu Dhabi, Dubai, Sharjah, Umm Al Quwain and Ras al Khaimah have exterminated many localities where T. keyserlingii had been previously reported (Fig 7) and have artificially fragmented the species’ distribution range. The genomic inferences in this study enable the re-location of specimens from current areas of occupancy, from destroyed habitats or breeding centres, to areas where the probability of future persistence is predicted to be higher. Setting up of a solid and continuous education and awareness campaign will help to ensure constant, but increasingly necessary, actions for the preservation of T. keyserlingii in Arabia. In addition, the protection of T. keyserlingii will indisputably simultaneously promote protection of other co-occurring species. The primary habitat of T. keyserlingii consists of sand sheets, low undulating sand dunes and sandy plains with relatively dense vegetation dominated by the grasses Pennisetum divisum and Panicum turgidum [26]. These vegetated coastal sandy areas of the UAE represent a unique ecosystem in urgent need of protection, that is the natural habitat of a large community of geckos, as well as of many other animals both vertebrate and invertebrate [26,28].

Supporting information

[Figure omitted. See PDF.]

S1 Fig. The mPTP results of Teratoscincus.

Results inferred from the concatenated mitochondrial dataset. Sample codes correlate to specimens in S1 Table and in Figs 2 and 5.

https://doi.org/10.1371/journal.pone.0244150.s001

(DOCX)

S2 Fig. Bayesian values vs. Number of clusters as resulted from DAPC analysis.

Analysis based on 26 individuals of T. keyserlingii from the UAE used for the ddRAD analyses. The Bayesian values do not improve as the clusters increase from 1 to 10 implying K = 1.

https://doi.org/10.1371/journal.pone.0244150.s002

(DOCX)

S3 Fig. Bayesian genetic clustering of the 26 individuals of T. keyserlingii from the UAE using the program STRUCTURE.

(A) The optimal clustering was identified using Evanno’s method (ΔK plotted against K) at K = 3. (B) The division or runs is also given (similar clustering in any subset of the total 15 runs per K).

https://doi.org/10.1371/journal.pone.0244150.s003

(DOCX)

S1 Table. Information on the Teratoscincus specimens included in the analyses and related GenBank accession numbers.

Localities of specimens from the UAE are indicated next to the country’s name and in Fig 2. Specimen code abbreviations: [CAS] California Academy of Sciences, San Francisco, USA; [CN] Institute of Evolutionary Biology, Barcelona, Spain; [JFBM] J. F. Bell Museum of Natural History, University of Minnesota, USA; [KU] University of Kansas Biodiversity Institution, Lawrence, USA; [MVZ Herp] Museum of Vertebrate Zoology, University of California, Berkeley, USA; [TG] Tony Gamble collection; [ZMMU] Zoological Museum of Moscow State University, Moscow, Russia. (*) Teratoscincus keyserlingii specimens used for the ddRADseq analyses (n = 26).

https://doi.org/10.1371/journal.pone.0244150.s004

(DOCX)

S2 Table. Data on the gene fragments used in this study.

Information including the primers used, with their orientation, sequences, references, and PCR conditions.

https://doi.org/10.1371/journal.pone.0244150.s005

(DOCX)

S3 Table. Pairwise uncorrected mitochondrial sequence divergence (p-distance) among and within Teratoscincus.

Values derived from the mitochondrial gene fragments of COI (below diagonal) and ND2 (above diagonal), and within each taxon or lineage (in bold; COI/ND2).

https://doi.org/10.1371/journal.pone.0244150.s006

(DOCX)

S1 Appendix. Shapefile of the terrestrial protected areas of the UAE.

https://doi.org/10.1371/journal.pone.0244150.s007

(ZIP)

S2 Appendix. ddRAD data generated in this study.

https://doi.org/10.1371/journal.pone.0244150.s008

(RAR)

Acknowledgments

We wish to thank Eng. Dawoud AbdulRahman AlHajri (Director General: Dubai Municipality), Khalid Mohammed Sharif Al-Awadhi (Executive Director: Dubai Health, Safety and Environment Sector), Eng. Alya Abdulrahim Abdulla Amin Alharmoudi (Director: Dubai Environment Department), Aisha Almurr Al Muhery (Section Manager: Dubai Natural Resources Conservation) for their help. We also wish to thank for the support of her Excellency Hana Saif Al Suwaidi, Chairperson of the Environment and Protected Areas Authority, Sharjah and Paul Vercammen, Operations Manager of the Breeding Centre for Endangered Arabian Wildlife. We thank Lauren Scheinberg (California Academy of Sciences, San Francisco, California, USA), Carol L. Spencer (Museum of Vertebrate Zoology, Berkeley, California, USA), and Luke Welton (University of Kansas, Lawrence, Kansas, USA) for their help with the shipment of tissue samples, and to Bernat Burriel-Carranza for his help generating the maps.

Citation: Tamar K, Els J, Kornilios P, Soorae P, Tarroso P, Thanou E, et al. (2021) The demise of a wonder: Evolutionary history and conservation assessments of the Wonder Gecko Teratoscincus keyserlingii (Gekkota, Sphaerodactylidae) in Arabia. PLoS ONE 16(1): e0244150. https://doi.org/10.1371/journal.pone.0244150

References

1. Trabucco A, Zomer RJ. Global aridity index (global-aridity) and global potential evapotranspiration (global-PET) geospatial database. CGIAR Consortium for Spatial Information; 2009. Available at: http://www.csi.cgiar.org.

2. Brito JC, Godinho R, Martínez‐Freiría F, Pleguezuelos JM, Rebelo H, Santos S, et al. Unravelling biodiversity, evolution and threats to conservation in the Sahara‐Sahel. Biol Rev. 2014; 89, 215–231. pmid:23848599

3. Durant SM, Pettorelli N, Bashir S, Woodroffe R, Wacher T, De Ornellas P, et al. Forgotten biodiversity in desert ecosystems. Science. 2012; 336, 1379–1380. pmid:22700901

4. Durant SM, Wacher T, Bashir S, Woodroffe R, De Ornellas P, Ransom C, et al. Fiddling in biodiversity hotspots while deserts burn? Collapse of the Sahara’s megafauna. Divers Distrib. 2014; 20, 114–122. https://doi.org/10.1111/ddi.12157.

5. Brito JC, Tarroso P, Vale CG, Martínez‐Freiría F, Boratyński Z, Campos JC, et al. Conservation Biogeography of the Sahara-Sahel: additional protected areas are needed to secure unique biodiversity. Divers Distrib. 2016; 22(4), 371–384. https://doi.org/10.1111/ddi.12416.

6. Ficetola GF, Bonardi A, Sindaco R, Padoa‐Schioppa E. Estimating patterns of reptile biodiversity in remote regions. J Biogeogr. 2013; 40(6), 1202–1211. https://doi.org/10.1111/jbi.12060.

7. Ghabbour SI. Threats to biodiversity in Arab countries. In: Reviews in Ecology: Desert Conservation and Development. Baraket HN, Hegazy AK editors. Cairo: Metropole Press; 1997. pp. 129–157.

8. Ward D. Biology of Deserts. Oxford: Oxford University Press; 2009.

9. Cox NA, Mallon D, Bowles P, Els J, Tognelli MF. The Conservation Status and Distribution of Reptiles of the Arabian Peninsula. Cambridge, UK and Gland, Switzerland: IUCN, and Sharjah, UAE: Environment and Protected Areas Authority: 2012.

10. Vine P. Natural Emirates: wildlife and environment of the United Arab Emirates. London: Trident Press Ltd; 1996.

11. Edgell HS. Arabian Deserts: Nature, Origin and Evolution. Dordrecht: Springer; 2006.

12. Frankham R, Ballou JD, Briscoe DA. Introduction to conservation genetics. Cambridge: Cambridge university press; 2002.

13. Butchart SH, Walpole M, Collen B, Van Strien A, Scharlemann JP, Almond REA, et al. Global biodiversity: indicators of recent declines. Science. 2010; 328(5982), 1164–1168. pmid:20430971

14. Groombridge B. Global biodiversity: status of the earth's living resources. London: Chapman and Hall; 1992.

15. Hunter M. Jr Benchmarks for managing ecosystems: are human activities natural?. Conserv Biol. 1996; 10(3), 695–697.

16. Burkey TV. Extinction rates in archipelagos: implications for populations in fragmented habitats. Conserv Bio. 1995; 9, 527–541. https://doi.org/10.1046/j.1523-1739.1995.09030527.x.

17. Noss RS, Csuti B, Groom MJ. Habitat fragmentation. In: Groom MJ, Meffe GK, Carroll CR, editors. Principles of conservation biology; 2012. pp. 213–250.

18. Uetz P, Freed P, Hošek J. The Reptile Database. http://www.reptile-database.org. (Accessed February 2020).

19. Roll U, Feldman A, Novosolov M, Allison A, Bauer A, Bernard R, et al. The global distribution of tetrapods reveals a need for targeted reptile conservation. Nat Ecol Evol. 2017; 1, 1677–1682. pmid:28993667

20. Todd BD, Willson JD, Gibbons JW. The global status of reptiles and causes of their decline. In: Sparling DW, Linder G, Bishop CA, Krest S. editors. Ecotoxicology of Amphibians and Reptiles. Boca Raton, USA: CRC Press; 2010. pp. 47–67.

21. Böhme M, Collen B, Baillie JEM, Bowles P, Chanson J, Cox N, et al. The conservation status of the world’s reptiles. Bio Conserv. 2013; 157, 372–385. https://doi.org/10.1016/j.biocon.2012.07.015.

22. Meiri S, Chapple DG. Biases in the current knowledge of threat status in lizards, and bridging the ‘assessment gap’. Biol Conserv. 2016; 204, 6–7. https://doi.org/10.1016/j.biocon.2016.03.009.

23. Szczerbak NN, Golubev ML. Gecko fauna of the USSR and contiguous regions. New York: Society for the Study of Amphibians & Reptiles; 1996.

24. Minton SA. A contribution to the herpetology of West Pakistan. Bull Am Mus Nat Hist. 1966; 134, 28–184.

25. Anderson SC. The lizards of Iran. Society for the Study of Amphibians and Reptiles. Ohio: Oxford; 1999.

26. Gardner AS. The Amphibians and Reptiles of Oman and the UAE. Frankfurt am Main: Chimaira; 2013.

27. Soorae PS, Mohamed AA, Els J, O’Donovan D, Pedo V Jr, Alzahlawi N, et al. Wonder Gecko Conservation Acton Plan. December 2017 Workshop Report. Environment Agency‐ Abu Dhabi and Conservation Planning Specialist Group; 2018.

28. Burriel-Carranza B, Tarroso P, Els J, Gardner A, Soorae P, Mohammed AA, et al. An integrative assessment of the diversity, phylogeny, distribution, and conservation of the terrestrial reptiles (Sauropsida, Squamata) of the United Arab Emirates. PLoS One. 2019; 14(5), e0216273. pmid:31048886

29. Arnold EN. The scientific results of the Oman flora and fauna survey 1975. Little-known geckos (Reptilia: Gekkonidae) from Arabia with descriptions of two new species from the Sultanate of Oman. J Oman Studies, Special Report. 1977; 1, 81–110.

30. Boulding EG. Genetic diversity, adaptive potential, and population viability in changing environments. In: Carroll SC, editor. Fox conservation biology: evolution in action. New York: Oxford University Press; 2008. pp 199–221.

31. QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation; 2020. http://qgis.osgeo.org.

32. USGS, Earth Resources Observatory and Science (EROS) Center. Earthshots: Satellite images of environmental change–Dubai, United Arab Emirates. 2015. Available from: http://earthshots.usgs.gov/node/76.

33. Gamble T, Bauer AM, Greenbaum E, Jackman TR. Evidence for Gondwanan vicariance in an ancient clade of gecko lizards. J Biogeogr. 2008; 35(1), 88–104. https://doi.org/10.1111/j.1365-2699.2007.01770.x.

34. Skipwith PL, Bauer AM, Jackman TR, Sadlier RA. Old but not ancient: coalescent species tree of New Caledonian geckos reveals recent post-inundation diversification. J Biogeogr. 2016; 43(6), 1266–1276. https://doi.org/10.1111/jbi.12719.

35. Zheng Y, Wiens JJ. Combining phylogenomic and supermatrix approaches, and a time-calibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species. Mol Phylogenet Evol. 2016; 94, 537–547. pmid:26475614

36. Sambrook J, Fritsch EF, Maniatis T. Molecular cloning: a laboratory manual. 2nd ed. Cold spring harbor laboratory press; 1989.

37. Katoh K, Misawa K, Kuma KI, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002; 30(14), 3059–3066. pmid:12136088

38. Kearse M, Moir R, Wilson A, Stones‐Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics, 2012; 28, 1647–1649. pmid:22543367

39. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2017; 34(3), 772–773. pmid:28013191

40. Silvestro D, Michalak I. raxmlGUI: A graphical front-end for RAxML. Org Divers Evol. 2012; 12, 335–337. https://doi.org/10.1007/s13127-011-0056-0.

41. Ronquist F, Teslenko M, Mark P, van der Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012; 61, 1–4. pmid:21828081

42. Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol Biol Evol. 2016; 33(7), 1870–1874. pmid:27004904

43. Heled J, Drummond AJ. Bayesian inference of species trees from multilocus data. Mol Biol Evol. 2010; 27, 570–580. pmid:19906793

44. Macey JR, Wang Y, Ananjeva NB, Larson A, Papenfuss TJ. Vicariant patterns of fragmentation among gekkonid lizards of the genus Teratoscincus produced by the Indian collision: a molecular phylogenetic perspective and an area cladogram for Central Asia. Mol Phylogenet Evol. 1999; 12(3), 320–332. pmid:10413626

45. Macey JR, Fong JJ, Kuehl JV, Shafiei S, Ananjeva NB, Papenfuss TJ, et al. The complete mitochondrial genome of a gecko and the phylogenetic position of the Middle Eastern Teratoscincus keyserlingii. Mol Phylogenet Evol. 2005; 36, 188–193. pmid:15904865

46. Nielsen SV, Bauer AM, Jackman TR, Hitchmough RA, Daugherty CH. New Zealand geckos (Diplodactylidae): cryptic diversity in a post-Gondwanan lineage with trans-Tasman affinities. Mol Phylogenet Evol. 2011; 59(1), 1–22. pmid:21184833

47. Metallinou M, Arnold EN, Crochet PA, Geniez P, Brito JC, Lymberakis P, et al. Conquering the Sahara and Arabian deserts: systematics and biogeography of Stenodactylus geckos (Reptilia: Gekkonidae). BMC Evol Biol. 2012; 12(1), 258. https://doi.org/10.1186/1471-2148-12-258.

48. Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. 2014. Available at: http://beast.bio.ed.ac.uk/Tracer.

49. Kapli P, Lutteropp S, Zhang J, Kobert K, Pavlidis P, Stamatakis A, et al. Multi-rate Poisson Tree Processes for single-locus species delimitation under Maximum Likelihood and Markov Chain Monte Carlo. Bioinformatics. 2016; 33, 1630–1638. https://doi.org/10.1093/bioinformatics/btx025.

50. Lemey P, Rambaut A, Drummond A, Suchard M. Bayesian phylogeography finds its roots. PLoS Comput Biol. 2009; 5, e1000520. pmid:19779555

51. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012; 29, 1969–1973. pmid:22367748

52. Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One. 2012; 7, e37135. pmid:22675423

53. Eaton DAR. PyRAD: assembly of de novo RADseq loci for phylogenetic analyses. Bioinformatics, 2014; 30, 1844–1849. pmid:24603985

54. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. A versatile open source tool for metagenomics. PeerJ, 2016; 4, e2584. pmid:27781170

55. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004; 32, 1792–1797. pmid:15034147

56. Nei M. Molecular Evolutionary Genetics. New York: Columbia University Press; 1987.

57. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989; 123(3), 585–95. pmid:2513255

58. de Meeûs T, Goudet J. A step-by-step tutorial to use HierFstat to analyse populations hierarchically structured at multiple levels. Infect Genet Evol. 2007; 7, 731–735. pmid:17765664

59. Weir BS, Cockerham CC. Estimating F-Statistics for the Analysis of Population Structure. Evolution, 1984; 38 (6), 1358–1370. pmid:28563791

60. Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008; 24: 1403–1405. pmid:18397895

61. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010; 11, 94. pmid:20950446

62. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000; 155, 945–959. pmid:10835412

63. Wang J. The computer program structure for assigning individuals to populations: easy to use but easier to misuse. Mol Ecol Resour. 2017; 17, 981–990. pmid:28028941

64. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005; 14(8), 2611–2620. pmid:15969739

65. Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015; 15(5), 1179–1191. pmid:25684545

66. Nazarov RA, Radjabizadeh M, Poyarkov NA Jr, Ananjeva NB, Melnikov DA, Pouyani ER. A new species of frog-eyed gecko, genus Teratoscincus strauch, 1863 (squamata: sauria: sphaerodactylidae), from central Iran. Russ J Herpetol. 2017; 24(4), 291–310.

67. Heinicke MP, Jackman TR, Bauer AM. The measure of success: geographic isolation promotes diversification in Pachydactylus geckos. BMC Evol Biol. 2017; 17(1), 9. pmid:28077086

68. Searle MP. Geology and tectonics of the Karakoram Mountains. Chichester, UK: John Wiley and Sons; 1991.

69. Le Fort P. Granites in the tectonic evolution of the Himalaya, Karakorum and southern Tibet. In: Shackleton RM, Dewey JF, Windley BF, editors. Tectonic Evolution of the Himalayas and Tibet. Philos T R Soc A. 1998; 326, 281–299.

70. Hildebrand PR, Noble SR, Searle MP, Waters DJ, Parrish RR. Old origin for an active mountain range: Geology and geochronology of the eastern Hindu Kush, Pakistan. Geol Soc Am Bull. 2001; 113(5), 625–639. https://doi.org/10.1130/0016-7606(2001)113<0625:OOFAAM>2.0.CO;2.

71. Molnar P, Tapponnier P. Cenozoic tectonics of Asia: effects of a continental collision. Science, 1975; 189(4201), 419–426. pmid:17781869

72. Tapponnier P, Mattauer M, Proust F, Cassaigneau C. Mesozoic ophiolites, sutures, and arge-scale tectonic movements in Afghanistan. Earth Planet Sci Lett. 1981; 52(2), 355–371.

73. Abdrakhmatov KY, Aldazhanov SA, Hager BH, Hamburger MW, Herring TA, Kalabaev KB, et al. Relatively recent construction of the Tien Shan inferred from GPS measurements of present-day crustal deformation rates. Nature. 1996; 384(6608), 450–453. https://doi.org/10.1038/384450a0.

74. Tapponnier P, Molnar P. Active faulting and Cenozoic tectonics of the Tien Shan, Mongolia, and Baykal regions. J Geophys Res: Solid Earth. 1979; 84(B7), 3425–3459.

75. Windley BF, Allen MB, Zhang C, Zhao ZY, Wang GR. Paleozoic accretion and Cenozoic redeformation of the Chinese Tien Shan range, central Asia. Geology, 1990; 18(2), 128–131.

76. Trifonov VG, Ivanova TP, Bachmanov DM. Recent mountain building of the central Alpine-Himalayan Belt. Geotectonics, 2012; 46(5), 315–332. https://doi.org/10.1134/S0016852112050068.

77. Sborshchikov IM, Savostin LA, Zonenshain LP. Present plate tectonics between Turkey and Tibet. Tectonophysics, 1981; 79(1–2), 45–73.

78. Shabanian E, Siame L, Bellier O, Benedetti L, Abbassi MR. Quaternary slip rates along the northeastern boundary of the Arabia-Eurasia collision zone (Kopeh Dagh Mountains, Northeast Iran). Geophys J Int. 2009; 178, 1055–1077. https://doi.org/10.1111/j.1365-246X.4183.x.2009.

79. Akbarpour M, Shafiei S, Sehhatisabet ME, Damadi E. A new species of frog-eyed gecko, genus Teratoscincus Strauch, 1863 (Squamata: Sphaerodactylidae), from southeastern Iran. Zool Middle East. 2017; 63, 296–302. https://doi.org/10.1080/09397140.2017.1388490.

80. Macey JR, Schulte JA II, Ananjeva NB, Van Dyke ET, Wang Y, Orlov N, et al. A molecular phylogenetic hypothesis for the Asian agamid lizard genus Phrynocephalus reveals discrete biogeographic clades implicated by plate tectonics. Zootaxa. 2018; 4467, 1–81. pmid:30313432

81. Avise JC. Phylogeography: the history and formation of species. Harvard university press; 2000.

82. Tamar K, Carranza S, Sindaco R, Moravec J, Trape JF, Meiri S. Out of Africa: Phylogeny and biogeography of the widespread genus Acanthodactylus (Reptilia: Lacertidae). Mol Phylogenet Evol. 2016; 103, 6–18. pmid:27395778

83. Šmíd J, Moravec J, Gvoždík V, Štundl J, Frynta D, Lymberakis P, et al. Cutting the Gordian Knot: Phylogenetic and ecological diversification of the Mesalina brevirostris species complex (Squamata, Lacertidae). Zool Scr. 2017; 46(6), 649–664. https://doi.org/10.1111/zsc.12254.

84. Simó-Riudalbas M, Tamar K, Šmíd J, Mitsi P, Sindaco R, Chirio L, et al. Biogeography of Mesalina (Reptilia: Lacertidae), with special emphasis on the Mesalina adramitana group from Arabia and the Socotra Archipelago. Mol Phylogenet Evol. 2019; 137, 300–312. pmid:31029747

85. Carranza S, Arnold EN. A review of the geckos of the genus Hemidactylus (Squamata: Gekkonidae) from Oman based on morphology, mitochondrial and nuclear data, with descriptions of eight new species. Zootaxa. 2012; 3378, 1–95. http://dx.doi.org/10.11646/zootaxa.3378.1.1.

86. Šmíd J, Carranza S, Kratochvíl L, Gvoždík V, Nasher AK, Moravec J. Out of Arabia: a complex biogeographic history of multiple vicariance and dispersal events in the gecko genus Hemidactylus (Reptilia: Gekkonidae). PLoS One. 2013; 8, e64018. pmid:23724016

87. Badiane A, Garcia-Porta J, Cervenka J, Kratochvíl L, Sindaco R, Robinson MD, et al. Phylogenetic relationships of Semaphore geckos (Squamata: Sphaerodactylidae: Pristurus) with an assessment of the taxonomy of Pristurus rupestris. Zootaxa. 2014; 3835, 33–58. pmid:25081434

88. de Pous P, Simó-Riudalbas M, Els J, Jayasinghe S, Amat F, Carranza S. Phylogeny and biogeography of Arabian populations of the Persian horned viper Pseudocerastes persicus (Duméril, Bibron & Duméril, 1854). Zool Middle East. 2016; 62(3), 231–238. https://doi.org/10.1080/09397140.2016.1202896.

89. Lambeck K. Shoreline reconstructions for the Persian Gulf since the last glacial maximum. Earth Planet Sci Lett. 1996; 142, 43–57.

90. Singhal S, Huang H, Title PO, Donnellan SC, Holmes I, Rabosky DL. Genetic diversity is largely unpredictable but scales with museum occurrences in a species-rich clade of Australian lizards. P Roy Soc B: Biol Sci. 2017; 284, 20162588. https://doi.org/10.1098/rspb.2016.2588.

91. Royal CD, Novembre J, Fullerton SM, Goldstein DB, Long JC, Bamshad MJ, et al. Inferring genetic ancestry: opportunities, challenges, and implications. Am J Hum Genet. 2010; 86, 661–673. pmid:20466090

92. Pritchard JK, Wen X, Falush D. Documentation for structure software: Version 2.3. 2009. Available from http://www.ccg.unam.mx/~vinuesa/tlem09/docs/structure_doc.pdf.

93. Falush D, Stephens M, Pritchard JK. Inference of population structure: Extensions to linked loci and correlated allele frequencies. Genetics. 2003; 164, 1567–1587. pmid:12930761

94. Weir BS, Goudet J. A Unified Characterization of Population Structure and Relatedness. Genetics, 2017; 206(4), 2085–2103. pmid:28550018

95. Andersen LW, Fog K, Damgaard C. Habitat fragmentation causes bottlenecks and inbreeding in the European tree frog (Hyla arborea). P Roy Soc B-Biol Sci. 2004; 271, 1293–1302. https://doi.org/10.1098/rspb.2004.2720.

96. Osborne AJ, Negro SS, Chilvers BL, Robertson BC, Kennedy MA, Gemmell NJ. Genetic Evidence of a Population Bottleneck and Inbreeding in the Endangered New Zealand Sea Lion, Phocarctos hookeri. J Hered. 20161; 107, 392–402. pmid:26995741

97. Putnam AS, Ivy JA. Kinship-Based Management Strategies for Captive Breeding Programs When Pedigrees Are Unknown or Uncertain. J Hered. 2014; 105, 303–311. pmid:24143031

98. Theodorou K, Couvet D. On the expected relationship between inbreeding, fitness, and extinction. Genet Sel Evol. 2006; 38(4), 371–87. pmid:16790228

99. Thévenon S, Couvet D. The impact of inbreeding depression on population survival depending on demographic parameters. Anim Conserv. 2002; 5(1), 53–60. https://doi.org/10.1017/S1367943002001075.

100. Mawdsley JR, O’Malley R, Ojima DS. A review of climate‐change adaptation strategies for wildlife management and biodiversity conservation. Conserv Biol. 2009; 23(5), 1080–1089. pmid:19549219

101. Utter F, Epifanio J. Marine aquaculture: genetic potentialities and pitfalls. Rev Fish Biol Fisher. 2002; 12, 59–77. https://doi.org/10.1023/A:1022644021870.

102. Frankham R. Genetic adaptation to captivity in species conservation programs. Mol Ecol. 2008; 17, 325–333. pmid:18173504

103. Edwards T, Cox EC, Buzzard V, Wiese C, Hillard LS, Murphy RW. Genetic assessments and parentage analysis of captive bolson tortoises (Gopherus flavomarginatus) inform their ‘‘rewilding” in New Mexico. PLoS One. 2014; 9, e102787. pmid:25029369

104. Valbuena-Ureña E, Soler-Membrives A, Steinfartz S, Alonso M, Carbonell F, Larios-Martín R, et al. Getting off to a good start? Genetic evaluation of the ex situ conservation project of the Critically Endangered Montseny brook newt (Calotriton arnoldi). PeerJ, 2017; 5, e3447. pmid:28626615

105. Snyder NF, Derrickson SR, Beissinger SR, Wiley JW, Smith TB, Toone WD, et al. Limitations of captive breeding in endangered species recovery. Conserv Biol. 1996; 10(2), 338–348.