Final datasets included genetic data for 44 specimens of D. galeatus, primarily composed of samples from the Australian Biological Tissues Collection (ABTC) at the South Australian Museum (SAMA) or new and more recent samples collected for this study and lodged at the Museum and Art Gallery of the Northern Territory (NTM) (Table S1). DNA was extracted using a Puregene DNA isolation kit (Gentra Systems, Minneapolis, MN, U.S.A.) following manufacturer protocols for DNA purification from solid tissue.

Our mtDNA dataset included 17 samples of D. galeatus, comprising 11 newly sequenced samples and six sequences already available on GenBank. For new samples, a fragment of the mitochondrial genome, including the 3’ end of the NADH dehydrogenase subunit 2 (ND2) gene and the tRNA genes was amplified and sequenced using the forward primers 5’- AAGCTTTCGGGGCCCATACC -3’ and the reverse primer 5’- CTAAAATRTTRCGGG­A­TCGAGGCC -3’. The mtDNA alignment comprised of 837 bp from the ND2 gene. These were aligned using the MUSCLE algorithm (Edgar 2004) and subsequently checked by eye for missense mutations and correct reading frames. Pairwise sequence divergences (uncorrected p distances) were calculated in MEGA 6.06 (Kumar et al. 2018). Lineage relationships within the D. galeatus complex were visualised by generating phylogenetic networks of mtDNA using the Neighbor-Net algorithm as implemented in SplitsTree v4.10 (Huson and Bryant 2006). Apomorphic diagnostic single nucleotide polymorphisms were visually identified from the ND2 alignment across members of the D. galeatus complex following methods outlined elsewhere (Renner 2016; Donnellan et al. 2023). SNP sites were numbered in reference to the first base pair of ND2 on an annotated GenBank accession of D. galeatus (AY369009). Polarity of apomorphic states was estimated relative to an exemplar of the closest lineage to the D. galeatus complex (D. tessellatus AMS R143855, JQ173631).

To infer the timing and pattern of divergences in the D. galeatus complex in a broader context we used BEAST v.2.6 (Bouckaert et al. 2014). This program assumes each tip represents an evolutionarily independent population, so we reduced sampling to exemplars of the most genetically distinctive populations of Diplodactylus, including four main lineages in the D. galeatus complex. Third codon sites were removed from BEAST analyses to reduce the potential for saturation to confound branch length estimation. Topology and timeframes of divergence were estimated using the relaxed-clock model and Yule speciation prior for 50 million generations, sampling every 10,000 generations, with the first 20% of trees discarded as burn-in. Data were partitioned by codon and the HKY+G model was applied to each partition. There are no fossil calibrations for the age of any Australian diplodactyloid geckos, so we used two normal secondary age priors taken from Brennan and Oliver (2017), specifically: a) the age of the note subtending D. galeatus and the 2n = 28 and 2 n = 30 races of D. tessellatus (mean 12 mya, sigma 2), and b) a broad normal prior on the crown age of Diplodactylus (mean 20 mya, sigma 4).

Given the overall morphological similarity between the main mtDNA lineages identified (see results), we undertook additional genetic analyses of population differentiation and relationships using SNP data generated by Diversity Array Technology (DArT). This method uses restriction-enzyme mediated genome reduction (Jaccoud 2001) prior to library construction and Next-Generation-Sequencing to sequence the most informative representations of genomic DNA sampling. The method is an alternative to whole genome sequencing, and has proven valuable for testing hypotheses of evolutionary independence and comparing genetic diversity across populations (Melville et al. 2016; Unmack et al. 2017; Gruber et al. 2018). We sampled 43 D. galeatus, as much as possible including exemplars from across the geographic range of each major mtDNA lineage. SNPs were called on data from the D. galeatus complex alone.

We used the R package ‘dartR’ (Gruber et al. 2018) to apply additional filtering to the SNP data in R version 4.1.2 (R Core Team 2022). To ensure the data quality, we removed secondary SNPs where they were restricted to a single sequenced tag using the gl.filter.secondaries() function, filtered loci by call rate (98%) using the gl.filter.callrate() function, and removed any monomorphic loci using the gl.filter.monomorphs() function. To visualise the divergence between samples and identify any potential hybrids we generated a distance-based principal coordinates analysis (PCoA) using the genetic distance matrix with the gl.pcoa() function in dartR (Gruber et al. 2018). We assessed divergence between population clusters in the PCoA based on the proportion of loci with fixed allelic differences. Fixed differences occur when two populations share no alleles at a single site. We tested for absolute fixed differences with simulation (1000 replications) using the gl.fixed.differences() function in dartR. We estimated phylogenetic relationships based on the 180 filtered SNP data in IQ-TREE version 1.6.12 (Nguyen et al. 2015) using the GTR+G model and fast bootstrapping (– bb 1000).

We conducted clustering analyses in dartR using STRUCTURE v2.3.4, which uses Bayesian Inference to determine the number of distinct genetic clusters (K) in the SNP dataset. Analyses used a range of cluster values (K = 1–5), both with and without admixture among populations. We also ran clustering analyses in ADMIXTURE v1.3 (Alexander et al. 2009) that, like STRUCTURE, uses SNP data to analyses population structure, except using a Maximum Likelihood function. ADMIXTURE v1.3 is available as a binary execution file for Linux or Mac, so we used Oracle VirtualMachine to run the binary via an Ubuntu command line. To evaluate which number of clusters is most representative of our dataset, we selected K with the highest log-likelihood in STRUCTURE, and K with the lowest Cross Validation Error in ADMIXTURE.

To understand patterns of morphological variation across populations we examined all Diplodactylus galeatus specimens held at SAMA (the type specimen was examined based on photographs) and NTM. The following characters were measured in mm (to one decimal place) using vernier callipers for all specimens except likely juveniles (<40 mm SVL): SVL – snout-vent length; TrunkL – Trunk length, between axilla and groin; TailL – tail length, from cloaca to tip (only original tails were used in analyses); TailW – tail width, at widest point; HeadL – head length, measured obliquely from tip of snout to angle of lower jaw; HeadW – head width, at widest point; HeadD – head depth, measured behind eyes and top of head; ArmL – arm length, from elbow to base of hand; LegL – knee to base of foot; IO – inter-orbital distance; OrbL – lower anterior to upper posterior of eye socket; NarEye – nare-eye distance, nostril to anterior corner of eye socket; SnEye – tip of snout to lower anterior corner of eye socket; IntNar – internarial distance between nostrils; Ros – rostral scale height; RosCre – proportion of rostral crease relative to rostral scale height; MenL – mental scale length, from mouth to posterior edge; MenW – mental scale width, at anterior edge; and Ear – width of external ear opening. We present mean and standard deviations and the range of values for each population.

We undertook body size corrections for each measurement to account for allometric growth using the allom() function in the R package ‘GroupStruct’ (Thorpe 1975; Chan and Grismer 2022). To test for differences in morphology between the three populations we ran univariate analysis of variance (ANOVA) and Tukey’s post-hoc tests for SVL and each allometry-corrected measurement. To further examine potential divergence in morphology across the populations we ran a PCA using the built-in R function prcomp() across 11 of the allometry-corrected morphometric measures (TrunkL, HeadL, HeadW, HeadD, ArmL, LegL, OrbL, SnEye, IntNar, Rostral and Ear).

To assess whether populations in the D. galeatus complex are associated with relatively less arid areas (i.e., putative refugia), we plotted climate values at species occurrence locations using annual precipitation (mm) and annual potential evapotranspiration (mm) grids from Australian Bureau of Meteorology (http://www.bom.gov.au/climate/maps/average). These climate variables form the basis of the global precipitation/evapotranspiration bioclimatic ratio for measuring aridity (Whitford and Duval 2020). Resolution of rainfall data was 0.05 degrees (approximately 5 km grain size) and evapotranspiration was 0.1 degrees (approximately 10 km grain size). We also extracted these values at 1000 points randomly selected across the two IBRA (Interim Biogeographic Regionalisation’s of Australia) regions inhabited by the D. galeatus species complex (Stony Plains and MacDonnell Ranges) as well as the intervening Finke IBRA region. All spatial data handling was completed in ArcMap 10.6.1 (ESRI) and plotting undertaken in R (R Core Team 2022).

As outlined elsewhere (Oliver et al. 2020) we consider populations to be species when they show evidence from multiple independent data sources for a history of evolutionary independence (i.e., the Generalised Lineage Concept sensu de Queiroz 2007). Under this framework we recognise as species lineages that satisfy two or more of the following criteria: (i) statistically-supported reciprocal monophyly in nDNA phylograms; (ii) evidence from SNPs for discrete population structure and lack of geneflow; or (iii) diagnostic morphological characters based on post-hoc analyses of groups that satisfy (i) or (ii). We consider that sustained lack of gene flow, supported by multiple lines of evidence and comprehensive sampling of geographic space and genes, is sufficient to delimit species.

Phylogenetic analyses of the mtDNA data highlight four main divergent lineages within the D. galeatus complex (Fig. 2A), with p distances among these ranging from 0.83–0.114 (Tables 1, S2). Multiple samples spanning the distributions of three of these four lineages were included in analyses, and in all cases lineages were monophyletic and showed much lower within-lineage divergences (max within lineage p distance was 0.045). One lineage corresponds to D. galeatus and occurs in a broad area of northern South Australia. A second lineage occurs in isolated breakaways and ranges on the border between South Australia and the Northern Territory, and corresponds to D. fyfeisp. nov. The final two lineages, occurring in the northern and southern MacDonnell Ranges, show the lowest levels of mtDNA divergence observed between major lineages (0.082–0.084), and correspond to the northern and southern ESUs of D. tjoritjarinyasp. nov., respectively. In the ND2 alignment (Table 2), 16 diagnostic (apomorphic) sites distinguish D. galeatus, seven diagnostic sites distinguish D. fyfeisp. nov. and 22 distinguish D. tjoritjarinyasp. nov.

In phylogenetic analyses of the species-level alignment of ND2 data (Fig. 2B) the monophyly of the D. galeatus complex was strongly supported, while the pattern of relationships among the four main lineages of D. galeatus complex were poorly resolved – with the exception of the sister clade pairing of the north and south lineages of D. tjoritjarinyasp. nov. The four main lineages were inferred to have diverged from each other over a relatively short timeframe centred on the Pliocene (5.3–2.4 mya).

After filtering the DArT-called SNP data (43 individuals; 254,287 loci; 35.4% missing data) there were 27,975 SNP sites from all 43 individuals in the D. galeatus species complex (1.78% missing data). Principal co-ordinates analysis with individuals as entities and loci as attributes revealed three tight clusters (Fig. 3A). Two of these corresponded to D. galeatus (n = 16) and D. fyfeisp. nov. (n = 12) and were mainly separated along axis 2 (32.6% of total variation). The other cluster and a single outlier corresponded to north (n = 14) and south (n = 1) MacDonnell Ranges populations of D. tjoritjarinyasp. nov. separated along axis 1 (37.6% of total variation). The percentage of loci with fixed differences between the well-sampled lineages (n >10) ranged from 7.9% (D. tjoritjarinyasp. nov. vs D. galeatus) to 10.4% (D. tjoritjarinyasp. nov. vs D. fyfeisp. nov.) and was significant in all cases (Table 3).

Population cluster analyses in STRUCTURE and ADMIXTURE recovered the optimal number of clusters as K = 3, splitting the data into north, central, and southern populations (Figs 3C, S1). The grouping of a single sample from Watarrka National Park (D. tjoritjarinyasp. nov. south) was consistently problematic. It grouped with the northern population in STRUCTURE without admixture, however in STRUCTURE and ADMIXTURE analyses that allowed admixture it consistently appeared admixed between all three populations (Fig. 3C).

The PCA of allometry-corrected morphometric characters revealed substantial overlap between the three species but some evidence of diverging characters in D. fyfeisp. nov. and D. tjoritjarinyasp. nov. (Fig. 4).

Mean morphometric values showed little variation across the three populations, with the exception of ear opening width (Fig. 4; Tables S3, S4). Analysis of allometry-corrected ear opening width confirmed that this character differed across populations (F₂ = 86.19, p < 0.001), with the opening on D. tjoritjarinyasp. nov. being significantly smaller than both D. fyfeisp. nov. (p = 0.001) and D. galeatus (p <0.001), and D. fyfeisp. nov. was also significantly smaller than D. galeatus (p = 0.009). Allometry-corrected rostral scale height differed across populations (F₂ = 12.49, p < 0.001), with rostral height significantly larger in D. fyfeisp. nov. than in both D. galeatus (p < 0.001) and D. tjoritjarinyasp. nov. (p = 0.001). Standardised arm length (ArmL) differed across populations (F₂ = 8.67, p < 0.001), with D. tjoritjarinyasp. nov. having significantly shorter arm length compared with D. galeatus (p = 0.046) and D. fyfeisp. nov. (p < 0.001).

There was also variation in pattern and colouration between D. tjoritjarinyasp. nov. and the two southern species. In D. tjoritjarinyasp. nov., the pale dorsal blotches always extend ≤1/4 down the torso of the animal (median = 1/8 down sides) when viewed in lateral profile, versus usually extending ≥1/4 down the torso in the southern populations (median D. galeatus = 1/3, median D. ­fyfeisp. nov. = 1/4) (Fig. 5). Animals with a continuous vertebral stripe accounted for 9.7% of specimens across the southern populations but were absent from D. tjoritjarinyasp. nov. The pale dots in the lateral regions of D. tjoritjarinyasp. nov. are small (≤3 scales in diameter) and usually not arranged in longitudinal rows (Fig. 5E,F). In contrast, D. galeatus and D. fyfeisp. nov. always have at least a few larger pale dots (>3 scales in diameter; Fig. 5A–D) and these are typically arranged in a mid-lateral row of ‘portholes’. The background dorsal and upper lateral colouration in D. galeatus and D. fyfeisp. nov. is typically a dark red (Fig. 5A–D) versus pinkish-red or reddish-brown in D. tjoritjarinyasp. nov. (Fig. 5E,F).

The main lineages in the Diplodactylus galeatus species complex inhabit distinct bioclimatic spaces, with D. galeatus occupying drier areas with less evapotranspiration, D. tjoritjarinyasp. nov. occurring in areas with both higher rainfall and evapotranspiration, and D. fyfeisp. nov. occupying an intermediate space but still with far less rainfall than D. tjoritjarinyasp. nov. (Fig. 6A). The three bioregions are positioned along a latitudinal gradient with the two occupied bioregions at either end and the uninhabited Finke bioregion in the middle (Fig. 6B). This latitudinal gradient was also reflected in the aridity measures, with the southern Stony Plains the most arid region, the northern MacDonnell Ranges the least arid, and the Finke intermediate between the two (Fig. 6C). The areas inhabited by D. fyfeisp. nov. and D. galeatus were more arid than the uninhabited regions to the north of their respective distributions.

Both the SNP and mtDNA datasets strongly support the reciprocal monophyly of the main isolated populations in the D. galeatus complex, providing evidence of historical independence. Evolutionary independence was further supported by fixed differences across thousands of loci between these populations. The observed mitochondrial genetic differences (p distances) between the three main isolate populations (9–12%) exceeds several Diplodactylus species pairs including: D. capensis and D. granariensis (8.3%), D. capensis and D. nebulosus (8.4%), D. granariensis and D. nebulosus (8.5%), and D. polyophthalmus and D. lateroides (8.5%) (Doughty and Oliver 2013). However, it is also worth noting that some widespread populations that are considered conspecific in other Diplodactylus also show divergences of this level (e.g., Diplodactylus conspicillatus; Oliver et al. 2009), emphasising that caution must be exercised when using mtDNA divergence yardsticks. While there is evidence that the three main isolate populations are diverging morphologically, there was overlap between populations for most characters. Diplodactylus galeatus and D. fyfeisp. nov. are usually diagnosable on the basis of standardised rostral scale height but are indistinguishable based on colour or pattern. Diplodactylus tjoritjarinyasp. nov. is readily diagnosed on the basis of pattern and usually diagnosable on ear width.

In this case, while morphological datasets and mtDNA only offer partial support for species recognition, with additional evidence for reciprocal monophyly and long-term absence of geneflow from multiple datasets, we argue that all three main isolate populations are recognisable as species under the generalised lineage concept.

We do not recognise the southern MacDonnell Ranges population as a separate species given poor coverage of genetic sampling, associated uncertainty about where these lineages might meet, lack of morphological divergence and slightly lower inferred mtDNA differentiation. Analyses of population structuring based on SNP data also struggled to place this population, likely because it is both somewhat divergent and poorly sampled. Given current knowledge, we split D. tjoritjarinyasp. nov. into northern and southern Evolutionary Significant Units (ESUs; Moritz 1994). The northern ESU includes all of the northern part of the MacDonnell Ranges IBRA region east and west of Alice Springs. The southern ESU is confirmed from the George Gill Range and presumably includes other sandstone ranges of the southern MacDonnell Ranges IBRA region but this needs confirmation with further genetic sampling, particularly near Hermannsburg.

Following Kaiser et al. (2013), position statements from the Australian Society of Herpetologists (ASH 2016), and in accordance with a large number of active herpetofaunal taxonomists (Wüster et al. 2021), we do not consider selected nomenclatural acts in self-published works after 1 January 2000, even if these may have priority under the rules of the International Code of Zoological Nomenclature.

To avoid as much repetition as possible in the descriptions and diagnoses of the taxa in this morphologically conserved group, we provide brief diagnoses for the species complex first, and then present more focused diagnoses and comparisons in the species accounts.

The distribution of the D. galeatus species complex spans two regions with differing climates and thus provided a unique opportunity to examine two potential mechanisms of population isolation in an arid zone lineage – climatic refugia versus ecological specialisation.

With a distribution entirely restricted to the relatively mesic uplands of the MacDonnell Ranges IBRA region, D. tjoritjarinyasp. nov. outwardly fits the climate refugia model. However, the closely related D. galeatus and D. fyfeisp. nov. are restricted to the truly arid Stony Plains IBRA region. The Finke IBRA region, which spans the gap between D. tjoritjarinyasp. nov. and the southern species as well the western gap between the two southern species, is less arid than the Stony Plains and thus does not pose an obvious climatic barrier for the species complex. Microhabitat use in the D. galeatus species complex also points to a distribution pattern decoupled from climatic refugia. Specifically, none of the D. galeatus complex are known to use deep gorges or rock crevices that likely to afford a microclimatic buffer from extreme aridity. Instead, they tend to be associated with sparsely vegetated rocky slopes. Within this habitat they have been recorded using small invertebrate (e.g., spider) burrows as diurnal shelter and their occurrence may be more limited by the presence of small rocks on which to perch when active rather than the availability of daytime shelter (Fig. 1B). Our observations also provide no evidence that activity is limited by warm conditions and we have found them active in the evening at temperatures of >34°C. These data suggest that evolutionary retreat from intensifying aridity, in and of itself, does not explain vicariance in this gecko clade.

The alternative hypothesis for the isolation of the D. galeatus species complex into three main populations is that, like many arid zone geckos (Heinicke et al. 2017), they have specialised to use a formerly more widespread substrate that has been sundered by environmental change. Diplodactlyus tjoritjarinyasp. nov. and D. fyfeisp. nov. are separated by about 250 km of the Finke IBRA region. Potentially suitable dissected tableland habitat occurs along the Finke River, but these tablelands are generally small and isolated (Wells 1969). One possibility is that suitable stony habitat linking these isolated tablelands has been covered by the sand dunes which are a dominant contemporary landscape feature of this region. However, the formation of these sand dunes has been dated to ca. 1 mya (Fujioka et al. 2009) which postdates our estimates for timing of isolation of these geckos across this area.

An alternative possible explanation of vicariance in the D. galeatus complex is that landform erosion has removed dissected tablelands that formerly linked the MacDonnell Ranges with the Beddome Range. In support of this hypothesis, Fujioka et al. (2005) found that global cooling ca. 4 mya initiated the stripping of soil mantles from silcrete tablelands in northern South Australia, which was followed by a period of active tableland erosion commencing 3 mya and slowing during the past 2 mya (Fujioka et al. 2005). Loss of connecting dissected tableland also potentially explains the isolation of D. fyfeisp. nov. from D. galeatus, with the area between about Oodnadatta and Mt Dare now mostly comprising unsuitable pavement gibber desert (Brandle 1998). This period of active erosion 2–3 mya is consistent with our estimated timing of divergence of the D. galeatus species complex during the Plio-Pleistocene. For further testing of the hypothesis of Plio-Pleistocene loss of tableland habitat, comparative genetic analyses of two threatened species with similar distributions could be invaluable. Namely, the pygopod gecko Ophidiocephalus taeniatus and the plant Acacia latzii both have relictual distributions restricted to dissected tablelands across disjunct areas (Nano et al. 2007; McDonald et al. 2012; Pedler et al. 2014).

We have been unable to find any commentary on what the drivers of the afore-mentioned widespread erosion in the Pliocene was. Paleoclimatic data from the Nullarbor deserts of southern Australia show that the early Pliocene corresponds with a mesic reversal (Sniderman et al. 2016). These data could suggest the intriguing idea that for the D. galeatus complex at least, landform erosion associated with a relatively mesic phase may have set in train isolation and divergence across the arid zone. However, current estimates for the timing of this mesic pulse (~3.5 mya) predate estimates for the onset of major landform erosion in the stony desert (2–3 mya). One potential explanation that would unify these two signals is that rates of erosion accelerated with the loss of stabilising vegetation at the end of the Pliocene wet pulse. This remains speculative and there is considerable uncertainty about processes, however, the emergent patterns shown by the D. galeatus complex certainly differ from the classic mesic refugial model of vicariance for Australia’s central uplands. This further emphasises how outwardly comparable patterns of distribution in taxa from arid ranges may be linked to a varying suite of climatic and ecological histories (Oliver and McDonald 2016).

Species in the D. galeatus complex show only slight morphological variation and SNP data were critical to testing the hypothesis of evolutionary independence. Work on geckos in the genus Gehyra from the AAZ in Western Australia has revealed a remarkable plethora of morphologically conservative taxa largely associated with stable upland areas (Doughty et al. 2018), and work on dragons in the genus Tympanocryptis has also revealed complexes of parapatric cryptic species in the AAZ (Melville et al. 2014). It seems certain that comprehensive genetic analyses will continue to reveal additional instances of cryptic speciation and overlooked short-range endemics in the AAZ, especially in substrate-specialist taxa associated with stable landforms such as mountains and other upland areas. In this context we are aware of at least two other deeply divergent and as yet undescribed species of Diplodactylus that are also associated with upland landforms around the periphery of Australia’s arid zone (Oliver pers. obs.).

The ongoing dismemberment of widespread taxa into complexes of more restricted range endemics also necessitates re-evaluation of conservation statuses. Recent attempts to locate D. tjoritjarinyasp. nov. in areas with dense Buffel grass (Cenchrus ciliaris) near Alice Springs have failed, suggesting this invasive species impacts habitat suitability for the gecko (P. McDonald pers. obs.). However, Buffel grass is a minor floristic component across most suitable habitats for D. tjoritjarinyasp. nov. (including both ESU’s; P. McDonald pers. obs.) and is thus unlikely to pose a serious conservation risk to the species, at least in the medium term. The extent of occurrence of D. fyfeisp. nov. (EOO = 654 km²) is also very small. Fortunately, our data suggest it is likely well adapted to climatic extremes, and it also occurs in an area that is sparsely inhabited and has not been subject to widespread habitat destruction or disturbance.