Persistence and dispersal in a Southern Hemisphere glaciated landscape: the phylogeography of the spotted snow skink (Niveoscincus ocellatus) in Tasmania

The aim of this research was to identify the effects of Pleistocene climate change on the distribution of fauna in Tasmania, and contrast this with biotic responses in other temperate regions in the Northern and Southern Hemisphere that experienced glacial activity during this epoch. This was achieved by examining the phylogeographic patterns in a widely distributed Tasmanian endemic reptile, Niveoscincus ocellatus. 204 individuals from 29 populations across the distributional range of N. ocellatus were surveyed for variation at two mitochondrial genes (ND2, ND4), and two nuclear genes (β-globin, RPS8). Phylogenetic relationships were reconstructed using a range of methods (maximum parsimony, Bayesian inference and haplotype networks), and the demographic histories of populations were assessed (AMOVA, Tajima’s D, Fu’s Fs, mismatch distributions, extended Bayesian skyline plots, and relaxed random walk analyses). There was a high degree of mitochondrial haplotype diversity (96 unique haplotypes) and phylogeographic structure, where spatially distinct groups were associated with Tasmania’s Northeast and a large area covering Southeast and Central Tasmania. Phylogeographic structure was also present within each major group, but the degree varied regionally, being highest in the Northeast. Only the Southeastern group had a signature of demographic expansion, occurring during the Pleistocene but post-dating the Last Glacial Maximum. In contrast, nuclear DNA had low levels of variation and a lack of phylogeographic structure, and further loci should be surveyed to corroborate the mitochondrial inferences. The phylogeographic patterns of N. ocellatus indicate Pleistocene range and demographic expansion in N. ocellatus, particularly in the Southeast and Central areas of Tasmania. Expansion in Central and Southeastern areas appears to have been more recent in both demographic and spatial contexts, than in Northeast Tasmania, which is consistent with inferences for other taxa of greater stability and persistence in Northeast Tasmania during the Last Glacial Maximum. These phylogeographic patterns indicate contrasting demographic histories of populations in close proximity to areas directly affected by glaciers in the Southern Hemisphere during the LGM.


Background
The glacial-interglacial oscillations of the Pleistocene have strongly influenced the distribution and evolution of many species [1]. Initial interest in the effects of these oscillations largely focused on temperate species within Europe [2]. The classic phylogeographic view is that these climatic changes forced a diverse range of temperate European taxa to become geographically restricted to a series of spatially discontinuous, low latitude, Mediterranean glacial refugia [1,3,4]. During glacial intervals these refugia accumulated genetically divergent lineages, and later provided the source populations for interglacial range expansions. As recolonising populations are expected to represent only a subset of the diversity associated with glacial refugia (due to founder events and high density blocking, see [5]), phylogeographic methods can locate refugia by assessing patterns of genetic diversity across a species' geographical range [2].
Recently, an increasing number of phylogeographic studies have questioned the simplicity and universality of this classical southern Europe refugia model, describing 'cryptic' northern refugia [6][7][8] and 'refugia within refugia' [8,9] for a range of temperate European species. However, withstanding this growing complexity several facts remain clear-for a large proportion of taxa the LGM had a strong biogeographic influence, and that for most species a limited number of large southern refugia supported the majority of genetic diversity throughout the Pleistocene. While such generalisations can be made for Europe and other Northern Hemisphere regions, the impact of Pleistocene climate change on biota in the Southern Hemisphere remains comparatively understudied [10][11][12][13][14][15].
During glacial periods, ice was absent or discontinuous across most of the Southern Hemisphere continents with the exception of Antarctica. This has meant that many austral studies have focused on the effects of Pleistocene changes in aridity, in tropical or ice-free temperate regions, rather than the effects of glaciers directly (e.g., [16][17][18][19]). Furthermore, Southern Hemisphere studies of the impacts of glacial and periglacial activity on species distributions have concentrated on New Zealand, Patagonia, and Antarctica (e.g., [20][21][22][23][24]), while Tasmania-the focus of Australia's most extensive Pleistocene glacial activity [25]has been comparatively neglected. Yet, unlike Patagonia and New Zealand, Tasmania has not experienced recent tectonic activity, which has the ability to confound interpretations regarding the influence of glaciations on species distributions and gene flow [22,24]. Although not presently glaciated, Tasmania possessed glaciers at least five times during the Pleistocene, and these glaciers were most extensive (covering up to 7000 km 2 ) during the early (~1.8 Myr) and middle Pleistocene (>130 kyr) [26][27][28][29][30] (Fig. 1). During the LGM, smaller ice caps formed at high altitudes, and local summer temperatures are estimated to have been 6-8°C cooler than present averages [25,28,31] (Fig. 1).
Plant studies indicate a wide range of phylogeographic patterns in Tasmania, but those of several Eucalypt species [32] have led to the suggestion that species responses to Pleistocene glacial activity in Tasmania may be similar to classical European patterns-whereby species retreated into a few large and distinct refugia during the LGM, and the genetic diversity of recolonising populations is limited [33]. However, these patterns contrast strongly with those typically observed in other parts of Australia and the Southern Hemisphere, which have suggested greater numbers of refugia, including many micro-refugia (small pockets of suitable habitat in an otherwise uninhabitable landscape) [34], and genetic structuring reflecting events which pre-date the LGM [17,18,22,24,35]. Consequently, Tasmania represents an important region for further research to better develop our understanding of the effects of Pleistocene glaciations on species distributions [36]. In this context, lizards have recently been highlighted as good models for use in phylogeographic studies [37]. As ectotherms, lizards are sensitive to changes in climate, which may be manifested in alterations of species distributions [38]. Lizards also typically have low mobility, and phylogeographic patterns from historic events will be retained for longer periods [39].
This study investigated the phylogeography of a Tasmanian endemic reptile, the spotted snow skink (Niveoscincus ocellatus). This species is an excellent model for investigating the impacts of Pleistocene glacial cycles in Tasmania for two reasons: it has a wide geographic distribution across Tasmania, and it has clear distribution restrictions associated with climate and habitat type. Niveoscincus ocellatus ranges from sea level to high elevations, including previously glaciated regions ( Fig. 1) [40,41], but it is not currently found at altitudes abovẽ 1200 m, beyond which it is sharply replaced by biennially reproducing, alpine specialist species: N. greeni and N. microlepidotus [42,43]. Mean temperature for the warmest month is presently~10°C at 1200 m elevation, yet the corresponding temperature would have been~4°C during the LGM [44]. It is expected that the historical distribution of N. ocellatus will have been strongly regulated by Pleistocene climate change, such that there will be genetic signatures of Pleistocene refugia and subsequent range expansion. Furthermore, N. ocellatus only occurs amongst rocky outcrops [45], which suggests that gene flow will be generally low [46] (but see [47,48]) and genetic signatures of historic distributions could be well preserved.

Mitochondrial DNA sequence variation and phylogeographic relationships
A total of 95 unique haplotypes were identified from the 204 N. ocellatus sampled, with 183 (12.9 %) of 1420 characters variable and 147 (10.4 %) parsimony-informative.  [28]. Black lines demarcate groups defined on the basis of mtDNA variation, with the dashed black line indicating the ambiguous placement of the Lake Mackenzie population either within the Southeast group or on its own. Elevation data were obtained from Geoscience Australia Tree topologies were very similar between maximum parsimony and Bayesian tree-building methods, only differing in relationships at shallow phylogenetic levels, while the composition of major clades were consistent across analyses; therefore only the Bayesian inference tree is presented ( Fig. 2; parsimony tree Additional file 1: Figure S1). Four major clades are evident: a 'Northeastern' clade, a 'Southeastern' clade, a 'Northwestern' clade, and a clade of three Lake Mackenzie haplotypes (Fig. 2). The relationships among these four clades were uncertain and received low topological support, but with the exception of the three Lake Mackenzie haplotypes each clade received posterior probability greater than 0.95. The topology obtained from BEAST analysis using a coalescent tree prior and a strict, externally calibrated, molecular clock was compatible with that described above, and suggests that Fig. 2 Bayesian inference tree for Niveoscincus ocellatus based on 1420 bp of ND2 and ND4 mitochondrial DNA sequence. Branch lengths are scaled relative to the scale bar. Posterior probabilities >0.95 are indicated by the black dots, with values at key nodes also specified above branches at nodes. Bootstrap values from parsimony analysis, where they exceeded 70 %, are also listed below branches. Numbers in parentheses indicate the number of individuals from a site exhibiting that haplotype. The branch leading to the outgroup N. greeni was truncated to aid visualisation, and likewise that leading to N. pretiosus was removed. Grey boxes highlight regionally monophyletic areas within the major clades the major lineages diverged during the last 2 Myr (Additional file 1: Figure S2).
SAMOVA recovered a group corresponding to the Northeastern clade when testing for two populations (Φ CT = 0.50, P = 0.00), and when testing for three populations recovered groups similar to that inferred by the phylogeny (Φ CT = 0.53, P = 0.00): Northeastern, Southeastern (but with inclusion of both Lake McKenzie and Dove Lake), and Northwestern (Mt. Oakleigh only). Given no appreciable increase in Φ CT with increasing number of populations in SAMOVA, three regional groups were defined for further population genetic analyses. The Northeastern and Northwestern (Mt Oakleigh and Dove Lake) clades each corresponded to respective groups, but the Lake Mackenzie clade (three haplotypes) was combined into the same group as the Southeastern clade (hereafter, the Southeastern group) as two of the Lake Mackenzie haplotypes clustered monophyletically among Southeastern haplotypes, and this locality is geographically confluent with Southeastern clade localities. The same three groups are readily distinguished by large numbers (>19) of mutations in the TCS network (Additional file 1: Figure S3). With respect to these groups, AMOVA indicated that spatial structuring of mitochondrial genetic variation was significant at each hierarchical level-among populations, among populations within groups, and among groups ( Table 1). The Northeastern group had the highest degree of spatial structuring, where most localities (excluding Ben Lomond) were individually monophyletic (Fig. 2). Structuring was lower within the Southeastern group, where only a few populations of contiguous regions formed monophyletic clades (Table 1; Fig. 2).
Haplotype diversity was higher in the Southeastern group than the Northeast (Table 2). In contrast, nucleotide diversity was equal within the Northeastern and Southeastern groups ( Table 2). This indicates that while Northeastern individuals were less likely to have unique haplotypes, the nucleotide differences between individuals, on average, were not smaller than those among Southeastern individuals.

Nuclear genetic variation
18 N. ocellatus individuals from 11 localities were initially sequenced for the nuclear genes β-globin and RPS8. Polymorphism at the nucleotide level was much lower for nDNA than mtDNA, representing four characters out of 613 (0.653 %) for RPS8, and 13 out of 656 (1.98 %) for β-globin. Subsequently, β-globin was successfully sequenced for 166 individuals in total (some individuals could not be resolved owing to the presence of heterozygosity for multiple length variants). A total of 55 unique alleles were identified, with 39 (5.8 %) of 670 characters variable, and 30 (4.5 %) parsimony informative.
The TCS network for β-globin revealed limited phylogeographic structure and a lack of consistency with mtDNA relationships (Fig. 3). A common allele was observed in all three mitochondrially-defined regions, at 18 out of 27 localities where nuclear data were obtained. The only suggestion of geographic structuring was for alleles 13-35 (excluding 15, 16, and 31-34), which appear restricted to the western part of the island, and alleles 49-53 which were only observed at Lost Falls, Coles Bay, and Bicheno. During SAMOVA values of Φ CT had plateaued already at two groups, reflecting a western group (Mt Oakleigh, Lake St Clair, Lagoon of Islands, and Strathgordon) and the remainder. There was significant but weak population genetic structuring in the nuclear data when populations were grouped according to the inferred mtDNA groups (Table 1). However, in contrast to mtDNA, β-globin structure among populations, and among populations within groups, was more than an order of magnitude greater than that among groups (Table 1). Also in contrast to mtDNA, there was no suggestion of greater β-globin population genetic structuring in the Northeast than the Southeast. However, heterozygosity was higher in the Southeast than the Northeast, consistent with patterns of mtDNA haplotype diversity (Table 2).
Demographic histories of regional groups The Southeastern group was the only group with significant mtDNA Tajima's D and Fu's Fs, and an inability to reject the null hypothesis of population expansion from the mismatch distribution (Fig. 4). The same result was observed when repeating these analyses while excluding Lake Mackenzie from the Southeast group. The Northwest group was not analysed due to the low sample size for this group. Given differences in sample size between the Northeast and Southeast groups (number of individuals per site often large for Northeast sites), we also randomly subsampled these sites to six individuals, but the results were qualitatively identical (non-significant Tajima's D and Fu's Fs, and significant mismatch distribution; haplotype diversity lower, and nucleotide diversity similar to the Southeast group). Extended Bayesian Skyline Plots conducted on multilocus data also indicated a stronger signal of recent population growth in the Southeast than the Northeast during the last 1.5 Myr, with the Northeast experiencing moderate growth and subsequent decline, and consequently no net change across the plot (Fig. 5). Similarly, the relaxed random walk analysis indicated more frequent and extensive movement in the Southeast than the Northeast (Fig. 6). Discussion Extensive genetic sampling of N. ocellatus across its range revealed strong phylogeographic structure at multiple spatial scales. There was broad-scale structure based on mtDNA, with distinct groups corresponding to Tasmania's Northwest, Northeast and Southeast (the latter extending onto the Central Plateau). Additionally, phylogeographic structuring was evident within Northeast and Southeast groups, although the degree of structuring as well as the demographic history of each region differed. The Northeastern group had the highest levels of internal structure and coalescent analyses did not provide evidence for any substantial net demographic expansion; the populations appear to have persisted in isolation at or near their sampling sites during Pleistocene glaciations without substantial change in population size. This contrasts with the lower spatial structuring and evidence of greater demographic and spatial expansion in the Southeastern group. While populations within the Northwestern group also appeared to have some level of genetic structure, low sampling currently makes it difficult to reliably infer the biogeographic history of populations from this region. The nuclear β-globin locus provided corroborating  evidence for the three groups suggested by mtDNA, but only within a population genetic (AMOVA)-rather than phylogeographic-framework. Estimated dates for the divergence of lineages from different regions and demographic expansion were firmly in the Pleistocene, but predate the LGM. The observed phylogeographic patterns suggest that Pleistocene climate-but not the LGM specifically-has been highly influential in shaping the distribution of N. ocellatus across Tasmania.

Demographic history
In reptiles, physiological processes essential for survival and reproduction are tightly linked to thermal opportunities [38,49]. Therefore, the phylogeographic patterns and the demographic histories of N. ocellatus populations are likely to reflect aspects of Tasmania's climatic history. During the Pleistocene, high altitude areas including Ben Lomond, Mt Wellington, Mt Field, and the Tasmanian Central Plateau experienced repeated glacial and periglacial activity [44]. While temperatures in lowland areas also fell during glacial periods, these regions were comparatively warm, experiencing temperatures that were as warm as, or in some cases warmer, than alpine regions inhabited by N. ocellatus today [50]. The contrast of phylogeographic patterns in N. ocellatus between the highly structured Northeastern group and the comparatively less structured Southeastern group likely reflects historic differences in persistence of N. ocellatus populations at or near the sampling sites during the colder, more biologically limiting, glacial periods. The observed diversity of haplotypes and the individually monophyletic status of most sites within the Northeastern group provide support for the climatic suitability of northeast Tasmania for populations of temperate species during or immediately following Pleistocene glaciations, and not just the LGM. Additional support is provided by phylogeographic studies in two plant species, Eucalyptus regnans [51] and Nothofagus cunninghamii [52], as well as vegetation modelling based on palynological and ecological data, which predicted refugia in this region for rainforest, grassy woodland and eucalypt species during the LGM [50]. In contrast, the non-monophyly of N. ocellatus at Ben Lomond is consistent with recolonisation from multiple, low elevation sources.
The monophyletic population structure found in some populations in the Southeastern group including Coles Bay, Fortescue Bay and the combined Russell River and Margate population, indicates that in some of these lowland regions, populations may have either persisted in isolation for longer periods than other Southeast populations, or reflect founder effects from peripheral range expansion, the latter suggested by relaxed random walk analysis (Fig. 6). Regardless, there was a higher level of movement throughout the region presently occupied by the Southeastern group; the genetic structure of populations within this region is comparatively low, with fewer mutations separating individuals from different populations, and with similar or even identical haplotypes shared amongst sites. These phylogeographic patterns are consistent with a history of genetic mixture and demographic instability of populations within the Southeast. Populations likely retreated from temporarily unsuitable locations during glacial periods, becoming conjoined and admixed with other retreating populations or possibly Fig. 6 Spatial projection of the relaxed random walk analysis of Niveoscincus ocellatus, based on the maximum clade credibility tree (black lines). Coloured areas reflect the 80 % Highest Posterior Density of the distribution of ancestral branches, with dark blue representing the oldest distribution, and red the youngest with persisting populations. Range expansion would have occurred when the climate became less severe and more habitats became available, consistent with the observed pattern of demographic and spatial expansion within the Southeastern group. Expansion probably occurred simultaneously from multiple refugial populations, as indicated by the high levels of genetic diversity and non-monophyly of sites such as Lake Mackenzie.
The patterns of glacial response inferred for other Tasmanian taxa bear some similarities with classical patterns in the Northern Hemisphere, where post-glacial recolonisation has often been rapid, and from a single or few source populations at lower latitudes and altitudes, leading to decreases in local genetic diversity [3]. For example, two large lowland refugia in the Storm Bay and Great Oyster Bay (Fig. 1) areas have been inferred for Eucalyptus globulus [32]. However, for N. ocellatus relaxed random walk analysis suggested persistence or early establishment in Northern Central Tasmania, and Northeast Tasmania immediately thereafter. This bears more similarity to temperate European species that persisted during glacials in 'cryptic high latitude refugia' , including amphibians and reptiles [53][54][55][56]. Furthermore, often what has represented single large refugia in some species (and has contributed to the idea of the classical European patterns) represents multiple glacial refugia for reptiles (e.g., [8,9,57,58]), and may also apply here for N. ocellatus.
While molecular clocks are limited in the temporal accuracy they provide [59,60], estimated dates for lineage divergence and demographic expansion of the Southeastern clade fall well within the Pleistocene period, but predate the LGM (17-20 ka) [25]. The comparatively shallow divergence among the mtDNA clades (cf. [17]), as well as the low sequence variability in the nDNA, and lack of congruent phylogeographic structure across mitochondrial and nuclear loci, is also consistent with recent (but pre-LGM) diversification within N. ocellatus. The larger effective size of the nDNA markers confers slower lineage sorting relative to mtDNA, which is the favoured explanation for the phylogeographic discordance of these markers (see also [8,61]).

Phylogeographic breaks
The location of phylogeographic breaks between the Northwestern, Northeastern and Southeastern clades in N. ocellatus do not correspond with any large-scale barriers in the current landscape or match the distribution of any previously reported biogeographic regions [62][63][64]. However, similar genetic patterns were observed from an RFLP survey of mtDNA variation in N. metallicus [65], where significant haplotype frequency differences were evident between the Southeast, Northeast, and Northwest, but sampling distribution, sample sizes, and marker variability in that study were too low for robust comparisons with our patterns. In contrast, there is no evidence of phylogeographic breaks within other Tasmanian reptiles [66][67][68], and in one instance levels of divergence to mainland conspecifics suggest colonisation after the LGM, but prior to the flooding of Bass Strait [69]; it is important to note that these taxa have more northern distributions than N. ocellatus, and are also non-saxicoline, and therefore it might not necessarily be expected that they share phylogeographic histories with N. ocellatus.
Niveoscincus ocellatus is a habitat specialist with low dispersal capacity [45,70,71], and these traits are known to promote population fragmentation, as relatively minor breaks in habitat distribution over short distances can often represent effective barriers to gene flow [72,73]. If connectivity between populations is low enough it is possible that even geographically close populations may become highly divergent through time. This may account for the divergence between the Coles Bay and Bicheno populations, where although sites are separated by just 30 km, rock coverage is limited and sandy beaches and dense heath dominate this region.
While poor connectivity between populations can explain why dispersal might be rare, it cannot be seen to prevent dispersal completely in N. ocellatus given that it currently occupies a wide range of locations not presently connected by rocky habitats. Waters et al. [5] highlight a number of systems where phylogeographic breaks have been linked to density-dependent effects, such as intra-specific competitive exclusion, without current geographic barriers to gene flow. These examples demonstrate that with sufficient levels of intra-specific competition, invaders (at low densities) will be blocked by established (high density) resident populations. This exclusion can operate actively, whereby residents are more likely to win territorial disputes (see Olsson and Shine [74] for an example of this in a congeneric species), restricting invader access to resources and mates. Additionally, this competition can occur at the genetic level, where even if invaders occasionally reproduce successfully, without a selective advantage their genes can be effectively excluded by random genetic drift [5]. There is also evidence for adaptation to local environmental conditions in this species [41,75,76], which might further impede gene flow.
Density-dependent effects can maintain barriers to gene flow at a variety of spatial scales [5]. This means that the same processes that maintain the phylogenetic breaks between the major clades of N. ocellatus can also explain the maintenance of structure observed within clades such as those in the Northeast, where unique haplotypes are almost always associated with just one population. In addition, because competitive exclusion is density-dependent, it can only explain the maintenance of phylogeographic breaks in demographically stable areas. However, if climatic conditions change and there is an increase in the density of invaders, or an absence of residents in newly available suitable habitat, genes are more likely to move throughout a landscape, and result in a lack of phylogeographic structuring.

Conclusions
This study indicates that populations of N. ocellatus became established and persisted at or near most contemporary sampling locations in Northeast Tasmania, and some low elevation sites in Southeast Tasmania, throughout the LGM and earlier Pleistocene glaciations. Other regions of central and Southeast Tasmania were uninhabited during the Pleistocene glaciations and were not recolonised from a single refugial population. Our findings contribute to the emerging picture of complexity in phylogeographic histories and of responses to Pleistocene glaciations across species and regions globally. Future studies should compare phylogeographic patterns in similar taxa (e.g., N. metallicus) to verify the inferences of this study, and also examine the fine-scale spatial distribution of genetic variation between areas identified as Northeast and Southeast clades, to better understand the origin and maintenance of these groups (e.g., [77]).

Sampling
Tissue samples were collected from 204 specimens representing 29 localities ( Fig. 1; Additional file 1: Table S1). Sampling was spread across the species' known latitudinal, longitudinal and altitudinal range (Fig. 1). Field samples included some museum collections (Additional file 1: Table S2) to ensure comprehensive sampling across the species range. Niveoscincus greeni and N. pretiosus were employed as outgroups.

DNA extraction, amplification and sequencing
Total genomic DNA was extracted from tail or liver tissue samples following the Qiagen DNeasy Blood and Tissue Extraction Kit (Qiagen, Hilden, Germany). A portion of two mitochondrial genes, ND2 (c. 547 bp) and ND4 (c. 873 bp), and two nuclear genes, β-globin (exons 2-3, c. 656 bp) and RPS8 (40S Ribosomal protein S8 intron 3, c. 613 bp), were amplified and sequenced, as they exhibited useful levels of intraspecific variability in previous population studies of squamates (e.g., [17,78]). The primers used were L4473 and ND2r102 for ND2 [79,80], ND4l and tRNA-Leu for ND4 [81], LC17 and LC18 for RPS8 [17], and Bglo1CR and Bglo2CR for β-globin [82]. A subset of individuals were initially sequenced in both directions for all fragments, but then only L4437 and tRNA_Leu were employed for subsequent mtDNA sequencing. For the nuclear fragments, only β-globin exhibited sufficient variation for further analysis. Most individuals were sequenced using Bglo1CR, but in instances of heterozygotes for length variation, Bglo2CR was also used for sequencing. PCR reactions were 21 μL and com-prised~20 ng of DNA template, 1 x reaction buffer, 1.9 mM MgCl 2 , 0.14 mM dNTPs, 0.24 mM of each primer, and 0.16 units of Taq DNA polymerase (Bioline MangoTaq). For ND2, ND4 and RPS8 amplification was carried out with an initial denaturation at 95°C for 3 min, followed by 35 cycles of 95°C denaturation for 25 s, 55°C annealing for 30 s, and 72°C extension for 90 s, with a final extension at 72°C for 5 min after cycling. Amplification for β-globin followed the same protocol but with annealing at 60°C.
DNA sequence chromatograms were edited and aligned in Geneious version 5.6.4 (Biomatters). Sequence data for ND2 and ND4 were concatenated, providing 1420 bp per individual. TCS v1.21 [83] was used to assign individuals to haplotypes, with attention to ensure that sequences that shared missing data at a given nucleotide position were not assigned to different haplotypes. A single representative of each mitochondrial haplotype was then chosen for phylogeny reconstruction. Given heterozygous nucleotide sites in β-globin, data analyses were conducted on inferred haplotype states rather than using a single sequence with ambiguous nucleotide codes for heterozygous sites. PHASE 2.1.1 [84] was used to phase alleles from heterozygous individuals, with SEQPHASE [85] employed for input and output file processing. DNA sequences analysed are available from GenBank (accession numbers: KJ858058-KJ858494, KP277213-KP277500) and also in Dryad [86].

Phylogenetic analysis
Maximum parsimony and Bayesian tree-building methods were used for mtDNA data. Maximum parsimony trees were constructed in PAUP* v4.0 [87] using equal weighting of character-state changes, random stepwise sequence addition, and tree-bisection reconnection (TBR) branch swapping. The maximum number of recovered trees was set to 5000, and node support was estimated using 500 bootstrap replicates of the dataset. Bayesian analyses were performed with MrBayes 3.2.2 [88]. Data were partitioned by gene region and analysed under substitution models suggested by JModelTest2 [89,90] using the Bayesian Information Criterion. MrBayes analyses were performed using the default prior probability distributions for model parameters and duplicate MCMC runs were conducted with four chains of 5,500,000 generations, with a tree sampling frequency of 500 generations. Three of the chains were heated according to 'Temp = 0.1'. Stationarity and adequate mixing of chains was assessed using Tracer 1.5 [91] to establish the attainment of asymptotes for LnL and substitution model parameters, and to ensure effective sample sizes greater than 200. Convergence was determined when the average standard deviation of split frequencies between runs was lower than 0.02. The first 25 % of the samples were discarded as 'burn in' given attainment of stationarity and convergence. Phylogenetic trees were not reconstructed for β-globin given shallow levels of genetic variation. TCS networks were reconstructed for both β-globin and mtDNA using PopART (http://popart.otago.ac.nz/index.shtml).

Population structuring and demographic histories
Populations were assigned to a regional group based on patterns in the mitochondrial DNA phylogeny and Spatial Analysis of Molecular Variance (SAMOVA; [92]). SAMOVA employed pairwise sequence difference and a range of values for the number of genetic groups, implemented using SAMOVA2 (http://cmpg.unibe.ch/software/samova2/). Population genetic statistics were then calculated for these groups using Arlequin 3.1 [93]. Haplotype diversity and nucleotide diversity were calculated for each group. Population structuring was quantified using Analysis of Molecular Variance (AMOVA) based on haplotype identity. AMOVA calculated hierarchical fixation indices among populations (Φ ST ), among populations within groups (Φ SC ), and among groups (Φ CT ), and P-values were estimated from 1000 permutations of haplotypes. Tajima's D [94] and Fu's Fs [95] were used to test for a signature of demographic expansion for each regional group, assuming selective neutrality of DNA variation. Tests were based on 1000 simulated samples using a coalescent algorithm. Additionally, mismatch distribution analyses were performed to test for exponential demographic expansion [96].
Changes in effective population size were reconstructed using the multilocus Extended Bayesian Skyline Plot (EBSP) [97]. Substitution models were derived from jModelTest2, and a linear model of population size change was implemented. A strict molecular clock was assumed given the shallow levels of divergence, with a mitochondrial mean rate of 1.52 % divergence per million years, with a standard deviation of 0.5 % divergence, based on mtDNA calibrations from other squamates [69]. Analyses were performed using BEAST 1.8.2 [98] with 250 million generations, and sampling every 25,000 generations. Operators were adjusted in accordance with recommendations in the EBSP tutorial at http://beast. bio.ed.ac.uk/tutorials to improve mixing. Replicate runs were tested for convergence visually using Tracer v1.5 [99] and to ensure high effective sample sizes (ESS).
A relaxed random-walk analysis was also employed to reconstruct the spatiotemporal pattern of movements. This approach has similarities to relaxed molecular clock analyses in that it relaxes the constraint of constant movement rates through time [100]. A log-normal prior of movement rates was employed. Analyses were conducted on mtDNA data using BEAST with a coalescent tree prior and strict clock. Locations of individuals were jittered by 0.1°latitude and longitude where they were collected at the same site. Replicate analyses of 250 million generations were run to ensure convergence, and a maximum clade credibility tree was constructed using TreeAnnotater v1.8.2 [98]. Visualisation of the spatiotemporal tree was achieved using SPREAD v1.0.6 [101].

Availability of supporting data
The data set(s) supporting the results of this article are available in the Dryad repository, http://datadryad.org/ review?doi=doi:10.5061/dryad.ff32k.

Additional files
Below is the link to the electronic supplementary material.
Additional file 1: Table S1. Sampling details of specimens. Table S2. Museum specimens. Table S3. Sampling locations of haplotypes illustrated in Fig. 3. Table S4 Sampling locations of haplotypes indicated in Figure S3. Figure S1. Maximum parsimony tree for mtDNA variation. Figure S2. Maximum clade credibility tree from BEAST analysis. Figure S3. Minimum spanning network for mtDNA variation.