Genetic Divergence among Regions Containing the Vulnerable Great Desert Skink (Liopholis kintorei) in the Australian Arid Zone

Knowledge of genetic structure and patterns of connectivity is valuable for implementation of effective conservation management. The arid zone of Australia contains a rich biodiversity, however this has come under threat due to activities such as altered fire regimes, grazing and the introduction of feral herbivores and predators. Suitable habitats for many species can be separated by vast distances, and despite an apparent lack of current geographical barriers to dispersal, habitat specialisation, which is exhibited by many desert species, may limit connectivity throughout this expansive region. We characterised the genetic structure and differentiation of the great desert skink (Liopholis kintorei), which has a patchy, but widespread distribution in the western region of the Australian arid zone. As a species of cultural importance to local Aboriginal groups and nationally listed as Vulnerable, it is a conservation priority for numerous land managers in central Australia. Analysis of mitochondrial ND4 sequence data and ten nuclear microsatellite loci across six sampling localities through the distribution of L. kintorei revealed considerable differentiation among sites, with mitochondrial FST and microsatellite F′ST ranging from 0.047-0.938 and 0.257-0.440, respectively. The extent of differentiation suggests three main regions that should be managed separately, in particular the southeastern locality of Uluru. Current genetic delineation of these regions should be maintained if future intervention such as translocation or captive breeding is to be undertaken.


Introduction
The Australian arid zone occupies 70% of the continent's landmass and supports an extraordinary biodiversity, including among the world's richest assemblages of lizards [1], [2]. Despite a longstanding recognition of the conservation value of this region, relatively few studies have described patterns of genetic structuring across whole species distributions [3]. Characterisation of genetic structure across a landscape is valuable to inform conservation because low dispersal of the species. Here we use mtDNA sequence data (ND4) and ten microsatellite loci to characterise genetic structure and divergence across the range of L. kintorei. We aim to identify any regions where the genetic distinctiveness of L. kintorei heightens the conservation value and influences the management options.

Ethics Statement
All animals were handled in accordance with a protocol considered and approved by Macquarie and Charles Darwin University Animal Ethics committee recommendations (ARA 2008/025 and ARA 2011/037). Sample collection was licensed by the Northern Territory Parks and Wildlife Commission and the South Australian Department of Environment and Heritage.

Laboratory procedures
Whole genomic DNA was extracted from tissue using a modified salting-out protocol [26]. For each sample the mitochondrial gene, NADH dehydrogenase subunit 4 (ND4), was targeted because previous work on reptiles has shown useful levels of variation for intraspecific studies, though its mutation rate is slow enough to allow inference of deeper divergence due to longterm isolation [12,27]. In addition, individuals were genotyped at ten microsatellite loci to characterise fine-scale population dynamics and more recent divergence.
Microsatellite PCRs were carried out in 10 μL volumes containing~50 ng of DNA. A -29 M13 sequence was added to the 5' end of each forward primer to allow for the incorporation of a complimentary M13 fluorescent-labelled tag, following the protocol of Schuelke [29]. Four tetranucleotide microsatellite loci were developed concurrent to this study (BX6, CKD, FQR, J3F; see S1 Text for microsatellite design). In addition to these, we utilised six previously developed markers (Est1, Est2, Est9, Est12 [30]; Ecu2, Ecu3 [31]). All microsatellite loci were amplified with identical reaction conditions: 2 uL 5x GoTaq Flexi Buffer (Promega), 2.5 mM MgCl 2 , 0.2 μM of each dNTP, 0.02 μM of forward primer, 0.1 μM reverse primer, 0.1 μM of fluoro-labelled tag (FAM, VIC, NED, or PET) and 1 U Taq Polymerase (Promega). Thermocycling began with an initial denaturation for 3 min at 94°C, followed by five touchdown cycles with 94°C denaturation for 30 sec, annealing temperatures (60°C, 58°C, 56°C, 54°C, 52°C) for 30 sec, and 72°C extension for 45 s. An additional 35 cycles were carried out at an annealing temperature of 50°C, followed by a final 72°C extension step for 10 min. PCR products were visualized by electrophoresis on 2% agarose gel. All PCR purification, sequencing and fragment separation was performed by Macrogen (Korea).

Data analysis
ND4 sequences were checked by eye and aligned with ClustalW, implemented in MEGA 5.0 [32], and submitted to GenBank (Accession numbers KM035773-KM035789). DNA sequences were then translated into amino acid sequences using the vertebrate mitochondrial code. No premature stop codons were observed, indicating that all sequences are true mitochondrial copies. Haplotype and nucleotide diversities were calculated in DnaSP [33].
Microsatellite alleles were visualised and scored using Peak Scanner 1.0 (Applied Biosystems). To ensure amplification and scoring consistency, at least 10% of samples at each locus were independently rerun and genotyped. Summary statistics, including exact tests for Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were conducted in GenAlEx 6.4 [37] and GENEPOP 4.2 [38]. Effective population size (Ne) estimates were calculated utilizing the approximate Bayesian framework implemented in ONeSAMP v1.2 [39]. Due to prohibitively small sample sizes at Docker River and Warburton, these sampling localities were excluded from population-level analyses.
When calculating F ST analogues from highly polymorphic data such as microsatellites, within-population variance can often approach the level of the total variance, resulting in very low F ST values even when the populations share no alleles [40,41]. Following Hedrick [40] and Meirmans [41], pairwise fixation index values calculated from microsatellite data (hereafter F 0 ST ) were standardised using the program RECODEDATA 0.1 [41]. STRUCTURE v.2.3 [42] analysis was used to assess genotypic clustering and assignment probabilities. We examined values of K = 1-8 (double the number of sample sites included in the analysis), with 10 replicate runs for each, 10 5 MCMC iterations burn-in and 10 4 main iterations. Hubisz et al. [43] developed a new model for STRUCTURE, which allows the use of sample-site information. This is different to the initial models including location priors, in that it adds power to analyses, but can disregard site information when true clustering is uncorrelated with sampling locations. We used the 'admixture' model with correlated allele frequencies, and repetitions were run with and without location information. The number of genetic clusters (K) was determined using the ΔK method of Evanno et al. [44].
Discriminant analysis of principal components (DAPC) was used to describe the genetic relationship between sampling localities. DAPC is a multivariate analysis that first uses principal components analysis (PCA) to transform data into uncorrelated components. These components are then analysed using a linear discriminant method, minimising within-group variance while maximising among-group variance [45]. Furthermore, this analysis does not assume HWE and LD, which are often violated when working with natural, small and fragmented populations [45].
DAPC was carried out in the R package adegenet [46], implemented in R 2.12 (R development core team 2013; www.r-project.org), with K selected using the find.clusters function and Bayesian Information Criterion (BIC). We also ran DAPC using sample locations as groups (K = 4) to assess the differentiation of our sample sites. PCA was performed in R using the dudi. pca function in the package ade4 [47]. Missing data were replaced with the mean (the origin of the X-and Y-axes, as in Horne et al. [48]). Determining the number of principal components (PCs) to retain as predictors for the discriminant analysis requires a balance between the statistical power of more PCs, and the stability of assignments, though there is no strict rule. Retaining too many PCs with respect to sample size can result in over-fitting the data. This trade-off can be assessed using the a.score function in the R package adegenet [46]. Analyses were carried out retaining a conservative 13 PCs, the optimal number suggested by a.score, given our relatively small dataset.

Summary statistics
Mitochondrial sequence data. Mitochondrial ND4 sequences of 585 bp were successfully amplified from 72 individuals sampled from the six localities. Sequences contained 23 (3.9%) variable sites, of which 18 (3.1%) were parsimony informative, revealing a total of 17 unique haplotypes (Fig 2). Haplotype and nucleotide diversity over all samples were 0.908 and 0.009, respectively, and the variance for both was < 0.0002 (Table 1). All haplotypes sampled at Uluru were unique to that locality. Newhaven and Sangster's Bore shared haplotypes with only each other, and Docker River and Warburton both shared haplotypes with Watarru. One sample from Warburton was unique to that locality, although small sample sizes may be the reason for this.
Microsatellite loci. One of the four remaining sample localities was out of HWE (Newhaven, P < 0.05; Table 2), which may be due to a Wahlund effect from some spatial structure [49]. Following Holm-Bonferroni sequential correction [50], 10 out of 180 locus x locus tests for LD (45 per sampling locality) were significant, all of which were for different locus pairs. The presence of some LD is unsurprising, given that it can be common in threatened species  Genetic differentiation between localities. All comparisons of genetic differentiation, except one (Newhaven-Sangster's Bore; P ND4 = 0.108), were high and significant (Table 3), indicating very low connectivity between localities. For ND4, overall F ST = 0.50, P < 0.00001. Pairwise population differentiation for both ND4 (F ST ) and microsatellites (F 0 ST ) was substantial (Table 3; Juke's Cantor distances between localities are given in Table 4). Newhaven and Sangster's Bore were the least differentiated from each other, with low and not-significant F ST for ND4, though microsatellite differentiation was relatively high. Uluru was the most differentiated from all other localities in all comparisons. STRUCTURE analysis yielded a best-fit value of K = 2 (S1 Fig) without location information, and K = 4 with location information. When K = 2 was considered, one cluster comprised the Uluru samples, and the other clumped Newhaven, Sangster's Bore and Watarru together (S1 Fig). When location information was used, STRUCTURE gave a best-fit value of K = 4 (Fig  3a), and there was clear delineation of sample sites with high assignment probabilities (Fig 3b).
A similar result was found in the DAPC: the lowest BIC value was indeed K = 2 (BIC = 124) when no location information was used, however the BIC score for K = 4 was only slightly  higher (BIC = 126; a difference in BIC of 0-2 is considered weak [51]). For K = 2, Uluru samples were largely separate from the rest (S1 Fig). However, when samples were grouped by location in the DAPC, there was considerable genetic separation between all four groups, with the largest separation along the first and second discriminant factors lying between Uluru and the other localities (Fig 3c). The optimal number of principal components retained for the DAPC accounted for 54% of the overall variation, and the first two discriminant factors accounted for 91% of this variation in the discriminant analysis. Only plots for K = 4 are given here, but see S1 Fig for plots of K = 2.

Discussion
We show genetic partitioning among regions containing the Vulnerable lizard Liopholis kintorei. Each of the localities from which L. kintorei was sampled contained similar levels of genetic variation, and individuals from the Uluru-Kata-Tjuta region were most genetically distinctive (Tables 1, 3 and 4). Our estimates of genetic divergence, in addition to environmental differences experienced among regions, indicate that each of these are reservoirs of important genetic variation and point to the risk of outbreeding depression should interbreeding occur.
Outbreeding depression is lowered reproductive fitness in generations subsequent to crossing of individuals from genetically differentiated parts of their distribution [16]. A decision tree developed by Frankham et al. [16] assists conservation managers to assess the probability of its occurrence, and thus decide whether populations should be kept separate. Applying this to the data for Liopholis kintorei indicates a risk of outbreeding depression should individuals be translocated. The third part of the tree suggests that if sites have been isolated from each other for 500 years or more, there is a high risk of outbreeding depression and they should remain separated. A crude estimate of divergence times based on a commonly cited mitochondrial calibration of 1.3-2% sequence divergence per million years [13,52], suggests that the Uluru lineage may have split from the others between 350 kya and 1.31 million years ago. The level of genetic divergence between sampling localities at ND4 (< 2%; Table 4) was below that found in a closely related species, Liopholis inornata (within-species divergence up to 6.1% [12]); however, L. inornata has a much broader distribution and was sampled across a wider area, that may explain the higher within-species divergence reported. When considering clades of L. inornata sampled across similar geographic scales, the level of differentiation within these two Liopholis species was similar.
Furthermore, there are environmental differences between the sites that may contribute to localised adaptation, another factor that flags the possibility of outbreeding depression from translocation according to Frankham et al. [16]. Similarly, Crandall et al. [15] propose evaluating 'ecological exchangeability' along with estimates of genetic divergence to decide on the parts of species distributions to be managed separately. Given the latitudinal range over which L. kintorei is distributed, it is not surprising that there are environmental differences across our sampling localities. For example, mean monthly minimum and maximum temperatures are consistently and significantly cooler at Uluru than the northernmost locality of Sangster's Bore (see S2 Text for climate data and statistical tests). Average annual rainfall at Sangster's Bore is 479 mm compared with 320 mm at Uluru (t 32 = 2.09, P = 0.045), though rainfall patterns differ with Uluru receiving more winter rain. Sangster's Bore experiences more days above 35°C and 40°C than Uluru (respectively, 169.6 vs. 109.3 days 35°C and 52.9 vs. 32.1 days 40°C; climate data from Australian Bureau of Meteorology, www.bom.gov.au; data not available for Watarru). Given their ectothermic physiology, life history traits in reptiles have been demonstrated to be linked to altitudinal or latitudinal variability in climate [52][53][54]. As such, the significant climatic differences between our sampling localities are likely to be relevant to L. kintorei, and may have led to adaptive differences in, for example, thermal tolerance or seasonal activity.
The habitat in which L. kintorei is found in these areas varies also. The Sangster's Bore and Newhaven localities occur along and adjacent to palaeodrainage lines, with the species' preferred habitat in semi-saline spinifex plains dominated by soft spinifex (Triodia pungens) and inland tea tree (Melaleuca glomerata). At Watarru, L. kintorei were found within open mulga woodland or Eremophila and woolybutt (Eragrostis eriopoda) grass shrubland, and at Uluru they occur on sand plains and flat swales dominated by either hard (Triodia basedowi) or soft spinifex (T. pungens and T. schinzii) [21].
Uluru was the most genetically distinct in all analyses, and shared no haplotypes with any other site; Watarru was also highly differentiated from the other areas, but shared haplotypes with the other two southwestern localities (Docker River and Warburton). More sampling at Docker River and Warburton is required to determine the true extent of differentiation among these localities. Sangster's Bore and Newhaven to the north were highly differentiated based on the microsatellite data set (F 0 ST = 0.285), but were not significantly differentiated at mitochondrial ND4 (F ST = 0.047). This discrepancy, taking into account the different inheritance and mutation rates for these genetic markers [55] might indicate a contemporary barrier to dispersal but higher levels of historical gene flow. As a result, we recognise three main delineations among localities for conservation management: one to the north (Sangster's Bore and Newhaven), one to the southeast (Uluru), and one in the southwest (Watarru, Docker River, Warburton).
Given the isolation of localities implied by an apparently patchy distribution [20,22] and the genetic differentiation among localities investigated here, genetic diversity, if lost, may not be replenished by migration. Effective population sizes are substantially lower than actual sizes in wild populations, with the ratio between them often approximating 0.1 [56], and up to 0.5 [57]. Census size estimates for L. kintorei are estimated to be low (<500 at Uluru [20]), and genetic Ne estimates in this study were all low (138-232, though these estimates should be treated with caution due to small sample sizes). The Vulnerable status of this species and low estimated population sizes suggest that genetic diversity and viability may be eroded rapidly over time. The threatening processes attributed to the decline of great desert skinks have not been removed, and consequently this is likely to continue. If genetic erosion is allowed to proceed, this can render localised parts of the distribution vulnerable to inbreeding and inbreeding depression [58]. Translocations to bring about a so-called 'genetic rescue' have been demonstrated to dramatically reverse the effects of inbreeding depression [17,59]. In the case of L. kintorei, unique haplotypes at some localities and high differentiation estimates suggest that if this need eventuates, parts of the distribution selected for translocations need to be carefully chosen to avoid the risk of outbreeding depression.
Knowledge of connectivity combined with landscape management of biological processes is needed to conserve biodiversity [15]. Conservation management for L. kintorei should prioritise the preservation of suitable habitat, in particular addressing recent and localised changes in fire regimes and predation pressure to reduce the risk of further localised population declines and thus erosion of genetic diversity. While further sampling needs to be conducted at Watarru, Docker River and Warburton, the evidence suggests three main delineations for management: (1) Uluru to the southeast, (2) Newhaven and Sangster's Bore to the north, and (3) Watarru, Docker River and Warburton to the southwest. Uluru in particular should be considered separately for management, and this distinctiveness should be recognised if intervention such as translocation or captive breeding is to be undertaken.