Effective biodiversity conservation planning starts with the genetic characterization within and among focal populations, in order to understand the likely impact of threats for ensuring the long-term viability of the species. The Wonder Gecko, Teratoscincus keyserlingii, is one of nine members of the genus. The species is distributed in Iran, Afghanistan and Pakistan, with a small isolated population in the United Arab Emirates (UAE), where it is classified nationally as Critically Endangered. Within its Arabian range, anthropogenic activity is directly linked to the species’ decline, highly localized and severely fragmented populations. Here we describe the evolutionary history of Teratoscincus, by reconstructing its phylogenetic relationships and estimating its divergence times and ancestral biogeography. For conservation implications of T. keyserlingii we evaluate the genetic structure of the Arabian population using genomic data. The presentstudy supports the monophyly of most species and reveals considerable intraspecific variability in T. microlepis and T. keyserlingii, which necessitate broad systematic revisions. The UAE population of T. keyserlingii likely arrived from southern Iran during the Pleistocene and no internal structure was recovered within, implying a single population status. Regional conservation of T. keyserlingii requires improved land management and natural habitat restoration in the species’ present distribution, and the expansion of current protected areas, or the establishment of new areas with suitable habitats for the species, mostly in northern Abu Dhabi Emirate.

We collected ddRAD data for 26 samples of Teratoscincus keyserlingii from the UAE, following the protocol, adaptors, and indices of Peterson et al. (2012). 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.

We processed raw Illumina reads using the program iPyRAD v.0.7.8. 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, the filtered reads for each sample were clustered using VSEARCH v.2.4.3 and aligned with MUSCLE v.3.8.31. We assembled the ddRADseq 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. The 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). The 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.