Phylogeography of the Endangered Otago Skink, Oligosoma otagense: Population Structure, Hybridisation and Genetic Diversity in Captive Populations

Climatic cooling and substantial tectonic activity since the late Miocene have had a pronounced influence on the evolutionary history of the fauna of New Zealand's South Island. However, many species have recently experienced dramatic range reductions due to habitat fragmentation and the introduction of mammalian predators and competitors. These anthropogenic impacts have been particularly severe in the tussock grasslands of the Otago region. The Otago skink (Oligosoma otagense), endemic to the region, is one of the most critically endangered vertebrates in New Zealand. We use mitochondrial DNA sequence data to investigate the evolutionary history of the Otago skink, examine its population genetic structure, and assess the level of genetic diversity in the individuals in the captive breeding program. Our data indicate that the Otago skink diverged from its closest relatives in the Miocene, consistent with the commencement of tectonic uplift of the Southern Alps. However, there is evidence for past introgression with the scree skink (O. waimatense) in the northern Otago-southern Canterbury region. The remnant populations in eastern Otago and western Otago are estimated to have diverged in the mid-Pliocene, with no haplotypes shared between these two regions. This divergence accounts for 95% of the genetic diversity in the species. Within both regions there is strong genetic structure among populations, although shared haplotypes are generally evident between adjacent localities. Although substantial genetic diversity is present in the captive population, all individuals originate from the eastern region and the majority had haplotypes that were not evident in the intensively managed populations at Macraes Flat. Our data indicate that eastern and western populations should continue to be regarded as separate management units. Knowledge of the genetic diversity of the breeding stock will act to inform the captive management of the Otago skink and contribute to a key recovery action for the species.


Introduction
Geological and climatic processes have acted to dramatically alter the landscape of the South Island of New Zealand [1][2][3]. Throughout the Miocene the South Island was an eroded peneplain dominated by rainforest vegetation [4][5][6]. However, New Zealand straddles the boundary of the Indo-Australian and Pacific plates, and tectonic activity along the fault line in the South Island (Alpine Fault) that commenced during the Miocene, and intensified during the Pliocene, resulted in the formation of the Southern Alps (.3000 m) [2,3,7,8]. As New Zealand's climate cooled during the Pliocene-Pleistocene [9], the predominant vegetation transitioned from rainforest to grassland and created an expansive subalpine and alpine zone in the mountainous regions of the South Island [1,6,10]. Extensive glaciers throughout the South Island were evident during glacial maxima [9], and repeated glacial cycles during the Pleistocene resulted in the continual expansion and contraction of the distribution of the resident biota [3,11]. These processes acted to fragment the range of many species, and led to extensive speciation and phylogeographic structure within the South Island (reviewed in [3,11]).
Anthropogenic impacts, particularly since European settlement ,200 years ago, have resulted in the decline or extinction of many native species in the South Island [12,13]. The introduction of a suite of mammalian predators and competitors (e.g. rodents, cats, rabbits, stoats, ferrets, weasels) has been a principal cause of many population declines, with human activities (e.g. agriculture, mining, housing developments) leading to decreased connectivity between populations due to habitat degradation and fragmentation [11,[13][14][15]. Such impacts appear to have been particularly pronounced in the tussock grassland habitats of the central Otago region of the South Island [14,[16][17][18][19].
The Otago skink (Oligosoma otagense), listed as Nationally Critical, is one of the most endangered vertebrates in New Zealand [20]. It is a large-sized (up to 125 mm snout-vent length: SVL) viviparous skink species that is endemic to the montane tussock grasslands of the Otago region [21][22][23]. The Otago skink is diurnal and inhabits schist rock outcrops [22,23]. It exhibits a metapopulation structure, and habitat fragmentation has been reported to limit dispersal between rock outcrops [24,25]. The Otago skink has experienced a substantial population decline since European settlement, and is currently estimated to inhabit only 8% of its former range [20,22,26,27] (Figure 1). The species is currently restricted to two regions at the periphery of its former range: in eastern Otago between Macraes Flat and Sutton, and a few scattered populations between Lake Hawea and Lindis Pass in western Otago [22] (Figure 1). The Otago skink is still declining, with the local extinction of several populations occurring within the last 40 years [22,26]. The largest populations occur in the Macraes Flat region, although most extant populations are extremely small (,60 individuals) [22]. It is estimated that fewer than 5000 Otago skinks remain in the wild [22,24,26,28].
The conservation management of the Otago skink is focused on ,2400 ha of land (managed by the Department of Conservation) near the Macraes Flat township [22]. A mark-recapture study has been conducted in this region since the 1990s to monitor the Otago skink populations [22]. Long-term intensive predator control (i.e. trapping, baiting, shooting) has been conducted at Macraes Flat, and several mammal-proof exclosures have been constructed at the site in an effort to increase the size of the Otago skink population [22]. Regular surveys for the Otago skink are also conducted in the western Otago populations, and a reserve was recently established to encompass the Glenfoyle Station population [22]. There are currently ,100 Otago skinks in captivity, and several generations of offspring have been successfully bred as part of the captive breeding program [22,23]. Although the Otago skink has been managed as a distinct taxon for several decades, it was previously considered to comprise two subspecies (otagense, waimatense; [29]) or forms (form 'otagense', form 'waimatense'; [30]). Oligosoma otagense and O. waimatense (scree skink; up to 110 mm SVL) were only formally recognised as distinct species in 1997, based on body shape and colour pattern [31]. Preliminary allozyme work indicated that hybridisation may occur between the Otago skink and scree skink [32], and the two species have been observed to produce viable hybrids in captivity (D. Keall unpublished data, [33]).
Here we use mitochondrial sequence data to examine the phylogenetic history (ND2, ND4, Cytochrome b) and phylogeography (control region) of the Otago skink. In particular, we investigate the population structure within and between the two regions (eastern Otago, western Otago), determine the level of genetic diversity that is present within the captive breeding stock, and examine whether there is evidence for hybridisation between the Otago skink and scree skink.

Sampling
For our phylogeographic analyses, we obtained tissue samples from 63 Otago skinks, encompassing both the eastern (three sites at Macraes Flat: Falcon, Wildlife and Alistair's Gully; 45u269S, 170u259E) and western populations (Lindis Pass: 44u369S, 169u409E; Glenfoyle station: 44u439S, 169u199E) ( Figure 1). The Wildlife and Falcon sites are within the area that is intensively managed by the Department of Conservation, while Alistair's Gully is on a privately owned farm with grazing and no predator control [23]. These samples were obtained from existing collections (National Frozen Tissue Collection [NFTC], housed at Victoria University of Wellington, New Zealand; Te Papa Tongarewa, National Museum of New Zealand, Wellington) and tissues collected by the New Zealand Department of Conservation during population surveys of the species (1998)(1999)(2000)(2001)(2002)(2003). In addition, tissue samples were collected from 87 Otago skinks held in captivity (by zoos, nature parks, and private breeders) as part of the captive breeding program. As the scree skink is known to hybridise with the Otago skink in captivity [33], we also obtained tissue samples from across the range of the scree skink (7 samples from the NFTC and Te Papa Tongarewa).
A recent molecular phylogeny for the New Zealand skink fauna [34] indicates that the Otago skink is part of a monophyletic lineage that also includes O. waimatense, O. acrinasum, O. infrapunctatum, O. taumakae and O. pikitanga (also see [35][36][37]). Thus, samples of these five species were included, along with 15 representative O. otagense samples, in the broader phylogenetic analyses (Table 1). Oligosoma polychroma and O. oliveri, species from other lineages within the New Zealand skink fauna [34], were used as outgroups in the phylogenetic analyses.

DNA extraction, amplification and sequencing
Total genomic DNA was extracted from liver, muscle, toe or tail-tip samples using a modified phenol-chloroform protocol [38] or a Bio-Rad Aqua Pure Genomic DNA Extraction Kit (Bio-Rad, Hercules CA, USA). For the phylogeographic analyses (i.e. Otago and scree skink samples) we sequenced the mitochondrial control region (,500 bp). For the broader phylogenetic analyses, involving representative Otago skink samples and other species from the same lineage, we sequenced portions of three mitochondrial genes: ND2 (,600 bp), ND4 (,700 bp), and Cytochrome b (,700 bp). These regions were targeted because our previous work across several taxonomic levels in New Zealand skinks has indicated useful levels of variability [34,[39][40][41][42][43][44][45][46][47][48][49][50]. The primers used to amplify and sequence these regions are provided in Table 2. PCR and sequencing were conducted as outlined in Berry & Gleeson [39] (control region) or Greaves et al. [40] (ND2, ND4 and Cytochrome b). PCR products were purified using ExoSAP-IT (USB Corporation, Cleveland, Ohio USA). The purified product was sequenced directly using a BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) and then analysed on an ABI310 or ABI 3730XL capillary sequencer.
Sequence data were edited using CONTIGEXPRESS in VECTOR NTI ADVANCE v9.1.0 (Invitrogen), and aligned using the default parameters of CLUSTAL X v1.83 [51]. We translated all coding region sequences to confirm that none contained premature stop codons. Sequence data were submitted to GenBank [GenBank: JN999929-JN999994] ( Table 1).

Phylogenetic analyses
The concatenated ND2, ND4 and cytochrome b sequence data was used for the phylogenetic analyses. Maximum Likelihood (ML) and Bayesian tree building methods were used. We used MODELTEST 3.7 [52] to identify the most appropriate model of sequence evolution based on the Akaike information criterion (AIC) [53]. MODELTEST, conducted in PAUP* 4.0b10 [54], was also used to estimate base frequencies, substitution rates, the proportion of invariable sites (I) and the among-site substitution rate variation (G) [55]. These values were then used as settings in PhyML 3.0 [56] to generate a ML tree with 500 bootstraps. MRBAYES 3.1.2 [57] was used to complete Bayesian analyses. We used MODELTEST to determine the most appropriate model of sequence evolution for our dataset. We ran the Bayesian analysis for five million generations, sampling every 100 generations (i.e. 50,000 sampled trees). We ran the analysis twice, using four heated chains per run. We discarded the first 25% of samples as burn-in and the last 37,500 trees were used to estimate the Bayesian posterior probabilities. The program TRACER 1.5 [58] was used to check for chain convergence and mixing.
Bootstrap values (500 ML bootstraps) and Bayesian posterior probabilities were used to assess branch support. We considered branches supported by bootstrap values of 70% or greater [59], and/or posterior probability values greater than or equal to 0.95 [60] to be supported by our data.
We estimated the divergence time of the eastern and western Otago skink populations using an evolutionary rate of 1.3-1.63% sequence divergence per million years, based on mitochondrial DNA calibrations from other squamate reptile groups (1.3%, [61]; 1.42-1.63%, [62]; 1.55%, [63]; 1.62%, [64]; 1.63%, [65]). A strict molecular clock (0.0065-0.00815 substitutions per site per million years), implemented in BEAST v1.6.1 [66], was used to estimate the divergence times. The New Zealand skink lineage is estimated to have originated ,20 mya [34,67], and this information was used as the maximum age of the tree root. A GTR+I+G model of evolution was employed with a speciation (Yule) tree prior. The analysis was run twice, with 20 million generations per run (total 40 million generations after the two runs were combined using LOGCOMBINER v1.6.1). The output was viewed in TRACER to check that stationarity had been reached, and ensure that the effective sample size (ESS) exceeded 200 [66].

Phylogeographic analyses: molecular diversity and population structure
The control region dataset was used for the phylogeographic analyses. Estimates of genetic diversity within Otago skink populations (number of haplotypes, h; haplotypic diversity, Hd; number of polymorphic sites, S; nucleotide diversity, p) were calculated in DNASP v4.50 [68]. Tamura-Nei (TrN)-corrected genetic distances within and among populations were calculated in MEGA 4 [69]. Haplotype networks were created using TCS v1.21 [70]. Genetic differentiation among populations was estimated in ARLEQUIN v3.5 [71]. Pairwise W ST values (an analogue of Wright's fixation index F ST ) were calculated to estimate among population differentiation. We conducted hierarchical Analysis of Molecular Variance (AMOVA; [72]) to investigate the partitioning of genetic variation within and among populations and regions. Both tests used TrN genetic distances. Significance levels of all the estimated values were calculated by 10,000 permutations, and adjusted according to the Bonferroni correction procedure [73] for multiple pairwise comparisons as described by Holm [74].
We used Tajima's D [75], Fu's Fs statistic [76] (calculated in ARLEQUIN) and mismatch distributions to test for signatures of  [75]. A smooth bell shape signifies either population expansion or spatial range expansion, whereas a multimodal distribution represents a long history in situ [77][78][79][80]. To distinguish between these two types of distribution, a raggedness index (RI, sum of the squared difference between neighbouring peaks) and the sum of squared deviations (SSD) between the observed and expected mismatch were calculated using the methods of Schneider & Excoffier [81] in ARLEQUIN. The spatial expansion hypothesis (both RI and SSD) was tested using a parametric bootstrap approach (200 replicates).  (Figure 2). In contrast, the O. waimatense individuals from the remainder of the distribution (OWA2, 3, 6) formed a well-supported lineage (100 BS, 1.0 PP), which was genetically divergent (11.5% sequence divergence) from both O. otagense lineages ( Table 3). The mean genetic divergence among recognised species (excluding OWA1, 4-5, 7) was 9.4-13.4% (Table 3). The divergence between the eastern and western populations of O. otagense was 4.9% and estimated to have occurred 3.7 mya (95% highest posterior density [HPD] confidence interval: 2.8-4.6 mya).

Phylogeographic analyses: molecular diversity and population structure
The edited control region alignment comprised 419 characters, with 46 variable sites within the six haplotypes identified within the Otago skink (Table 4, Figure 3a). No haplotypes were shared between the eastern and western populations (Figure 3a). Four haplotypes were found to occur within the Macraes Flat region, with one present in all three populations (Haplotype B), another present at Falcon and Alistair's Gully (Haplotype A), and a further two (Haplotypes C and D) found exclusively at the Wildlife site (Table 4, Figure 3a). Two haplotypes were identified in the western region (Haplotypes E and F), with both occurring in the Lindis Pass population, but only Haplotype E was present in the Glenfoyle Station population (Table 4, Figure 3a). Consistent length variation was observed between the control region haplotypes in the eastern (419 bp) and western (415 bp) regions of the Otago skinks range. Both haplotype and nucleotide diversity was generally low within each population (Table 4).
Substantial genetic differentiation was evident between the two regions (10.6-11.5% sequence divergence, Table 5; pairwise W ST = 0.958), with 36 variable sites separating the eastern and western populations (Figure 3a). Although the level of genetic differentiation among populations within each region was lower (0.2-1.5% sequence divergence, Table 5), significant pairwise W ST values were found among all five populations ( Table 6). The AMOVA confirmed that most genetic variation was partitioned among regions (94.60%), rather than among populations in each region (3.09%), or within populations (2.31%). There was no strong evidence for recent range expansion in any of the five populations (Table 4).
Five haplotypes, each originating from the eastern region, were identified from the 87 captive Otago skinks (Tables 1 and 5;   Table 5). Three unique control region haplotypes were identified in the four scree skinks (Haplotype J: OWA5; Haplotype K: OWA7; Haplotype L: OWA1, OWA4) which appeared, based on our phylogenetic analyses, to derive from lineages that had experienced past introgression with the Otago skink (   and 3c). This introgression appears to have occurred with both the eastern (Haplotype L) and western (Haplotypes J and K) Otago skink populations (Figure 3c). None of the scree skink haplotypes are currently found in the wild or captive Otago skink populations ( Figure 3), with low to moderate levels of genetic differentiation to the Otago skink haplotypes present in the respective regions (Table 5).

Discussion
Our study indicates that the Otago skink exhibits substantial genetic divergence (9.7-13.4%) from its closest relatives in the    Miocene [34], corresponding to the commencement of tectonic activity along the Alpine Fault and the uplift of the Southern Alps [2,3,7,8]. At present, these six species have non-overlapping distributions in the South Island (Figure 1), which may indicate that the diversification of this lineage during the late Miocene occurred through allopatric speciation [34]. However, both our phylogenetic and phylogeographic datasets indicate that hybridisation has occurred historically (i.e. introgression) between the Otago skink and scree skink in the northern Otago-southern Canterbury region. Incomplete lineage sorting is unlikely to explain our results given the deep divergence between the two species, the restriction of the pattern to the southern end of the scree skinks range, and the ability for the two species to produce viable hybrids in captivity. Hardy [30] hypothesized that the divergence between these two species occurred during the Pleistocene, with the two species retreating to isolated refugia in Otago (Otago skink) and Canterbury (scree skink) during glacial maxima, and expanding their ranges during interglacials. Although the divergence between the two species occurred earlier in the late-Miocene [34], they appear to have hybridised when they came into secondary contact in the northern Otago-southern Canterbury region during Pleistocene interglacials. Our data suggests that the scree skink hybridised with both the eastern and western lineages of the Otago skink. This might have involved Otago skink populations that have since gone locally extinct (i.e. scree skinks were found to have haplotypes that are not currently present in the wild or captive Otago skink populations). Alternatively, the scree skink haplotypes may have diverged since the introgression event. Interestingly, we detected no evidence of scree skink haplotypes occurring in either the eastern or western lineages of the Otago skink. Although the two species do not currently occur in sympatry, they have the capacity to produce viable hybrids in captivity (D. Keall unpublished data, [33]). This not only has important implications for selection of breeding stock for the Otago skink captive breeding program, but also for the conservation management of both species in the wild (e.g. [82]). Future studies should use nuclear markers (e.g. nuclear genes, microsatellites) to further investigate the patterns of introgression between the Otago and scree skink.
Two main lineages were evident within the Otago skink, with the divergence between the populations in eastern Otago and western Otago occurring during the Pliocene (mean 3.7 mya, range 2.8-4.6 mya). Phylogeographic breaks of an equivalent age have been reported in this region for the sympatric grand skink (O. grande, 3.8 mya; [39]), McCann's skink (O. maccanni, 3.6 mya; [45]), green skink (O. chloronoton, 5.3 mya; [40]), and cryptic skink (O. inconspicuum, 2.7 mya; [50]). These divergences in South Island skink species are concordant with the Nevis-Cardrona fault system (marked by the Cardrona and Nevis rivers) that has been active since the Miocene and delineates a topographic boundary between eastern Otago (undulating grassland habitats) and western Otago (deeply-eroded rugged mountainous habitat) (reviewed in [83]). However, the extinction of the Otago skink from this boundary region makes it difficult to identify the exact position of this split. The consistent length variation in the control region haplotypes, and the lack of shared haplotypes between eastern and western  Otago, indicates that there has been no geneflow between the remnant Otago skink populations in the two regions for a substantial period of time. This is exemplified by our finding that the divergence between the eastern and western regions accounts for ,95% of the genetic variation within the Otago skink. However, given the recent extinction of the populations in central Otago it is unknown what pattern of geneflow (e.g. isolation by distance) was evident throughout the species range prior to its decline. Substantial genetic structuring was also evident within both the eastern (three Macraes Flat populations) and western (Lindis Pass, Glenfoyle station) populations of the Otago skink. Despite this, the presence of shared haplotypes within each region indicates that there has been recent geneflow among adjacent populations. This pattern of population structuring is consistent with that observed in the grand skink [39], a species that has similar habitat preferences (i.e. rocky outcrops in montane tussock grassland), and a near identical distribution (both historical and current distribution) and conservation status (Nationally Critical) [20,22,84]. Similar to the recommendations for the grand skink [39], we suggest that the eastern and western populations of the Otago skink should continue to be treated as separate management units (e.g. [85][86][87]). Although a reserve was recently established at Glenfoyle Station and Otago skink survey work is regularly conducted in the western populations, the populations in the Macraes Flat region are subject to the most intensive management by the Department of Conservation [22,84]. Our study indicates that these managed populations (Wildlife, Falcon) contain all of the haplotypes known from the eastern region, but only a small proportion of the species total genetic diversity.
Considerable genetic and haplotypic diversity was evident in the Otago skinks in the captive breeding program. All captive animals had haplotypes that originated from the eastern region. The initial founders of the captive breeding program were sourced from near Sutton in the 1980s-1990s (A. Hutcheon, personal communication), and this accounts for most captive animals having haplotypes (i.e. Haplotypes G, H, I) that originate from outside of the intensively managed populations at Macraes Flat. Given the recent population declines of the Otago skink in the Sutton region [22,26], the captive breeding stock may contain haplotypes that have since gone extinct in the wild. Although the recent additions to the captive breeding program have been sourced from the Macraes Flat region, only two (A, B) of the four haplotypes from the region are represented in the captive populations.
The Department of Conservation is in the process of establishing a separate captive breeding program for the western populations of the Otago skink (A. Hutcheon, personal communication), and continues to implement in-situ conservation management actions at the Glenfoyle station reserve [22]. As these populations represent a distinct management unit, the captive breeding stock should help to safeguard the full extent of genetic diversity that is evident in the species. Indeed, the current recovery plan aims to maximise the genetic diversity present in captivity populations to enable the reintroduction of individuals to the wild [22]. However, as discussed previously by Connolly & Cree [23], the Otago skink breeding program should consider the phenotype, or 'quality' (i.e. morphology, parasite load, body condition, growth rate, sprint speed), of the captive animals, rather than just the genetic diversity represented in the breeding stock.

Conclusions
The Otago skink has experienced a drastic reduction in its distribution over the last ,200 years and now persists in two remnant regions in the periphery of its original distribution [22,26]. A consistent pattern in conservation biology has been that during range decline, most species are able to persist in the peripheral regions of their historical range after the local extinction of the populations from the interior of the geographic distribution [88,89]. Compared to the core of the species range, populations occurring at range margins are usually exposed to suboptimal environmental conditions (reviewed in [90]), and therefore remnant populations may not be representative of the species habitat preferences and population genetic structure (reviewed in [89]). This may be true in the Otago skink, as the haplotypes present in the captive animals and the introgressed scree skinks indicate that there has been a significant recent reduction in genetic diversity present across the species range. It may also indicate that the Otago skink currently persists in suboptimal or atypical habitats, and might benefit from translocations (within predator-proof exclosures) to areas near the core of the species' historic range. Hybridisation with the scree skink is another potential concern for the future conservation management of the Otago skink. Although there is currently limited opportunity for hybridisation between the two species as their ranges do not overlap, as many species are anticipated to shift their distributions as a result of climate change [91], there is the potential for hybridisation in the future if the scree skink expands its range southwards into the distribution of the Otago skink.