- Research article
- Open Access
Persistence and dispersal in a Southern Hemisphere glaciated landscape: the phylogeography of the spotted snow skink (Niveoscincus ocellatus) in Tasmania
BMC Evolutionary Biology volume 15, Article number: 121 (2015)
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.
The glacial-interglacial oscillations of the Pleistocene have strongly influenced the distribution and evolution of many species . Initial interest in the effects of these oscillations largely focused on temperate species within Europe . 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 ), phylogeographic methods can locate refugia by assessing patterns of genetic diversity across a species’ geographical range .
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–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–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–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–24]), while Tasmania—the focus of Australia’s most extensive Pleistocene glacial activity -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 km2) during the early (~1.8 Myr) and middle Pleistocene (>130 kyr) [26–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  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 . 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) , 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 . In this context, lizards have recently been highlighted as good models for use in phylogeographic studies . As ectotherms, lizards are sensitive to changes in climate, which may be manifested in alterations of species distributions . Lizards also typically have low mobility, and phylogeographic patterns from historic events will be retained for longer periods .
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 above ~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 . 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 , which suggests that gene flow will be generally low  (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. 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 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).
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 pre-date 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.
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 . 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 .
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  and Nothofagus cunninghamii , 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 . 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 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 . For example, two large lowland refugia in the Storm Bay and Great Oyster Bay (Fig. 1) areas have been inferred for Eucalyptus globulus . 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–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) . The comparatively shallow divergence among the mtDNA clades (cf. ), 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]).
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–64]. However, similar genetic patterns were observed from an RFLP survey of mtDNA variation in N. metallicus , 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–68], and in one instance levels of divergence to mainland conspecifics suggest colonisation after the LGM, but prior to the flooding of Bass Strait ; 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.  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  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 . 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 . 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.
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., ).
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 , LC17 and LC18 for RPS8 , and Bglo1CR and Bglo2CR for β-globin . 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 comprised ~20 ng of DNA template, 1 x reaction buffer, 1.9 mM MgCl2, 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  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  was used to phase alleles from heterozygous individuals, with SEQPHASE  employed for input and output file processing. DNA sequences analysed are available from GenBank (accession numbers: KJ858058–KJ858494, KP277213–KP277500) and also in Dryad .
Maximum parsimony and Bayesian tree-building methods were used for mtDNA data. Maximum parsimony trees were constructed in PAUP* v4.0  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 . 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  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; ). 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 . 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  and Fu’s Fs  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 .
Changes in effective population size were reconstructed using the multilocus Extended Bayesian Skyline Plot (EBSP) . 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 . Analyses were performed using BEAST 1.8.2  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  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 . 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 . Visualisation of the spatiotemporal tree was achieved using SPREAD v1.0.6 .
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.
Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405(6789):907–13.
Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23(10):564–71.
Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond B Biol Sci. 2004;359(1442):183–95.
Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F. Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998;7(4):453–64.
Waters JM, Fraser CI, Hewitt GM. Founder takes all: density-dependent processes structure biodiversity. Trends Ecol Evol. 2013;28(2):78–85.
Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Proc R Soc B Biol Sci. 2010;277(1682):661–71.
Schmitt T, Varga Z. Extra-Mediterranean refugia: the rule and not the exception. Fron t Zool. 2012;9:22.
Salvi D, Harris DJ, Kaliontzopoulou A, Carretero MA, Pinho C. Persistence across Pleistocene ice ages in Mediterranean and extra-Mediterranean refugia: phylogeographic insights from the common wall lizard. BMC Evol Biol. 2013;13(1):147.
Gómez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In: Phylogeography of southern European refugia. Dordrecht, The Netherlands: Springer; 2007: 155–188.
Beheregaray LB. Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol Ecol. 2008;17(17):3754–74.
Dubey S, Shine R. Geographic variation in the age of temperate-zone reptile and amphibian species: Southern Hemisphere species are older. Biol Lett. 2011;7(1):96–7.
Dubey S, Shine R. Are reptile and amphibian species younger in the Northern Hemisphere than in the Southern Hemisphere? J Evol Biol. 2012;25(1):220–6.
Tingley R, Dubey S. Disparity in the timing of vertebrate diversification events between the northern and southern hemispheres. BMC Evol Biol. 2012;12.
Turchetto-Zolet AC, Pinheiro F, Salgueiro F, Palma-Silva C. Phylogeographical patterns shed light on evolutionary process in South America. Mol Ecol. 2013;22(5):1193–213.
Fraser CI, Nikula R, Ruzzante DE, Waters JM. Poleward bound: biological impacts of Southern Hemisphere glaciation. Trends Ecol Evol. 2012;27(8):462–71.
Costa LP. The historical bridge between the Amazon and the Atlantic Forest of Brazil: a study of molecular phylogeography with small mammals. J Biogeogr. 2003;30(1):71–86.
Bell RC, Parra JL, Tonione M, Hoskin CJ, Mackenzie JB, Williams SE, et al. Patterns of persistence and isolation indicate resilience to climate change in montane rainforest lizards. Mol Ecol. 2010;19(12):2531–44.
Hoskin CJ, Tonione M, Higgie M, MacKenzie JB, Williams SE, VanDerWal J, et al. Persistence in peripheral refugia promotes phenotypic divergence and speciation in a rainforest frog. Am Nat. 2011;178(5):561–78.
Barlow A, Baker K, Hendry CR, Peppin L, Phelps T, Tolley KA, et al. Phylogeography of the widespread African puff adder (Bitis arietans) reveals multiple Pleistocene refugia in southern Africa. Mol Ecol. 2013;22(4):1134–57.
Ritchie PA, Millar CD, Gibb GC, Baroni C, Lambert DM. Ancient DNA enables timing of the Pleistocene origin and Holocene expansion of two Adélie penguin lineages in Antarctica. Mol Biol Evol. 2004;21(2):240–8.
Fraser CI, Nikula R, Spencer HG, Waters JM. Kelp genes reveal effects of subantarctic sea ice during the Last Glacial Maximum. Proc Natl Acad Sci. 2009;106(9):3249–53.
Wallis GP, Trewick SA. New Zealand phylogeography: evolution on a small continent. Mol Ecol. 2009;18(17):3548–80.
McGaughran A, Torricelli G, Carapelli A, Frati F, Stevens MI, Convey P, et al. Contrasting phylogeographical patterns for springtails reflect different evolutionary histories between the Antarctic Peninsula and continental Antarctica. J Biogeogr. 2010;37(1):103–19.
Sersic AN, Cosacov A, Cocucci AA, Johnson LA, Pozner R, Avila LJ, et al. Emerging phylogeographical patterns of plants and terrestrial vertebrates from Patagonia. Biol J Linn Soc. 2011;103(2):475–94.
Barrows TT, Stone JO, Fifield LK, Cresswell RG. The timing of the last glacial maximum in Australia. Quat Sci Rev. 2002;21(1):159–73.
Colhoun EA. Glaciations of the west coast range. Tasmania. Quatern Res. 1985;24(1):39–59.
Kiernan K. The extent of late Cenozoic glaciation in the Central Highlands of Tasmania, Australia. Arct Alp Res. 1990;341–354.
Colhoun EA, Hannan D, Kiernan K. Late Wisconsin glaciation of Tasmania. In: Papers and proceedings of the royal society of Tasmania. 1996. p. 33–45.
Augustinus PC, Macphail MK. Early Pleistocene stratigraphy and timing of the Bulgobac Glaciation, western Tasmania, Australia. Palaeogeogr Palaeoclimatol Palaeoecol. 1997;128(1):253–67.
Kiernan K, Lauritzen SE, Duhig N. Glaciation and cave sediment aggradation around the margins of the Mt Field Plateau, Tasmania. Aust J Earth Sci. 2001;48(2):251–63.
Mackintosh A, Barrows T, Colhoun E, Fifield L. Exposure dating and glacial reconstruction at Mt. Field, Tasmania, Australia, identifies MIS 3 and MIS 2 glacial advances and climatic variability. J Quat Sci. 2006;21(4):363–76.
McKinnon GE, Jordan GJ, Vaillancourt RE, Steane DA, Potts BM. Glacial refugia and reticulate evolution: the case of the Tasmanian eucalypts. Philos Trans R Soc Lond B Biol Sci. 2004;359(1442):275–84.
Byrne M. Evidence for multiple refugia at different time scales during Pleistocene climatic oscillations in southern Australia inferred from phylogeography. Quat Sci Rev. 2008;27(27):2576–85.
Rull V. Microrefugia. J Biogeogr. 2009;36(3):481–4.
Lessa EP, D’Elia G, Pardinas UFJ. Genetic footprints of late Quaternary climate change in the diversity of Patagonian-Fueguian rodents. Mol Ecol. 2010;19(15):3031–7.
Zhang Z-Y, Cashins S, Philips A, Burridge CP. Significant population genetic structuring but a lack of phylogeographic structuring in the endemic Tasmanian tree frog (Litoria burrowsae). Aust J Zool. 2014;62(3):238–45.
Camargo A, Sinervo B, Sites JW. Lizards as model organisms for linking phylogeographic and speciation studies. Mol Ecol. 2010;19(16):3250–70.
Angilletta MJ, Niewiarowski PH, Navas CA. The evolution of thermal physiology in ectotherms. J Therm Biol. 2002;27(4):249–68.
Araújo MB, Nogués-Bravo D, Diniz-Filho JAF, Haywood AM, Valdes PJ, Rahbek C. Quaternary climate changes explain diversity among reptiles and amphibians. Ecography. 2008;31(1):8–15.
Rawlinson P: Biogeography and ecology of the reptiles of Tasmania and the Bass Strait area. In: Biogeography and ecology in Tasmania. Springer; 1974: 291–338
Uller T, While GM, Cadby CD, Harts A, O’Connor K, Pen I, et al. Altitudinal divergence in maternal thermoregulatory behaviour may be driven by differences in selection on offspring survival in a viviparous lizard. Evolution. 2011;65(8):2313–24.
Wapstra E, Swain R. Geographic and annual variation in life-history traits in a temperate zone Australian skink. J Herpetol. 2001;194–203.
Olsson M, Shine R, Wapstra E, Ujvari B, Madsen T. Sexual dimorphism in lizard body shape: the roles of sexual selection and fecundity selection. Evolution. 2002;56(7):1538–42.
Colhoun EA. Periglacial landforms and deposits of Tasmania: Periglacial and Permafrost Research in the Southern Hemisphere. S Afr J Sci. 2002;98(1 & 2):55–63.
Melville J. Evolutionary correlations between microhabitat specialisation and locomotor capabilities in the lizard genus Niveoscincus. Aust J Zool. 2007;55(6):351–5.
Gardner MG, Godfrey SS, Fenner AL, Donnellan SC, Bull CM. Fine-scale spatial structuring as an inbreeding avoidance mechanism in the social skink Egernia stokesii. Aust J Zool. 2012;60(4):272–7.
Berry O, Tocher MD, Sarre SD. Can assignment tests measure dispersal? Mol Ecol. 2004;13(3):551–61.
Dennison S, Smith SM, Stow AJ. Long-distance geneflow and habitat specificity of the rock-dwelling coppertail skink, Ctenotus taeniolatus. Austral Ecol. 2012;37(2):258–67.
Sinervo B, Mendez-De-La-Cruz F, Miles DB, Heulin B, Bastiaans E, Villagrán-Santa Cruz M, et al. Erosion of lizard diversity by climate change and altered thermal niches. Science. 2010;328(5980):894–9.
Kirkpatrick JB, Fowler M. Locating likely glacial forest refugia in Tasmania using palynological and ecological information to test alternative climatic models. Biol Conserv. 1998;85(1):171–82.
Nevill PG, Bossinger G, Ades PK. Phylogeography of the world’s tallest angiosperm, Eucalyptus regnans: evidence for multiple isolated Quaternary refugia. J Biogeogr. 2010;37(1):179–92.
Worth JR, Jordan GJ, McKinnon GE, Vaillancourt RE. The major Australian cool temperate rainforest tree Nothofagus cunninghamii withstood Pleistocene glacial aridity within multiple regions: evidence from the chloroplast. New Phytol. 2009;182(2):519–32.
Ursenbacher S, Carlsson M, Helfer V, Tegelström H, Fumagalli L. Phylogeography and Pleistocene refugia of the adder (Vipera berus) as inferred from mitochondrial DNA sequence data. Mol Ecol. 2006;15(11):3425–37.
Ursenbacher S, Conelli A, Golay P, Monney J-C, Zuffi M, Thiery G, et al. Phylogeography of the asp viper (Vipera aspis) inferred from mitochondrial DNA sequence data: Evidence for multiple Mediterranean refugial areas. Mol Phylogenet Evol. 2006;38(2):546–52.
Hofman S, Spolsky C, Uzzell T, COGĂLNICEANU D, BABIK W, Szymura JM. Phylogeography of the fire‐bellied toads Bombina: independent Pleistocene histories inferred from mitochondrial genomes. Mol Ecol. 2007;16(11):2301–16.
Teacher AGF, Garner TWJ, Nichols RA. European phylogeography of the common frog (Rana temporaria): routes of postglacial colonization into the British Isles, and evidence for an Irish glacial refugium. Heredity. 2009;102(5):490–6.
Podnar M, Mayer W, Tvrtković N. Phylogeography of the Italian wall lizard, Podarcis sicula, as revealed by mitochondrial DNA sequences. Mol Ecol. 2005;14(2):575–88.
Canestrelli D, Cimmaruta R, Nascetti G. Population genetic structure and diversity of the Apennine endemic stream frog, Rana italica–insights on the Pleistocene evolutionary history of the Italian peninsular biota. Mol Ecol. 2008;17(17):3856–72.
Bromham L, Penny D. The modern molecular clock. Nat Rev Genet. 2003;4(3):216–24.
Ho SYW, Larson G. Molecular clocks: when times are a-changin’. Trends Genet. 2006;22(2):79–83.
Miraldo A, Hewitt GM, Paulo OS, Emerson BC. Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones. BMC Evol Biol. 2011;11(1):170.
Shiel R, Koste W, Tan L: Tasmania revisited: rotifer communities and habitat heterogeneity. In: Rotifer Symposium V: 1989. Springer: 239–245.
Mesibov R. Faunal breaks in Tasmania and their significance for invertebrate conservation. Mem Queensland Mus . 1994;36(1):133–6.
Vélez S, Mesibov R, Giribet G. Biogeography in a continental island: population structure of the relict endemic centipede Craterostigmus tasmanianus (Chilopoda, Craterostigmomorpha) in Tasmania using 16S rRNA and COI. J Hered. 2012;103(1):80–91.
McCoull CJ: Geographic variation and adaptation in the Tasmanian metallic skink (Niveoscincus metallicus) PhD thesis, University of Tasmania; 2001.
Keogh JS, Scott IA, Hayes C. Rapid and repeated origin of insular gigantism and dwarfism in Australian tiger snakes. Evolution. 2005;59(1):226–33.
Dubey S, Shine R. Evolutionary diversification of the lizard genus Bassiana (Scincidae) across Southern Australia. PLoS ONE. 2010;5(9):e12982.
Chapple DG, Hoskin CJ, Chapple SNJ, Thompson MB. Phylogeographic divergence in the widespread delicate skink (Lampropholis delicata) corresponds to dry habitat barriers in eastern Australia. BMC Evol Biol. 2011;11.
Chapple DG, Chapple SN, Thompson MB. Biogeographic barriers in south-eastern Australia drive phylogeographic divergence in the garden skink, Lampropholis guichenoti. J Biogeogr. 2011;38(9):1761–75.
Atkins N, Swain R, Wapstra E, Jones SM. Late stage deferral of parturition in the viviparous lizard Niveoscincus ocellatus (Gray 1845): implications for offspring quality and survival. Biol J Linn Soc. 2007;90(4):735–46.
Wapstra E, Uller T, While G, Olsson M, Shine R. Giving offspring a head start in life: field and experimental evidence for selection on maternal basking behaviour in lizards. J Evol Biol. 2010;23(3):651–7.
Henle K, Lindenmayer DB, Margules CR, Saunders DA, Wissel C. Species survival in fragmented landscapes: where are we now? Biodivers Conserv. 2004;13(1):1–8.
Ayre DJ, Minchinton TE, Perrin C. Does life history predict past and current connectivity for rocky intertidal invertebrates across a marine biogeographic barrier? Mol Ecol. 2009;18(9):1887–903.
Olsson M, Shine R. Ownership influences the outcome of male-male contests in the scincid lizard, Niveoscincus microlepidotus. Behav Ecol. 2000;11(6):587–90.
Pen I, Uller T, Feldmeyer B, Harts A, While GM, Wapstra E. Climate-driven population divergence in sex-determining systems. Nature. 2010;468(7322):436–8.
Cadby CD, Jones SM, Wapstra E. Geographical differences in maternal basking behaviour and offspring growth rate in a climatically widespread viviparous reptile. J Exp Biol. 2014;217(7):1175–9.
Phillips BL, Baird SJ, Moritz C. When vicars meet: a narrow contact zone between morphologically cryptic phylogeographic lineages of the rainforest skink, Carlia rubrigularis. Evol. 2004;58(7):1536–48.
Greaves SN, Chapple DG, Gleeson DM, Daugherty CH, Ritchie PA. Phylogeography of the spotted skink (Oligosoma lineoocellatum) and green skink (O. chloronoton) species complex (Lacertilia: Scincidae) in New Zealand reveals pre-Pleistocene divergence. Mol Phylogenet Evol. 2007;45(2):729–39.
Macey JR, Larson A, Ananjeva NB, Fang Z, Papenfuss TJ. Two novel gene orders and the role of light-strand replication in rearrangement of the vertebrate mitochondrial genome. Mol Biol Evol. 1997;14(1):91–104.
Sadlier RA, Smith SA, Bauer AM, Whitaker AH. A new genus and species of live-bearing scincid lizard (Reptilia: Scincidae) from New Caledonia. J Herpetol. 2004;38(3):320–30.
Forstner MRJ, Davis SK, Arévalo E. Support for the hypothesis of anguimorph ancestry for the suborder Serpentes from phylogenetic analysis of mitochondrial DNA sequences. Mol Phylogenet Evol. 1995;4(1):93–102.
Dolman G, Phillips B. Single copy nuclear DNA markers characterized for comparative phylogeography in Australian wet tropics rainforest skinks. Mol Ecol Notes. 2004;4(2):185–7.
Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–9.
Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68(4):978–89.
Flot JF. SeqPHASE: a web tool for interconverting PHASE input/output files and FASTA sequence alignments. Mol Ecol Res. 2010;10(1):162–6.
Cliff HC, Wapstra E, Burridge CP: Data from: Persistence and dispersal in a Southern Hemisphere glaciated landscape: the phylogeography of the spotted snow skink (Niveoscincus ocellatus) in Tasmania. Dryad. 2015. http://datadryad.org/review?doi=doi:10.5061/dryad.ff32k.
Swofford DL. PAUP*. Phylogenetic analysis using parsimony (* and other methods). In: Version 4. Sunderland, Massachusetts: Sinauer Associates; 2003.
Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19(12):1572–4.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.
Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52(5):696–704.
Rambaut A, Suchard MA, Xie D, Drummond AJ: Tracer v1.5. In.: Available from http://beast.bio.ed.ac.uk/Tracer; 2013.
Dupanloup I, Schneider S, Excoffier L. A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002;11(12):2571–81.
Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005;1:47.
Tajima F. The effect of change in population size on DNA polymorphism. Genet. 1989;123(3):597–601.
Fu Y-X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genet. 1997;147(2):915–25.
Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9(3):552–69.
Heled J, Drummond AJ. Bayesian inference of population size history from multiple loci. BMC Evol Biol. 2008;8.
Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian Phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29(8):1969–73.
Rambaut A, Drummond A: Tracer v1. 4. In.: Available from http://beast.bio.ed.ac.uk/Tracer; 2007.
Lemey P, Rambaut A, Welch JJ, Suchard MA. Phylogeography takes a relaxed random walk in continuous space and time. Mol Biol Evol. 2010;27(8):1877–85.
Bielejec F, Rambaut A, Suchard MA, Lemey P. SPREAD: spatial phylogenetic reconstruction of evolutionary dynamics. Bioinformatics. 2011;27(20):2910–2.
We are grateful to staff at the CSIRO Australian National Wildlife Collection, the South Australian Museum and the Queen Victoria Museum and Art Gallery as well as C. Jennings for assistance in accessing some of the samples used in this study. We would like to thank the volunteers who helped collect samples in the field: M. Caldwell, A. Demir, T. Feldmanis, J. Gruber, K. Munch, P. Richardson and E. K. Yuni, with particular thanks to R. De Paoli. We would also like to thank T. Wenner and A. Pracejus for valuable assistance in the laboratory, J. Kirkpatrick, G. Jordon, J. Worth for insights on biotic responses to glaciations in Tasmania, and C. Moritz and D. Rosauer for technical support in the production of the map included in this publication. Samples were collected with approval from the University of Tasmania Animal Ethics Committee (permit A12286). This research was funded by the Australian Research Council (FT110100597).
The authors declare that they have no competing interests.
CPB and EW conceived and designed the study. HBC and EW conducted fieldwork. HC and CB conducted laboratory work and data analyses. All authors contributed to the writing of the manuscript and approved the final version. All authors read and approved the final manuscript.
H.B.C. conducted this research as part of her BSc(Honours) degree, and is broadly interested in biogeography and conservation of terrestrial organisms. E.W. is a behavioural ecologist that concentrates his research on reptiles, and has ongoing interest in species responses to climate change. C.P.B. is a molecular ecologist with interests in historical biogeography, gene flow, and phylogenetics.
Below is the link to the electronic supplementary material.
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.
About this article
Cite this article
Cliff, H., Wapstra, E. & Burridge, C. Persistence and dispersal in a Southern Hemisphere glaciated landscape: the phylogeography of the spotted snow skink (Niveoscincus ocellatus) in Tasmania. BMC Evol Biol 15, 121 (2015). https://doi.org/10.1186/s12862-015-0397-y
- Glacial refugia