• Open Access

Higher genetic diversity is associated with stable water refugia for a gecko with a wide distribution in arid Australia


Correspondence: Paul E. Duckett, Department of Biological Sciences, Macquarie University, Balaclava Road, Sydney, NSW 2109, Australia.

E-mail: paul.duckett@mq.edu.au



Those locations and features providing refugia for species during unfavourable climatic phases may be important to identify and conserve to help protect biodiversity in the near future. During the Pleistocene, climates oscillated between glacial and interglacial periods. In the Northern Hemisphere, the impact from glacial ice sheets caused species to disperse to warmer southern refugia, but less is known about changes to species' distributions during these periods in the Southern Hemisphere.


In Australia, the climate alternated between highly arid and wetter periods during the Pleistocene. It has been proposed that arid refugia may be associated with the inland ranges (areas of higher relief), the mesic east, or areas that maintained favourable species-specific ecological conditions. We tested these hypotheses with a phylogeographical analysis of a widely distributed tree-dwelling gecko, Gehyra variegata (2n = 40a chromosomal race) throughout the central and eastern regions of arid Australia.


We generated a mtDNA sequence and microsatellite dataset by sampling 740 G. variegata40a throughout its known distribution. We also use species distribution modelling to predict the species' likely past, present and future distribution.


The majority of G. variegata40a lineages diverged during the Pleistocene, and those located in regions of arid Australia where mean annual water run-off is highest, displayed higher levels of genetic diversity in comparison locations with lower run-off. We also show that genetic diversity increased with proximity to water sources.

Main conclusions

It is likely that G. variegata40a contracted to refugia associated with stable water sources during Pleistocene arid phases. However, modelling suggests that unfavourable climate conditions will be present in this region by 2070. Therefore, Pleistocene refugia for G. variegata40a are unlikely to be refugia in the future. More generally, our results suggest that water run-off could be a useful predictor to identify favourable conditions for some arid species.


Global climate change has profoundly affected the natural distribution and abundance of biological diversity (Hughes, 2000; Walther et al., 2002; Thuiller et al., 2005). With global temperatures now increasing at unprecedented rates, the natural distribution of species is expected to be severely impacted [Thomas et al., 2004; Intergovernmental Panel on Climate Change (IPCC), 2007]. In recognition of this threat, conservation programmes are beginning to value and utilize the growing body of climate change literature to facilitate the preservation of biodiversity (Rouget et al., 2006; Vos et al., 2008). One important approach to help reduce the future impact on biodiversity is knowledge of historical processes that have influenced the evolutionary and colonization histories of species (Byrne et al., 2008). Molecular data can provide useful indications of historical processes because demographic change will be reflected in the distribution and abundance of genetic diversity. These processes are likely to have influenced connectivity, created population bottlenecks and/or opened up new ecological niches across a species' distribution. Thus, throughout a species' distribution, there is likely to be varying amounts of genetic diversity resulting from differences in demographic history. Identifying and conserving the locations of high intraspecific genetic diversity is one strategy needed to ensure that the evolutionary potential of a species is preserved for the future (Scoble & Lowe, 2010).

During the Pleistocene epoch, c. 2,588,000–11,700 yr bp (Gibbard et al., 2010), climatic conditions oscillated between dry and wet phases. How these conditions impacted on species throughout this period has been speculated since the 19th century (Darwin, 1859). In the last three decades, molecular techniques have become another valuable tool in biogeographic studies and have provided important insights into the origins and assembly of the world's biomes (Crisp, 2006). For example, evidence from the Northern Hemisphere, where most research has been conducted (Byrne et al., 2008), suggests that Pleistocene glacial oscillations caused the distribution of species to change considerably (Avise, 2000) and may have been responsible for the decline and extinction of many species (Guthrie, 2003). Those species surviving the glaciations responded by contracting into remnant habitat (refugia) during glacial maxima, and expanding during the interglacial periods that followed (Hewitt, 1996). This repeating pattern of range contractions and expansions, which mirrored the Pleistocene climate, has shaped the genetic structure of numerous species (Hewitt, 1996, 1999; Avise, 2000).

In the Southern Hemisphere, the Australian arid biome (defined by a moisture index of ≤ 0.4) has received relatively less phylogeographical analysis, and a number of questions relating to the diversification and the location of refugia for biodiversity in this biome remain unresolved (Byrne et al., 2008). The desertification of Australia began towards the end of the Miocene (c. 23–5 Ma) (Bowler et al., 2006), but the arid regions currently present in Australia have a much younger origin around the early Pliocene (c. 5–2.5 Ma) (Fujioka et al., 2005, 2009). Despite the recent origin, the arid biome is Australia's largest, in addition to being among the largest desert systems in the world, covering approximately 70% of mainland Australia (c. 7.5 million km2). Although the biome displays relative topographic homogeneity, habitat heterogeneity exists at both broader and smaller spatial scales. Species richness and ecological processes have been well described in arid Australia (Barker & Greenslade, 1982; Stafford Smith & Morton, 1990), and the region houses a rich and largely endemic lizard fauna (Cogger, 2000).

In Australia, the climatic oscillations during the Pleistocene did not result in substantial glacial and interglacial cycles, instead most regions experienced hyper aridity interspersed with more humid phases. This resulted in the formation of extensive dune fields during the hyper-arid periods, which created an inhospitable environment for most species due to the loss of water and vegetation. Some suggest that the waxing and waning between vegetated and denuded landscapes in Australia were analogous to movements of glacial ice in the Northern Hemisphere, whereby biota contracted into large-scale refugia (Hewitt, 2001; Byrne et al., 2008). However, the response varies among taxa, where both large-scale expansion and contraction events and more localized responses have been discovered (Byrne et al., 2008). As the climatic conditions in Australia continue to change, the location of past arid zone refugia are important areas to identify for conservation programmes. Potentially, these areas may be refugia again in the future. Furthermore, parts of a species' distribution in close proximity to refugia typically contain the highest levels of genetic diversity, and this diversity may be needed to maintain the evolutionary potential of species during periods of rapid change. Levels of genetic variation reflect effective population sizes and isolated parts of a species' distribution with few individuals, or recently founded areas, will have relatively low effective population sizes. Therefore, we predict that areas where a species has persisted over historical time-scales (refugia) will have consistently higher effective population sizes and as a consequence, higher levels of genetic variation.

In this study, we investigate the demographic history for the Gehyra variegata 2n = 40a chromosomal race (from here on referred to as G. variegata40a), a small gecko with an extensive distribution throughout the Australian arid zone. The species is commonly found within woodland habitats, underneath the exfoliating bark and woody crevices of its host trees, which it utilizes for thermoregulation, predator avoidance, egg incubation and foraging (Bustard, 1968). Due to this species dependence on vegetation (and indirectly the water that supports it), we hypothesize that the Pleistocene climate strongly influenced its distribution. We test whether levels of genetic variation are associated with environmental features in the Australian arid zone and whether lineage divergence in G. variegata40a corresponded with the Pleistocene period.

The inland MacDonnell ranges and mesic east have been proposed as potential refugia during the Pleistocene (Byrne et al., 2008). Additionally, the channel country, a large area (c. 300,000 km2) located in the central to eastern arid zone of Australia, may have also provided refugia for species during drier periods. This region receives the majority of its water via run-off, where rainfall in northern Queensland flows south–west until it reaches Lake Eyre in South Australia. This paleodrainage channel has been a permanent feature of the arid zone for at least 60,000 years (Magee & Miller, 1998). We predict that historical refugia will contain higher genetic variation and/or demographic signatures of population stability. The environmental features that we evaluate as potential refugia are the MacDonnell Ranges in central Australia, the mesic east and channel country. We also test whether water sources are associated with genetic variation at smaller spatial scales. We then use species distribution modelling (SDM) to assess whether conclusions drawn for our genetic data are congruent with the predicted distribution during the last glacial maxima (LGM) and to predict whether the historical refugia we identified are likely to represent future refugia for G. variegata40a.


Study species and sampling strategy

Gehyra variegata is a gecko with two distinct chromosome races (2n = 40a and 40b) and is widely distributed throughout the arid regions of Australia (King, 1979). The 2n = 40b race and the closely related Gehyra lazelli (Sistrom et al., 2009) are rock dwellers (Kitchener et al., 1988; Sarre, 1998). Here, we focus on the most widely distributed tree-dwelling G. variegata40a from the central and eastern arid regions of Australia, this race favours exfoliating bark and tree debris for diurnal retreat sites (Bustard, 1968; Henle, 1990; Duckett & Stow, 2012). We used morphological characteristics to avoid misidentification with G. lazelli, which is sympatric with the G. variegata40a in some areas (Sistrom et al., 2009). ND4 mitochondrial sequence data from each individual also confirmed the presence of a single species.

Our study sites were located throughout the arid and semi-arid zones of central and eastern Australia (Fig. 1), which closely matches the complete distribution of G. variegata40a (King, 1979). The arid zones were broadly categorized to allow patterns of genetic diversity to be associated with broad-scale environmental features across the species' distribution, including the inland ranges (greatest elevation), the mesic east (densest vegetation) and channel country (largest catchment). These features largely correspond to the MacDonnell Ranges, Brigalow Belt and Channel Country, respectively (Australian-Government, 2012). The woody vegetation in most of these landscapes are dominated by Acacia aneura (mulga) and Acacia cambagei (stinking gidgee or wattle), with the exception of the eastern regions where woodland is typically comprised of the genera Eucalyptus and Callitris. Between 2009 and 2011, we captured 740 individual geckos by searching potential retreat sites throughout its known distribution (King, 1979). Capture locations were recorded using GPS, and prior to release, a tail-tip tissue biopsy was taken and stored in 95% ethanol for analysis back in the laboratory.

Figure 1.

The mean annual run-off for regions within each of the Australian states and territories, and the location of sampling localities in this study (colour coded to their designated population). Numbers 1, 2 and 3 marks the location of the Inland ranges, Mesic east and the Simpson Desert, respectively. The highest levels of genetic variation were associated with the channel country where mean annual run-off is also higher (C3, 6, 9 and 12). C2, 4 and 7 represent likely suture zones. (Modified from source to standardize legends between states and territories). Please view the original source for a greater resolution within each state/territory – Australian Government Natural Resources Atlas, National Land and Water Resource Audit 2000 (www.anra.gov.au).

We obtained the mean annual surface water run-off statistics for multiple regions within each of the states and territories covered by this study from the Australian Government Natural Resources Atlas, National Land and Water Resource Audit 2000 (www.anra.gov.au) (Fig. 1). Historical measures were not available for Australia; therefore, contemporary measures may not reflect those of the Pleistocene period. However, a paleodrainage channel has persisted in this region for at least 60,000 years (Magee & Miller, 1998).

Laboratory procedures

Total DNA was extracted using a modified salting-out protocol (Sunnucks & Hales, 1996). All 740 individuals were genotyped at 16 microsatellite loci (Table 1), and sequences for the mtDNA ND4 region were obtained by polymerase chain reaction (PCR). The final reagent concentrations and the thermocycling conditions for PCR are outlined elsewhere (Forstner et al., 1995; Hoehn & Sarre, 2006; Duckett & Stow, 2010). To test our data quality, we regenotyped and sequenced five per cent of our sample. Individuals were randomly selected, genotyped at each microsatellite locus and resequenced at ND4. After ascribing genotypes and sequences to individuals, these data were then compared with the previously collected data to estimate error rates.

Table 1. Observed (HO) and expected (HE) heterozygosity for 16 microsatellites (Hoehn & Sarre, 2006; Duckett & Stow, 2010) across 740 individual Gehyra variegata40a samples
Locus N N A Range H E H O F IS
  1. FIS values marked* showed significant deviation from Hardy–Weinberg Equilibrium after adjustment for multiple tests (α < 0.05).

Mean 30 0.860.750.12*

Summary genetic data

For the microsatellite dataset, we calculated summary statistics including number of alleles (NA), observed heterozygosity (HO) and expected heterozygosity (HE) using genalex v6.0 (Peakall & Smouse, 2006). Measurements of FIS, linkage disequilibrium and the significance of any deviation from Hardy–Weinberg equilibrium (HWE) were calculated with fstat v2.9.3, correcting for multiple tests (Goudet, 2001). Estimates for null allele frequency for each locus was performed using the expectation maximization (EM) algorithm (Dempster et al., 1977) using FREENA (Chapius and Estoup, 2007). The sequence alignment and summary statistics for the mtDNA dataset were produced using mega v5 (Tamura et al., 2011) (Table 2).

Table 2. Summary statistics for 740 Gehyra variegata40a mtDNA ND4 sequences
G. variegata 40a x
Base pairs541
Haplotype diversity0.97
Conserved sites274
Variable sites267
Parsimony informative sites251
Singleton sites16
Alignment gaps/Missing data2
G+C/A+C content0.46/0.54
Transition/Transversion bias (R)5.16

Phylogenetic analyses

We used modeltest v3.7 (Posada & Crandall, 1998) to determine the best-fit model of DNA substitution among 88 tree construction models. Each model was assessed using the Akaike, Corrected Akaike, or Bayesian information criterion likelihood scores. Utilizing all sequences, both distance based and Bayesian methods were adopted to construct trees and compare the distribution of haplotypes without phylogenetic inferences. In mega v5, we used 10,000 bootstrap replicates to generate neighbour-joining and maximum parsimony trees. In MrBayes (Ronquist & Huelsenbeck, 2003), we used two runs, each with one cold and three heated chains for 2 × 107 generations, and with the sampling of likelihood parameters taken every 1 × 102 generation. The convergence of the chains was confirmed by analysis of the log likelihoods similarity, potential scale reduction factor (1.00 to < 1.01) and the standard deviation of split frequencies (≤ 0.01). The first 10% of trees were discarded as burn-in, and those remaining were used to construct a majority rule consensus tree. We identified and make reference to well-supported mtDNA clades that exhibited reciprocal monophyly (Moritz, 1994). Sequence divergence was quantified using mega v5, and consideration was given as to whether these clades represented individuals from discrete parts of the distribution (Ryder, 1986).

Divergence dating

It has become common practice in phylogenetics to date divergence events. To achieve absolute dating requires reliable calibration points from the fossil record (Ho & Phillips, 2009), which are rare for most datasets (Heads, 2005). However, it is widely accepted that reasonable approximations can be generated (Byrne et al., 2008). In this study, we estimate the divergence among lineages of G. variegata40a using beast v1.5.2 (Drummond & Rambaut, 2007). The calibration point for the time to most recent common ancestor (TMRCA) for all our samples was adopted from a phylogeny of 35 Gehyra species throughout Oceania, which dated species' divergences based on fossil record calibration and ND2 sequence data (Arnold & Poinar, 2008; Heinicke et al., 2011). For each clade, we used the mtDNA sequences for each unique haplotype and defined taxon sets that were reciprocally monophyletic based on the MrBayes phylogeny. We used the best-fit substitution model inferred from modeltest (General Time Reversible + Gamma/Invariant Sites), with an uncorrelated lognormal relaxed clock model with an estimate clock rate to allow for rate heterogeneity among lineages (Drummond et al., 2006). We used a coalescent tree prior, because all individuals were of the same species, and the choice of constant population size model was based on the results from Bayesian skyline plots (not shown). For the TMRCA prior, we specified a normal distribution with a mean ± SD of 8 ± 3 Ma, corresponding to the estimated time when G. variegata first diverged from its closest relative (Heinicke et al., 2011). Posterior estimates of parameters were obtained from two independent runs, each with 7.5 × 108 iterations with parameters taken every 1 × 103 iterations, with the first 20% discarded as burn-in. The convergence of the stationary distribution was assessed via analysing the posterior samples in tracer v1.5 to ensure all Effective Sample Size parameters were > 200 (Rambaut & Drummond, 2007). We consider the time estimates approximate, with our primary objective being to determine whether divergences correspond with the Pleistocene period.

Demographic history

We analysed mismatch distributions to infer the demographic history for each of the clades we identified by testing whether the mtDNA sequence data deviated from neutral expectations. Multimodal distributions suggest historically stable or amalgamated populations, in comparison with uni-modal populations that suggest expanding populations (Slater, 1987; Rogers & Harpending, 1992). We performed Tajima's D and Fu's Fs to test for recent directional selection (Tajima, 1989; Fu, 1997), and the McDonald–Kreitman test to test for natural selection (McDonald & Kreitman, 1991). Tajima's D is used to detect changes in population size or selection, where significantly negative values suggest expanding populations or positive selection, and positive values indicates contracting populations or balancing selection. Fu's Fs is used to detect population expansion or genetic hitchhiking. In both tests, where values approach zero, populations are inferred to be stable. We then tested for deviations from a sudden expansion model using Harpending's Raggedness Index, where significant values indicate population stability (Harpending, 1994). We also tested for evidence of recent admixture where a significant result may indicate past suture zones (Chakraborty, 1990). The significance of each test was assessed with 1 × 104 parametric bootstraps using arlequin v3.11 (Excoffier et al., 2005). To visualize the distribution and frequency of mtDNA haplotypes for each clade, we used tcs v1.21 (Clement et al., 2000).

Measures of genetic diversity

To quantify and compare the levels of genetic diversity within each clade, we used our microsatellite and mtDNA sequence datasets. Because larger sample sizes are expected to contain more genetic diversity, in each test, we account for differing sample sizes using rarefaction to the minimum number of individuals within a clade (N = 4). Standardized microsatellite diversity, referred to as allelic richness (AR), was calculated in fstat v2.9.3 (Goudet, 2001). Haplotype richness (HR), based on the inverse of Simpson's Index of Diversity, and Phylogenetic Distance (PD) was calculated in r v2.14.0 (R Core Development Team, 2011) with the Picante and Vegan libraries. To test whether higher genetic diversity was associated with broad-scale environmental features, for each clade, we estimated genetic diversity (PD, HR, and AR) and regressed these measures against the geographic distance from the landscape features being tested as potential refugia. To test any association between genetic variation and proximity to the MacDonnell ranges and the mesic east, geographic distances were measured from the central point of each clade to either the central point of the MacDonnell Ranges or the nearest point of contact with the mesic east bioregion. Mean annual run-off was also regressed against measures of genetic variation. Regressions were carried out using vassarstats (http://vassarstats.net). Our prediction was that genetic diversity will increase with proximity to locations that have acted as historical refugia.

Testing whether higher genetic diversity is associated with proximity to water

To test whether higher genetic diversity was associated with proximity to a water source over smaller distances (within a single ‘landscape’ type), we regressed standardized measures of genetic diversity generated from a total of 102 individuals with linear geographic distance from water sources in the Simpson Desert. The Simpson Desert covers a region with among the lowest average annual rainfall totals in Australia, and may approximate what were more widespread conditions during dryer phases of the Pleistocene. At this location, we sampled individuals in six separate areas (< 3 km2) to a maximum distance of c. 26 km from a single water source, Eyre Creek and Lake Nappanerica. Each of these locations contained a sufficient number of individuals to capture common alleles (> 5%). Regression and significance testing used vassarstats.

Species distribution modelling – past, present and future

If the highest levels of genetic diversity represent Pleistocene arid zone refugia for G. variegata40a, then these areas should be congruent with their predicted LGM distribution. We performed SDM with maxent 3.3.1 (Phillips et al., 2006; Phillips & Dudík, 2008). We modelled the likely present and past environmental suitability for G. variegata40a using a substantial occurrence record dataset, and eight climatic variables (Appendix S1) at a 2.5 arc min (c. 5 km2) resolution obtained from the WorldClim database (www.worldclim.org) (Hijimans et al., 2005). The climatic variables allow the characterization of the mean and variability of temperature and rainfall throughout the region and were selected over others due to their higher predictive ability in the model. The past distribution of G. variegata40a was projected to the LGM (c. 21,000 years before present) using the Paleoclimate Modelling Intercomparison Project Phase II (PMIP2) dataset. Future climate projections were sourced from Duckett et al. (2013). All predicted distributions were transformed into presence–absence maps. See Appendix S1 for a detailed description of the modelling procedure.


Summary statistics for the genetic data

The 16 primer pairs successfully amplified polymorphic loci with unambiguous alleles (Table 1), and revealed high levels of allelic variation, with NA = 8–54, HO = 0.20–0.86, HE = 0.24–0.95. Error rates estimated from repeat genotyping and sequencing, and checks of data entry, found errors at a low frequency (< 0.01), and these were corrected. All loci showed a significant homozygote excess across the complete dataset. However, analysis at smaller spatial scales (within 4 km2) rarely showed significant deviation from HWE. This indicates that the homozygote excess apparent when all samples are pooled is a likely Wahlund effect owing to genetic structure (Wahlund, 1928). We had a high rate of amplification success across all loci, and null alleles were at a low frequency (mean ± SD; 0.04 ± 0.02), and there was no evidence for linkage disequilibrium. Following the removal of ambiguously aligned nucleotides, 740 mtDNA ND4 sequences were trimmed to 541 bp, which revealed high levels of variation among haplotypes (Table 2).

Species phylogeny and estimated timing of divergence

All three model evaluation methods provided strong support for the generalized time-reversible mode of evolution with a gamma distribution and a proportion of invariant sites (GTR+G+I) (Tavare, 1986). All tree construction algorithms consistently clustered individuals on the same branches, and allowed the identification of 12 major clades, where a minimum of 2% sequence divergence between reciprocally monophyletic groups was found. In most cases, the clades contained individuals that were sampled from the same geographic region (Fig. 2). The basal divergences within our phylogeny pre-date the Pleistocene epoch (7,970,000–2,740,000 yr bp), whereas the majority of divergences among G. variegata40a lineages are congruent with a late Pleistocene influence (1,900,000–200,000 yr bp; Table 3).

Table 3. Summary statistics exploring the demographic history for each of the Gehyra variegata40a clades (C) we identified
CFu's FsP-valueTajima's DP-valueRaggednessP-valueChakraborty Obs/ExpP-valueTime to most recent common ancestor (TMRCA) (years)
  1. Values marked * showed significant deviation from neutral expectations (P < 0.02). Statistics are not provided for C1 due to a lack of haplotypic variation.

Figure 2.

A consensus Bayesian cladogram based on 541 bp of mtDNA ND4 sequence data. The 740 individual Gehyra variegata40a samples have been delineated into 12 major clades by reciprocal monophyly. Each line touching the cladogram edge represents a single sample and its designated clade. Note this type of tree only represents branching patterns and branch lengths do not represent time.

Contrasting demographic histories

Analysis of the number of synonymous and nonsynonymous substitutions in each clade detected a statistically significant excess of synonymous substitutions in clade 5 (P < 0.05), indicating a purifying selection on the ND4 gene. The assumption of neutrality could not be rejected in the other 11 clades. We infer contrasting demographic histories among clades (C1–12) (Table 3 and Appendix S2). Both clades 3 and 12, located within the channel country, displayed numerous haplotypes separated by relatively few mutational steps, while this might suggest recent population expansion this was not statistically supported with Fu's Fs or Tajimas D. Other channel country clades (C6 and 9) possessed haplotypes with a larger number of mutational steps between them relative to each of the other clades assessed. While these sorts of patterns are indicative of demographic stability, there was no statistical support for this. We found significant evidence for recent admixture in the south channel country, and the east (C2, 4, and 7, P ≤ 0.02). These suture zones are unlikely to represent past refugia and can be attributed to overlapping populations from past expansion events with subsequent isolation. We did not infer the demographic history for C1 due to the presence of only a single haplotype.

Higher genetic diversity is associated with the channel country

Our analysis of the microsatellite and mtDNA sequence datasets within C1–12 revealed that standardized measures of genetic variation were generally higher within C3, 6, 9 and 12 (Table 4). We show that the decline in mtDNA genetic diversity (PD and HR) is not strongly associated with increasing distance from either the ‘mesic east’ (R2 = 0.008 and 0.016; P = 0.388 and 0.345) or ‘inland ranges’ (R2 = 0.058 and 0.046; P = 0.223 and 0.249). Yet a weak association was inferred from the AR dataset (R2 = 0.273 and 0.245; P = 0.040 and 0.051). Here, the association visible in AR suggest increases in genetic diversity with proximity to the inland ranges, although this is largely driven by one point close to the mesic east. The clades that displayed the highest levels of genetic variation were located within the channel country of arid zone Australia, which is characterized by a higher mean annual run-off than any other sampling locality (Fig. 1). Across all 12 clades, we found a strong association between higher levels of mtDNA genetic diversity (PD and HR) and higher mean annual run-off (R2 = 0.459 and 0.546; P = 0.007, 0.003); a weaker association remained for AR, respectively (R2 = 0.230; P = 0.057) (Fig. 3).

Table 4. Standardized measures of genetic diversity for Gehyra variegata40a include allelic richness (AR), haplotype richness (HR) and phylogenetic diversity (PD) across each of the 12 clades (C)
Figure 3.

Regressions of standardized measures of genetic variation (PD, HR and AR) with the geographic distance (km) of collected samples from the Mesic East and the MacDonnell Ranges, and with Mean annual run-off (gigalitres per year). The strength (R2) and significance (P) of each regression is provided, a ‘trend’ line was included when the relationship was significant.

Higher genetic diversity is associated with localized water sources

The analysis of our Simpson Desert dataset (n = 102) shows that genetic diversity (PD, HR and AR) declines significantly with increasing distance from water sources in this region (Eyres Creek and Lake Nappanerica) (R2 = 0.786, 0.623 and 0.641; P = 0.009, 0.031 and 0.027, respectively) (Fig. 4). Within each of these areas, there is a visible increase in the stand size and frequency of Acacia species, which provides habitat for G. variegata40a. This suggests refugia may be associated with stable water sources, which facilitates the persistence of its habitat.

Figure 4.

Regressions of standardized measures of genetic variation with the geographic distance of collected samples from stable water sources found within the Simpson Desert. The strength (R2) and significance (P) of each regression are provided.

Species distribution modelling

Species distribution modelling revealed the likely distribution of those areas with favourable climatic conditions for G. variegata40a during the LGM, present and future. Modelling for past and present climates suggests the channel country has long been part of this species' distribution, which provides support for our genetic data. However, future projections suggest that favourable climatic conditions will not be retained in this region, and will contract and shift in an approximate south-east direction (Fig. 5).

Figure 5.

Presence–absence maps for the predicted distributions of Gehyra variegata40a. Favourable climatic conditions areas are shown in dark grey and unsuitable areas in light grey. The following distributions are presented: (a) the past [last glacial maxima (LGM)]; (b) the present (2010); (c) the future (2070 A2 scenario – using realistic population and CO2 increases and economic growth is favoured over environmental concerns). Each coloured dot represents the present location of sampled individuals and the clade they belong to. The black circle represents the approximate location of the channel country.


Phylogeographic approaches and SDM have rarely been applied to identify historical refugia in arid Australia. Genetic data from G. variegata40a sampled from several arid Australian landscapes and SDM indicate long-term persistence of this species in the channel country, a region that indirectly receives more water than other arid regions within Australia. While the channel country is inferred to have been part of the distribution of this species since at least the LGM and may have acted as a refugia throughout the Pleistocene, SDM for G. variegata40a did not predict favourable climatic conditions in the channel country in the near future.

Our data show broad-scale geographic structuring of lineages, where the estimated time of genetic divergences for the majority of clades were congruent with a peak in the amplitude for Pleistocene climatic oscillations; c. 600,000–200,000 yr bp (Byrne et al., 2008). The Pleistocene period has been associated with lineage splitting for marsupial, acacia and other lizard species in arid Australia (Blacket et al., 2001; Byrne et al., 2001; Chapple et al., 2004). The survival of G. variegata40a throughout hyper-arid Pleistocene events suggests the presence of refugia. These refugia could have been on the periphery of the expanding central Australian deserts. Thus, the species' current widespread distribution may reflect a major recolonization event since the (LGM). However, differences in levels of haplotype richness and, potentially, demographic histories across the contemporary distribution of G. variegata40a point to the presence of multiple localized refugia in the past.

Higher levels of genetic diversity were inferred for clades located in the channel country, in comparison with those clades from the inland ranges or mesic east. This provided some evidence for broader-scale refugia during the Pleistocene in channel country region. In the clades of the ‘mesic east’ and the ‘inland ranges’, potential locations for refugia (Byrne et al., 2008), we found comparatively lower levels of genetic diversity. This suggests that Pleistocene refugia might have been less extensive in these regions compared with the channel country for G. variegata40a. The channel country region receives water from monsoonal rains falling in the north of Queensland which then drain throughout (and supply) several hundred thousand square kilometres via an ancient paleodrainage channel towards Lake Eyre in the north-east of South Australia (Fig. 1). Outside of this region rainfall predominately flows towards the coast, and sporadic and localized rains supply much of the arid region. Thus, it would seem that local water availability and not local rainfall influences persistence. The importance of a stable water supply to G. variegata40a is also demonstrated by our genetic data at shorter spatial scales (Fig. 4). Here, we found that higher genetic diversity was associated with water presence, suggesting longer-term persistence occurs in closer proximity of water. Although again, the water sources are rarely filled by local rains and receive the vast majority of their water via rainfall in northern locations. During very arid phases, suitable vegetation for G. variegata40a is likely to be lost in areas other than those with a stable water supply. These data support the notion that water supply is likely to be associated with historical refugia.

The SDM performed in this study still made several commonly made assumptions which are yet to be fully evaluated in this area of research. For example, in some situations, assuming distribution equilibrium has been shown to result to inaccurately predict a species' range (Elith et al., 2010); however, this is more likely to be problematic for invasive species in a state of spread rather than the widely distributed and native G. variegata40a. We also assumed that in-situ behavioural and physiological responses to rapid climatic changes were unlikely. There may be some justification for this because the localized extinction of many lizard species has recently been documented (Sinervo et al. 2010). Our work also highlights the importance of considering the influence of climate outside the regions being evaluated by SDM. For G. variegata40a, water availability (run-off) may be an important predictor of species' presence and not localized rainfall per se, which currently is not explicitly captured in the models. This is relevant for many arid regions throughout the world where localized rains are infrequent yet water delivered by run-off is more common. For example, the Nile, Colarado and Yellow River all flow through desert ecosystems which themselves experience little rainfall. As climates continue to warm globally, incorporating run-off variables will likely improve the predictive ability of SDM for many species, especially given that climate models predict widespread decreases in water run-off throughout many countries in the Americas, Africa, Eurasia and Australasia (Arnell et al., 2011).

The phylogeography of G. variegata40a revealed that past refugia were associated with water supply and not necessarily where rains fall, and this is potentially the case for many other arid species throughout the world. The long-term persistence of stable water supplies in the channel country appear to have provided refugia for several aquatic species (Hamilton et al., 2005), and our data suggest that these regions may also be important for the persistence of woodland biodiversity. However, modelling of future climate suggests that the channel country region in Australia may not be suitable for our woodland species by 2070. Arid ecosystems are likely to change with the prediction of increasing temperatures and evaporation. Greater knowledge of how water run-off has influenced persistence throughout catchments in arid regions may improve our ability to predict the impacts of climate change and to identify and protect contemporary reservoirs of genetic variation.


This work was supported by Macquarie University. Lizards were handled in accordance with Macquarie University and South Australian Animal Ethics committee recommendations. Tissue collection was licensed by the New South Wales National Parks and Wildlife Service (S12353), Northern Territory Parks and Wildlife Commission (33155), Queensland Environment Protection Agency (WITK04616407) and South Australia Department of the Environment and Natural Resources (Z25499-4). We thank Macquarie University volunteers for their assistance during field collections.


Paul E. Duckett has research interests in marine ecology and biogeography with reference to the potential impacts of climate change. He is primarily interested in historical refugia and whether or not their location will be beneficial for the future of global biodiversity.

Adam Stow is a Senior Lecturer at Macquarie University, Sydney. With research interests in molecular ecology and conservation biology his primary focus is assessing gene flow, dispersal and genetic variability in human impacted environments.

Author contributions: P.E.D and A.S. conceived the ideas; P.E.D. collected the data; P.E.D analysed the data and P.E.D led the writing.