Reduced representation sequencing detects only subtle regional structure in a heavily exploited and rapidly recolonizing marine mammal species

Abstract Next‐generation reduced representation sequencing (RRS) approaches show great potential for resolving the structure of wild populations. However, the population structure of species that have shown rapid demographic recovery following severe population bottlenecks may still prove difficult to resolve due to high gene flow between subpopulations. Here, we tested the effectiveness of the RRS method Genotyping‐By‐Sequencing (GBS) for describing the population structure of the New Zealand fur seal (NZFS, Arctocephalus forsteri), a species that was heavily exploited by the 19th century commercial sealing industry and has since rapidly recolonized most of its former range from a few isolated colonies. Using 26,026 neutral single nucleotide polymorphisms (SNPs), we assessed genetic variation within and between NZFS colonies. We identified low levels of population differentiation across the species range (<1% of variation explained by regional differences) suggesting a state of near panmixia. Nonetheless, we observed subtle population substructure between West Coast and Southern East Coast colonies and a weak, but significant (p = 0.01), isolation‐by‐distance pattern among the eight colonies studied. Furthermore, our demographic reconstructions supported severe bottlenecks with potential 10‐fold and 250‐fold declines in response to Polynesian and European hunting, respectively. Finally, we were able to assign individuals treated as unknowns to their regions of origin with high confidence (96%) using our SNP data. Our results indicate that while it may be difficult to detect population structure in species that have experienced rapid recovery, next‐generation markers and methods are powerful tools for resolving fine‐scale structure and informing conservation and management efforts.


| INTRODUC TI ON
Population genetic analyses have become an integral part of conservation studies. These analyses allow the quantification of parameters relevant to endangered populations such as genetic diversity, effective population sizes, gene flow and the inference of past population histories (Beaumont, Zhang, & Balding 2002;Lopes & Boessenkool, 2010). Identifying management units and determining the degree of connectivity among demes are often the first steps in allocating conservation resources and developing management strategies (Moritz, 1999;Palsbøll, Bérubé, & Allendorf, 2007;Schwartz, Luikart, & Waples, 2007). For example, species with a relative lack of population structure that show little variation across their ranges benefit from management of all subpopulations as a single unit (Palsbøll et al., 2007). Conversely, species characterized by genetically-distinct, geographically isolated demes require management of distinct units to maximize maintenance of the species' evolutionary potential (Lowe & Allendorf, 2010;Palsbøll et al., 2007).
While describing population structure is crucial for species management, considering the underlying causes of this structure and discriminating between events of evolutionary significance and those with anthropogenic causes is also critical (e.g., rock wren; Weston & Robertson, 2015). For instance, if genetically distinct populations have diverged gradually as the result of glacial isolation or post-glacial colonization (e.g., kea; Dussex, Wegmann, & Robertson, 2014), it is often considered best to work toward maintaining that existing population structure. Conversely, for species that have undergone severe reductions in population size and genetic population fragmentation (e.g., mohua; Tracy, Wallis, Efford, & Jamieson, 2011), restoration of populations to better reflect the population structure prior to human disturbance may be more beneficial (Moritz, 1994;Palsbøll et al., 2007).
Recent advances in next-generation sequencing (NGS) allow the characterization of large numbers of molecular markers, increasing the resolution of comparisons of individuals within and between populations (Helyar et al., 2011). Single nucleotide polymorphisms (SNPs) have become the markers of choice for an increasing number of modern population genetic studies (Andrews, Good, Miller, Luikart, & Hohenlohe, 2016;Helyar et al., 2011), despite initial speculation on their relative lack of genetic diversity compared to individual microsatellites (Rosenberg, Li, Ward, & Pritchard, 2003).
While 100 SNPs may only be as informative as 10-20 microsatellites (Kalinowski, 2002), tens of thousands of SNPs can be discovered via NGS, even for non-model species (Ellegren, 2014;Seeb et al., 2011). The degree of definition afforded by such large marker datasets therefore statistically outperforms studies in which a few dozen microsatellites were used to describe population structure (Helyar et al., 2011;Liu, Chen, Wang, Oh, & Zhao, 2005).
A growing number of studies have also explored the specific applications of reduced representation sequencing (RRS) methods in evolutionary biology (e.g., RADseq; Andrews et al., 2016). One of these methods, referred to as Genotyping-by-sequencing (GBS; Elshire et al., 2011) has been used in a range of non-model species as a means of inferring genotypes and describing population genomic structure (Johnson et al., 2015;Skovrind et al., 2016). Species for which traditional markers (i.e., microsatellite, single mtDNA sequences) have shown limited resolution in the description of population structure (e.g. Dussex et al., 2016;Robertson & Gemmell, 2005) could benefit from RRS methods, due to the higher definition and overall performance of the resulting marker sets. For instance, species characterized by strong dispersal potential and high gene flow are interesting candidate species for re-evaluations of population structure analyses using RRS.
The New Zealand fur seal (NZFS, Arctocephalus forsteri) has undergone a rapid recolonization following near extinction. Prior to the arrival of humans, NZFS were distributed throughout the east and west coasts of New Zealand's North and South Islands (Smith, 1989(Smith, , 2005. Pre-historic Polynesian settlers hunted NZFS as a source of meat, which resulted in a rapid population decline (Emami-Khoyi et al., 2017;Smith, 1989Smith, , 2005 and the near elimination of the mainland fur seal population (Lalas & Bradshaw, 2001). By the time of European arrival in New Zealand in the late 18th Century, the NZFS mainland breeding range was confined to the south-western South Island (Smith, 1989). Sealing by European settlers for the acquisition of furs had a further devastating impact on NZFS and resulted in a depletion of the NZFS population (Lalas & Bradshaw, 2001). Overall, an estimated 1.5 million skins were harvested between 1792 and 1849 around Australia's southern coast, New Zealand, and at the adjacent subantarctic islands (Ling, 1999), suggesting that the number of NZFS before sealing may have been as high as 1-2 million (Baird, 2011). However, because records of the number of fur seals killed during that time are incomplete, inaccurate, or deliberately misleading to protect commercial interest of the day, it is unclear how many seals may have remained in New Zealand in the 1830s, when sealing stopped due to low profitability (Lalas & Bradshaw, 2001).
Since implementation of regulations in the late 19th century and the official cessation of commercial sealing in 1946 (Sorensen, 1969), the NZFS has rapidly recolonized much of its former range, and has re-established breeding colonies in the South Island of New Zealand (Baird, 2011;Bradshaw, Lalas, & Thompson, 2000;Crawley & Wilson, 1976). However, only one breeding colony, founded in the early 1990s, has been re-established in the North Island (Dix, 1993).
Population declines and fragmentation often lead to population differentiation via the effect of genetic drift when demographic units within a larger panmictic population become fixed for different alleles (Charlesworth, 2009;Frankham, 2005). Alternatively, genetic differentiation can arise via founder effects following the recolonization of a species' former range (Ramachandran et al., 2005;Templeton, 1980). However, previous studies on pinniped species with histories of intense exploitation have shown a lack of (e.g., Arctocephalus gazella, Hoffman, Grant, Forcada, & Phillips, 2011) or only moderate fine-scale structure (e.g., Arctocephalus tropicalis, Wynen et al., 2000;Halichoerus grypus, Fietz et al., 2016). Research based on modern and ancient DNA suggests that, while there was a shift in the genetic composition of NZFS in southern New Zealand resulting from the demographic decline, this decline may not have been long enough to induce a significant loss of genetic diversity as evidenced by comparable levels of genetic diversity in modern and ancient populations . This result is in stark contrast with the severe demographic decline and near range-wide extirpation in the 1800s. Moreover, in spite of potential founder effects resulting from recolonization from a few "refugia" colonies, studies based on microsatellite data have suggested that any population differentiation would have been transient, owing to the high vagility of the species combined with a rapid recolonization time (Figure 1, Dussex et al., 2016).
Studies aimed at assigning NZFS individuals to their colonies or broader regions of origin using microsatellite markers and classical assignment approaches have also had limited success (Dussex et al., 2016;Robertson & Gemmell, 2005). Genetic distinction of potential management units has thus been difficult. A prevailing view from the literature published on NZFS genetics to date is that more numerous and/or variable genetic markers could significantly aid in resolving structure and improving assignment probabilities (Dussex et al., 2016;Robertson & Gemmell, 2005)-using an RRS approach provides the opportunity to empirically test this prediction.
Today, NZFS are the most commonly caught marine mammal in New Zealand commercial trawl net fisheries (Thompson & Abraham, 2010 Thompson & Abraham, 2010) and males accounted for about 45% of reported bycaught individuals (Baird, 2011). Bycatch incidents seem to occur in clusters throughout New Zealand waters, mostly off the northern East (~40%) and West Coasts of the South Island and around offshore islands (Thompson & Abraham, 2010).
While it is possible that bycatch may affect some colonies more than others, there is a high degree of uncertainty on the localized impact F I G U R E 1 Extant New Zealand fur seal (A. forsteri) breeding colonies with full and empty dots depicting colonies sampled and not sampled in this study, respectively, with the number of samples analysed in this study and the year of first breeding from King (2005). Cape Palliser (CP) is currently the only known breeding colony in the North Island. The dashed line represents the genetic clusters based on Dussex et al. (2016) and Salis et al. (2016) microsatellite and mtDNA results, and the dotted line represents the approximate limit of the admixture zone between the two main clusters. Arrows depict the post-sealing recolonization scenario from hypothetical south-western South Island refugia depicted with a star inferred by ABC analyses (Dussex et al., 2016;this study) of commercial fisheries on the species. West Coast colonies in particular face the challenges of lower habitat accessibility (Smith, 1989(Smith, , 2005, climatic shifts that affect food availability  and potentially higher conflict with fisheries (Ministry for Primary Industries, 2014). Moreover, the possibility of sexual segregation and niche divergence as seen in the country's other native pinniped, the New Zealand sea lion (Phocarctos hookeri) (Leung, Chilvers, Nakagawa, Moore, & Robertson, 2012), could result in sexbiased mortality which would then negatively impact population dynamics. This means that in spite of an overall trend in steady growth, these effects may negatively affect colony viability. The ability to assign bycaught individuals to their colonies of origin could thus help identify the colonies most at risk of decline.
In this study, we use SNP markers discovered via GBS to examine the population structure in NZFS and compare the markers' performance with the traditional markers used in previous studies on the same species (Dussex et al., 2016;Lento, Haddon, Chambers, & Baker, 1997;Robertson & Gemmell, 2005). We also assess the power of GBS-generated SNPs to assign potential bycatch to their colony by estimating the proportion of pups that can be correctly assigned to their known colony and/or regions of origin using a GBS SNP dataset. We hypothesize that in spite of the high individual vagility and rapid recolonization potential of the species, GBS-generated SNPs will allow the detection of fine-scale population differentiation among colonies. Moreover, we hypothesize that the SNPs presented here will outperform previously used microsatellite markers and produce more confident assignment probabilities.

| Sampling
Between 1998 and 2003, we obtained 169 NZFS skin samples from eight mainland NZFS breeding colonies ( Figure 1). Only pups were sampled because they are a direct representation of each breeding colony, and to avoid bias caused by seasonal and yearly variation in adult dispersal. For each pup sampled, a small piece of skin was taken from the tip of a digit on a hind flipper using piglet ear notch pliers as per an established protocol (Majluf & Goebel, 1992) and stored in 70% ethanol.

| DNA extraction, library preparation, and sequencing
DNA was extracted from fur seal skin samples using a modified version of a lithium chloride DNA extraction protocol (Gemmell & Akayima, 1996). DNA quality and concentration were evaluated using both a Nanodrop 2000c spectrophotometer (Thermo Fischer Scientific) and a Fluorometric Quantitation Qubit system (Thermo Fischer Scientific). For library preparation, we used GBS as a sequencing preparation protocol (Elshire et al., 2011).
We used two 96-well plates for library construction. Prior to adapter ligation, an automated picogreen assay (Thermo Fisher Scientific) was performed to normalize the DNA concentration in all samples to 20 ng/μl. The restriction enzyme Pst1 was used due to its empirical success for sequencing of vertebrates and its capacity to produce libraries with high depth of coverage (De Donato, Peters, Mitchell, Hussain, & Imumorin, 2013). 17 μl H20, 2 μl CutSmart buffer, 0.5 μl NEB buffer, 0.5 μl Pst1 were added to wells containing DNA and adapters. A solution containing ligase buffer, ATP, and T4 ligase was then added to each well to aid in adhering adapters to sticky ends of the now-cut fragments. This combination of DNA, enzyme, adapters, and ligation solution was incubated at 22°C for 1 hr, and the temperature was then increased to 65°C for 30 min to inactivate T4 ligase (Elshire et al., 2011). 5 μl from each well were pooled together into one solution containing all fragments, each of which was identifiable by its respective barcode adapter. Fragments were amplified in three subsequent polymerase chain reactions and fragment length (in base pairs) was examined on a bioanalyser (Agilent 2100 bioanalyser, Santa Clara, CA, USA).
The amplified library was purified using a Pippin Prep (SAGE Science, Beverly, MA, USA) to select the DNA sequencing library in the size range of 150-500 bp. Finally, single-end sequencing (1 × 100) was performed on an Illumina HiSeq2500 using v4 chemistry, which yielded approximately 25 Gb of raw sequence data per lane.

| SNP calling and filtering
SNP calling was performed in Stacks v. 1.30 (Catchen et al., 2011). We used "process_radtags" to demultiplex the data and discard low quality data with a probability of sequencing error >0.10% (Phred score = 10). SNP loci were then called using the Stacks "denovo_map.pl" pipeline (i.e., without a reference genome), using a minimum stack depth (-m) of 5 and default parameters for mismatches between loci within samples (-M = 2) and when building the catalog (-n = 1). The "populations" program in Stacks was used to filter out monomorphic loci and sort individuals into separate files for each colony. 180,303 biallelic SNPs remained after this filtering. PLINK 1.9 (Purcell et al., 2007) was then used to exclude individuals with >90% missing data (-mind filter), SNPs with >20% missing call rates (-geno filter), and SNPs with minor allele frequency <5% (-maf filter).
SNPs significantly deviating (exact test) from Hardy-Weinberg equilibrium (HWE) (-hardy filter) and pairs of SNPs in linkage disequilibrium (LD) (-ld filter) in at least six of our eight populations were removed. For downstream analyses, we used the whole set of SNPs generated in PLINK, but also generated a separate subset of 5,000 and 1,000 random SNPs for computationally intensive analyses (e.g., STRUCTURE, DIYABC). Finally, the resulting. Ped files were converted into other file formats with the program PGDSPIDER 2.0 (Lischer & Excoffier, 2012), enabling their input into relevant population genetic analysis software.

| Genetic diversity and population structure
For our analyses of population structure, we evaluated genetic variation at the colony and genetic cluster levels. Three main clusters were defined based on previous NZFS genetic studies (Dussex et al., 2016;Salis et al., 2016)  Genetic diversity (H o , H e and F IS ) was estimated for each colony and region using Arlequin v. 3.5 (Excoffier & Lischer, 2010).
Pairwise F ST was calculated among the eight colonies in Arlequin v. 3.5 using Slatkin's linearized F ST (Slatkin, 1995) and tested for Isolation by distance (IBD) using a Mantel test implemented in GenAlEx (Peakall & Smouse, 2006). We then tested for correlations between estimates of H e and F ST derived from microsatellite loci by Dussex et al. (2016)  We used the Bayesian cluster methodology implemented in STRUCTURE 2.2 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) to make inferences regarding the number of population clusters in our sample of NZFS using a random subset of 5,000 SNPs generated with PLINK 1.9 (Purcell et al., 2007). We used an admixture model with correlated allele frequencies and performed 10 iterations (chain length 500,000 steps, burn-in = 200,000 steps) for each K (from 1 to 8), using the LOCPRIOR option to assist with clustering of faint-structure data (Hubisz, Falush, Stephens, & Pritchard, 2009). The number of distinct genetic clusters were inferred using the LnP(K) value of STRUCTURE 2.2 and the ∆K (Evanno, Regnaut, & Goudet, 2005) method implemented in STRUCTURE HARVESTER (Earl & vonHoldt, 2012). CLUMPAK (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015) was then used to estimate individuals' assignment coefficient (q) to each genetic cluster and to visualize the results. As a comparison, we also used the program ADMIXTURE (Alexander, Novembre, & Lange, 2009) to estimate the number of clusters for each K (from 1-8). This  . DAPC has been shown to perform as well as STRUCTURE and has the advantage of being unconstrained by the assumptions of HWE (Jombart et al., 2010). DAPC is similar to principal component analysis, but unlike PCA, which maximizes the total variation in the dataset, DAPC maximizes the variation among different groups and minimizes variation within groups (Jombart et al., 2010). This analysis was performed with prior information on individual colonies.

| Demographic history
Previous studies have established several aspects of NZFS population history using Approximate Bayesian Computation (ABC; Beaumont, Zhang, & Balding 2002) approaches with microsatellite (Dussex et al., 2016) and mitogenome data (Emami-Khoyi et al., 2017). Here, we did not attempt to reiterate the demographic model comparisons described in these studies, but rather to estimate demographic parameters based on SNP data for the most likely scenarios inferred in these studies. We thus designed a scenarios including  Figure S1, Table 1). We also designed a null model of constant effective population size though time followed by the same population divergence. The generation time for NZFS is estimated at 4-6 years in females and at 5-9 years in males (Dickie & Dawson, 2003). Thus, we assumed an average 7 year generation time for NZFS.
Because DIYABC is a computationally intensive method, analyses were performed using a subset of 1,000 randomly selected SNPs from our complete dataset. We generated 10 6 simulations for each scenario. As summary statistics, for single sample statistics we used the mean gene diversity across polymorphic loci and mean gene diversity across all loci (Nei, 1987). For two sample statistics, we used mean of nonnull F ST distances between the two samples across loci, the mean of F ST distances between the two samples across all loci, mean of nonnull Nei's distances between the two samples across loci, and mean of Nei's distances between the two samples across all loci (Nei, 1972). We used the standard Hudson's (2002) algorithm for selection of minimum allele frequency (MAF).
The posterior probability of each scenario was estimated using both the direct and logistic regression approaches (Fagundes et al., 2007). The ten thousand datasets (1%) with the smallest Euclidean distances were then retained to build posterior parameter distribution. To check model performance, we first empirically evaluated the power of the model to discriminate among scenarios (confidence in scenario choice). The approach consists of simulating pseudo-observed datasets with parameters drawn from the posterior parameter distribution of the considered scenario and positioning the summary statistics of the observed data in the summary statistic distribution of these pseudo-observed data. The scenario is then considered suitable if the observed data summary statistics are included in the confidence interval drawn from pseudo-observed data. We calculated statistical measures of performance and Type I and Type II error rates as a means of model checking (Cornuet, Ravigné, & Estoup, 2010;Excoffier, Estoup, & Cornuet, 2005).

| Outlier loci detection
To maximize assignment power to the colonies sampled, we attempted to identify colony-specific and genetic clusters-specific sets of outlier loci that showed high differentiation between colonies and genetic clusters. These analyses were performed for each of the eight colonies and for the three genetic clusters.
To correct for multiple testing, we transformed the p-values obtained to q-values (i.e., the false discovery rate FDR, analog to the p-value; Benjamini & Hochberg, 1995). Secondly, we used a Bayesian simulation-based test (Beaumont & Balding, 2004) that has been further refined and implemented in the software Bayescan 2.0 (Foll & Gaggiotti, 2008). We based our analyses on 10-pilot runs each consisting of 5,000 iterations, followed by 100,000 iterations with a burn-in of 50,000 iterations. We used a FDR of 0.05 as the threshold for outlier locus detection in this test.

| Population assignments
To assess the power of GBS-generated SNPs to assign potential NZFS bycatch to their colony or genetic cluster, we treated pups as unknown in origin and attempted to assign them back to their home colonies and genetic clusters. We first performed assignments of pups to the eight colonies using all 26,026 SNPs, a random subset of 5,000 SNPs and the 74 colony-specific outlier loci. Secondly, we performed assignments of pups to the three genetic clusters using all 26,026 SNPs, the same random subset of 5,000 SNPs and 18 genetic cluster-specific outlier loci (see Results). This made for a total of six analyses.
We used Genodive 2.0b23 (Meirmans & Van Tienderen, 2004) to assign pups to their colony or cluster. The program uses the leaveone-out (LOO) validation procedure in which an individual to be assigned is removed from its source population and treated as an individual of unknown population origin before calculation of the population allele frequency. The purpose of the LOO procedure is to remove a bias present when allele frequencies are calculated from the same individuals that are subsequently assigned. We used the home likelihood criterion (the likelihood that an individual comes from the population where it was sampled) because it is more appropriate when only part of all possible source populations have been sampled (Meirmans & Van Tienderen, 2004). We used a significance threshold of 0.05 and replaced zero frequencies by 0.005, as suggested by Meirmans and Van Tienderen (2004) and generated 10,000 permutations (i.e., number of datasets). Missing values were replaced by randomly picking alleles from the global allele pool.

| Data filtering
Out of a total of 169 genotyped individuals, we excluded two individuals for having more than 90% of loci missing. Out of a total 180,303 biallelic SNPs, we excluded 97,339 SNPs for being missing in at least 80% of individuals, and 56,937 SNPs with <5% minor allele frequency. We further excluded a total of 3,834 SNPs because they were either out of HWE, or in LD. After applying these filters, our dataset consisted of 167 individuals with data for 26,026 SNPs. TA B L E 1 Prior and posterior distributions of parameters for a demographic scenario including a Polynesian-and European-induced decline followed by recolonization and population divergence that obtained the highest posterior probability (PP = 1) when comparing two demographic scenarios in New Zealand fur seal (A. forsteri). Timing of events is in generations, assuming a generation time of 7 years (Dickie & Dawson, 2003). Type I and II errors were of 0.04 and 0.05, respectively

| Genetic diversity and population structure
Measures of heterozygosity were similar across colonies and regions (

| Demographic history
Demographic reconstructions using the Approximate Bayesian Computation (ABC) approach strongly supported a scenario of Polynesian-and European-induced declines followed by range-wide recolonization and population divergence with a posterior probability of 1. Our estimates supported severe bottlenecks with potential 10-fold and 250-fold declines in response to Polynesian and European hunting, respectively (Table 1). The observed data summary statistics for the preferred model were included in the confidence interval drawn from pseudo-observed data. Type I and II errors were of 0.04 and 0.05, respectively.

| Outlier loci detection
Using the Fdist2 method implemented in LOSITAN, we did not find any outlier loci under putative selection. The Bayescan method

| Population assignments
We attempted to assign pups treated as unknowns to colonies and identified clusters using the LOO test. Using either all 26,026 SNPs or the 5,000 randomly selected SNPs, we assigned 167 pups to their colonies of origin with, on average, 44.2% and 42.1% accuracy, respectively, though low assignment probabilities to NZ-North-East colonies CP (14%-21%) and OP (0%-7%) drove these averages down (Table 4a). We observed the highest assignment probabilities at the colony level for OBI with 96% and 93% correct assignments for the 26,026 SNPs and 5,000 SNPs, respectively. Using outlier loci increased our colony-level assignment probabilities by around 10%. Average assignment probabilities to colonies using 74 colonyspecific outlier loci where higher on average (53.3%) when considering all colonies and even higher (69%) when excluding CP and OP, which both showed the lowest assignment probabilities (Table 4a).
Individual assignment probabilities of pups to the three genetic clusters (NZ-North-West, NZ-North-East and NZ-South) were on average higher (64%-69%; Table 4b) than for assignments at the colony level. As expected, and in line with assignments probabilities to colonies, assignment probabilities were higher for the NZ-North-West (80%-97%) and NZ-South (77%-95%) clusters than for the highlyadmixed NZ-North-East cluster (13%-34%; Table 4b) across all SNP datasets. However, the assignment probability to the admixed colonies was highest (34%) when using the 18 outlier loci (Table 4b).
F I G U R E 3 DAPC for 167 New Zealand fur seal (A. forsteri) from eight breeding colonies based on 26,026 SNP loci F I G U R E 2 Individual clustering assignment in New Zealand fur seal (A. forsteri) (n = 167) from eight breeding colonies for the most likely number of clusters (K = 2) using 5,000 random SNPs in STRUCTURE

| Subtle population structure and near panmixia
Our panel of 26,026 SNPs illustrates the potential of RRS methods to identify previously undetected sources of variation in nonmodel species. Overall, our data strongly support a pattern of near panmixia in NZFS, but pairwise F ST values among range-wide colonies supported a weak but significant west to east genetic gradient. Bayesian clustering and DAPC results support the existence of two main genetic clusters (NZ-North-West and NZ-South), and further suggest that NZ-North-East colonies (CP and OP) are sites of genetic admixture between these two main clusters.
This observed structure was further supported by constructing a PCA from a matrix of estimated genomic relationships (Patterson, Price, & Reich, 2006), which in turn were calculated following the bioinformatics and statistical methods in Dodds et al. (2015). This approach allows the PCA to be formed when there are missing genotype calls and without the use of stringent filtering as the method directly accounts for the read depth and missing observations and thus gives further support to our results.

| Performance of SNPs loci
H e and H o values did not differ substantially among colonies (0.225-0.250), which is consistent with previous results based on microsatellite data (Dussex et al., 2016;Robertson & Gemmell, 2005) with H e and H o between (0.67-0.79). Higher estimates of H e derived from microsatellite markers may be a consequence of ascertainment bias caused by selecting the most polymorphic markers (Haasl & Payseur, 2011;Queirós et al., 2015). However, there was no significant correlation between H e or F ST estimates derived from microsatellite and the SNPs generated here. The lack of correlation between H e estimates could be due to low sample size (e.g., Fischer et al., 2017) while for pairwise F ST , it could also be due to an ascertainment bias in the selection of microsatellite loci (Haasl & Payseur, 2011;Ryynänen, Tonteri, Vasemägi, & Primmer, 2007).
TA B L E 4 Proportions of New Zealand fur seal (A. forsteri) pups correctly assigned to (a) their colonies or (b) genetic clusters inferred from Bayesian clustering analyses. Assignments were performed independently using all 26,026 SNPs, 5,000 randomly selected SNPs, 74 outlier loci identified at the colony level, and 18 outlier loci identified at the genetic cluster level In spite of this very subtle population structure, assignment testing results indicate that large SNP marker sets allow NZFS pups to be assigned to their geographic regions of origin with high confidence, but less so to specific colonies. The proportion of pups correctly assigned to genetic clusters (64%-69%) and to breeding colonies (42%-53%) were on average higher for all SNP datasets used here than for microsatellite markers (59.3% and 22.7% correctly assigned to clusters and colonies, respectively) (Dussex et al., 2016;Robertson & Gemmell, 2005). Furthermore, a random subset of 5,000 SNPs was roughly as effective as the entire 26,026 SNP dataset in achieving these assignments, while using outlier loci in some instances increased the proportion of pups correctly assigned to their colony or genetic cluster of origin. As expected, the proportion of correct assignments was lowest for the NZ-North-East cluster and its colonies (CP and OP), which is in line with the identification of a zone admixture between the NZ-North-West and NZ-South clusters.
Our higher assignment probabilities can likely be attributed to the higher statistical power of our dataset (Helyar et al., 2011;Morin, Luikart, & Wayne, 2004;Morin et al., 2009), although we were still unable to achieve confident assignment to specific colonies. One exception is that of OBI, which also had one if the highest proportion of correct assignments with microsatellite data (68.8%; Dussex et al., 2016). This result is in line with the suggestion that OBI may represent a refugium from sealing (Dussex et al., 2016;Salis et al., 2016;Smith, 2005) and thus a colony with a higher proportion of private alleles, thus increasing the proportion of pups correctly assigned. While colony-level assignment for species may be unfeasible with the tools presently available, our results suggest that outlier loci assignment techniques could enable confident assignment in the future (e.g., Nielsen et al., 2012). However, identifying highly differentiated outlier loci will be crucial for accurate assignments as genetic structure is predicted to become rapidly eroded due to the high vagility of the species.
As a whole, our measures of intra and interspecific population structure and assignments are consistent with what might be expected from the severe historic bottleneck and rapid range-wide recolonization NZFS experienced, a model which has been confirmed through Bayesian reconstruction here and in previous studies (Dussex et al., 2016;Emami-Khoyi et al., 2016).

| Demographic recovery and transient population structure
The subtle population structuring seen here for NZFS is most likely due to the high vagility of the species and its recolonization history (Dussex et al., 2016;Emami-Khoyi et al., 2017). Demographic reconstructions strongly supported a scenario of Polynesian-and European-induced declines followed by range-wide recolonization and population divergence with a posterior probability of 1, which is consistent with archaeological records (Smith, 1989) and the recent exploitation history of the species (Lalas & Bradshaw, 2001 Polynesian N e of ~570,000 and ~4,900,000, which is one to two orders of magnitude higher than previous estimates of historical N e (Dussex et al., 2016;Emami-Khoyi et al., 2017). Our contemporary N e estimates seem more realistic with an N e of 600-1,000 predicted following demographic recovery after cessation of sealing.
The estimate for the post-European surviving population seems extremely low, yet not unlikely based on the history of near extermination of this and other fur seal species in the 1830s (Lalas & Bradshaw, 2001). For instance, only 5 NZFS were reported in the Bounty Islands in the 1830s (Taylor 1982 (Crawley, 1990;Wilson, 1981). It is thus possible that gene flow from unsampled colonies in Foveaux Strait, or subantarctic islands, may have rapidly inflated the N e of modern populations.
While theory predicts that demographic decline can result in population structuring via the effects of genetic drift (Charlesworth, 2009;Frankham, 2005), NZFS as well as other otariids do not show marked population structuring (Dussex et al., 2016;Hoffman et al., 2011). This could be due to the fact that the demographic bottleneck in NZFS was too short in duration to induce a significant loss of genetic diversity and strong population differentiation via genetic drift and to a high dispersal rate contributing to range wide gene flow and homogenization of allelic distribution. NZFS most likely survived in refugia post-exploitation, and quickly recolonized after sealing became unprofitable (Dussex et al., 2016;Emami-Khoyi et al., 2017).
Therefore, the lack of marked population structuring is most likely due to the high vagility of the species and its recolonization history.
It is even possible that the subtle structure detected here is transient and that the whole species may represent a panmictic unit within a few generations because population size has increased dramatically in the years since official the cessation of seal hunting (Baird, 2011). Consequently, genetic assignments may become more difficult in the future and may require more traditional approaches, such as tagging and mark-recapture to assign bycatch to colonies (e.g., Bradshaw, Davis, Lalas, & Harcourt, 2000).

| CON CLUS ION
To date, few studies have compared the performance of microsatellite and SNPs for highly vagile species. Overall, while the genetic structure described here was not substantially different from that examined with microsatellite loci, a larger number of markers allowed for better assignments of pups to their genetic cluster and colony of origin. This was especially true for colonies that were outside of the zone of admixture.
Our results, however, suggest that because of the ascertainment bias associated with the selection of microsatellite loci, overestimation of population structuring is possible when using such markers. SNP loci are thus not only less biased but also offer more statistical power. Nevertheless, due to the species' recolonization history and high vagility, describing finer-scale structuring remains challenging.

DATA A R C H I V I N G S TAT E M E N T
Raw sequence data, processed data and scripts are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.63nm407.