The genetic legacy of extreme exploitation in a polar vertebrate

Understanding the effects of human exploitation on the genetic composition of wild populations is important for predicting species persistence and adaptive potential. We therefore investigated the genetic legacy of large-scale commercial harvesting by reconstructing, on a global scale, the recent demographic history of the Antarctic fur seal (Arctocephalus gazella), a species that was hunted to the brink of extinction by 18th and 19th century sealers. Molecular genetic data from over 2,000 individuals sampled from all eight major breeding locations across the species’ circumpolar geographic distribution, show that at least four relict populations around Antarctica survived commercial hunting. Coalescent simulations suggest that all of these populations experienced severe bottlenecks down to effective population sizes of around 150–200. Nevertheless, comparably high levels of neutral genetic variability were retained as these declines are unlikely to have been strong enough to deplete allelic richness by more than around 15%. These findings suggest that even dramatic short-term declines need not necessarily result in major losses of diversity, and explain the apparent contradiction between the high genetic diversity of this species and its extreme exploitation history.

www.nature.com/scientificreports www.nature.com/scientificreports/ species that were heavily hunted in the past. In particular, recent severe demographic declines or bottlenecks 15,16 can be inferred by comparing the observed genetic diversity of a contemporary sample with the diversity expected under alternative historical scenarios simulated based on the coalescent [17][18][19] . Applied in a comparative context, this approach has shown that the intensity of recent bottlenecks can be influenced by species-specific traits such as breeding ecology and mating system variation 20 . However, our understanding of within-species variation in demographic histories remains limited because most studies of individual species are conducted with suboptimal sample sizes of both individuals and loci, and test for the presence or absence of bottlenecks rather than quantifying their intensity 21 . In order to understand species-wide responses to anthropogenic exploitation, studies are needed that combine exhaustive geographical sampling with modern computational methods capable of deriving quantitative estimates of demographic parameters, such as approximate Bayesian computation (ABC) 22,23 .
The Antarctic fur seal (Arctocephalus gazella) provides an excellent case to evaluate the genetic consequences of exploitation. This pinniped species is polygynous 24 , highly site faithful 25,26 and has a generation time of around ten years 27 . Despite breeding on all of the major sub-Antarctic islands (Fig. 1), it was driven to the brink of extinction by the 18 th and 19 th century sealing industry 28,29 . Harvesting begun shortly after South Georgia was discovered by Captain James Cook in 1775 and 'reckless extermination was the only method of seal-hunting (...) so that the first in the field at a new sealing ground was sure of an immense booty, and late-comers as likely as not would go empty away' 30 . Sealing reached its peak at South Georgia in 1800 when over 100,000 seals were harvested 31 . Shortly afterwards, with the discovery of the other sub-Antarctic islands, a wave of sealing spread across the entire sub-Antarctic region (Fig. 2), culminating in an estimated 1.7 million Antarctic fur seals being taken 32 .
By the turn of the 20 th century, Antarctic fur seals were considered all but extinct and the species was placed under legal protection. For the best part of three decades and despite several dedicated sub-Antarctic expeditions, very few individuals were sighted at the traditional breeding grounds 28,[33][34][35] with the exception of Bouvetøya where a small population of an uncertain fur seal species was observed in the late 1920s 36 . However, the Discovery expedition of 1936 observed a small breeding population of Antarctic fur seals with 12 pups at South Georgia 37 , which grew rapidly in the 1960s and 1970s to reach an estimated five million individuals by 1999 38 . As breeding www.nature.com/scientificreports www.nature.com/scientificreports/ populations were not observed at the other sub-Antarctic islands until the 1950s to the 1980s [39][40][41][42][43][44] , several authors have speculated that the former geographic distribution of this species was likely recolonized from South Georgia 28,38,[45][46][47] . However, a number of studies have reported genetic differences among some of these populations [48][49][50][51] , pointing towards a more complex picture in which more than one relict population may have survived sealing.
In this study, we genotyped over 2,000 Antarctic fur seals at 39 microsatellite loci to investigate the genetic legacy of extreme exploitation in a marine mammal. In order to determine how the species as a whole was affected by sealing, samples from all major breeding locations were analysed to determine the most likely number of relict populations and ABC was used to reconstruct their recent demographic histories. Contrary to expectations, evidence was found for at least four genetically distinct populations, all of which experienced severe demographic declines. Nevertheless, high levels of neutral genetic variability were retained in comparison to other otariid species. Taken together, our results support the assertion that only the strongest bottlenecks lead to major losses of diversity and highlight the importance of relict populations to species recoveries and the maintenance of genetic variation.

Results
Genetic data were generated from 2,000 Antarctic fur seals sampled from eight populations covering the circumpolar distribution of the species (see Methods for details). The 39 microsatellite loci carried between 2 and 26 alleles (mean = 11.44) and none deviated significantly from Hardy-Weinberg equilibrium in more than two out of eight sampling localities after correction for multiple tests (see Supplementary methods and results). The genotyping error rate determined by re-genotyping 96 samples was low (0.003% per allele or 0.005% per reaction). www.nature.com/scientificreports www.nature.com/scientificreports/ population structure. Bayesian analyses of the microsatellite dataset using STRUCTURE produced consistent results across different clustering solutions ( Supplementary Fig. S1). ∆k, which provides an indication of the uppermost hierarchical level of structure, peaked at k = four, while the average log likelihood of the data also gradually increased until k = four, after which a plateau was reached ( Supplementary Fig. S2). The clustering solution for k = four resolved the South Shetland Islands, South Georgia and Bouvetøya as three distinct genetic clusters and Kerguelen, Heard and Macquarie Islands as a fourth cluster, while animals from Marion and Crozet Islands showed mixed membership to the third and fourth clusters (Fig. 1). Further increases in k only introduced additional admixture ( Supplementary Fig. S1). PCA revealed a similar pattern, with Bouvetøya separating from Kerguelen, Heard and Macquarie Islands along PC1 and PC2, Marion and Crozet Islands occupying an intermediate position, and the South Shetlands and South Georgia separating along PC3 ( Supplementary  Fig. S3). Pairwise F st values between populations were significant for all but two comparisons (Supplementary  Table S1). Consequently, our genetic data are consistent with four relict populations comprising (i) the South Shetland Islands, (ii) South Georgia, (iii) Bouvetøya; and (iv) Kerguelen, Heard and Macquarie Islands, plus a further genetically admixed population (v) comprising Marion and Crozet Islands. Subsequent analyses therefore focused on these five populations, although similar results were obtained when sampling locations rather than genetic clusters were analysed ( Supplementary Fig. S4).
Genetic diversity. Allelic richness (A r ), standardised across species as the average number of alleles per ten individuals, was high relative to other otariids, with the Antarctic fur seal being ranked fourth highest out of 13 species (Supplementary Fig. S5) for which comparable microsatellite data were available 20 . In relative terms, genetic diversity did not vary a great deal across the geographical distribution of the Antarctic fur seal (Fig. 3a Bottleneck inference. Two complementary coalescent-based approaches were used to infer the extent to which these populations experienced recent demographic declines. First, the proportion of loci with heterozygosity-excess (prop het-exc) provides an indicator of recent bottlenecks because rare alleles are preferentially www.nature.com/scientificreports www.nature.com/scientificreports/ lost during bottlenecks but their loss has little impact on heterozygosity 15 . Recent bottlenecks therefore generate a transient excess of heterozygosity relative to a population at equilibrium with an equivalent number of alleles, which can be quantified through coalescent simulation. A significant excess of heterozygosity relative to expectations was detected in all five genetic populations identified by STRUCTURE (Table 1 and Supplementary  Table S1). The magnitude of prop het-exc was also similar across populations ( Fig. 3b) suggestive of more or less equally strong declines.
Second, ABC was used to select between a bottleneck and a non-bottleneck model as well as to estimate posterior distributions of the model parameters. ABC was able to distinguish between the two models, with 87% of the simulations being classified correctly under the bottleneck model and 94% of the simulations being classified correctly under the non-bottleneck model ( Supplementary Fig. S6). Posterior probabilities supported the bottleneck model in all five populations and the goodness of fit to the empirical data was consistently better for the bottleneck model (Supplementary Table S3). As another indicator of model quality, posterior predictive checks 23,52 showed that the bottleneck model was able to reproduce the relevant observed summary statistics across all populations ( Supplementary Fig. S7).
As a further evaluation step, we used cross-validation to estimate prediction errors for all model parameters. Cross-validation in ABC uses randomly selected simulations and treats them as pseudo-observed datasets. For each of these simulations, the ABC model was used to estimate a set of summery statistics. Because we also know the true summary statistics, we can calculate the prediction error as where θ i is the estimated parameter value and θ i is the true parameter value of the i th simulated data set. If the estimated values are closer to the true values than expected by chance (i.e. they contain information about the true values), the prediction error will be below one. Under the bottleneck model, the prediction error was 0.457 for the bottleneck effective population size (N e bot), indicating that posterior estimates for this parameter contain information about the underlying true values ( Supplementary Fig. S8). Although in principle the posterior estimates for N e hist also had an acceptable prediction error (0.611), a plot of the cross-validation results ( Supplementary Fig. S8) revealed a systematic underestimation of the true values, which is why the posteriors should be interpreted with caution. By contrast, prediction errors for the other parameters were substantially higher (Supplementary Table S4), with visual inspection of the cross-validation results revealing a large amount of scatter and a poor fit of the predicted to the true values. Consequently, subsequent analyses focused on the posteriors for N e bot, while the N e hist estimates are also shown in Supplementary Fig. S9.
The N e bot estimates were similar in magnitude across populations (Fig. 3c), which is again suggestive of recent bottlenecks of comparable strength. However, the posterior distribution of N e bot was skewed towards slightly smaller values (maximum density ~150) in the eastern population comprising Kerguelen, Heard and Macquarie Islands, whereas maximum densities were closer to ~200 elsewhere. Similar results were obtained when the sampling locations were analysed separately ( Supplementary Fig. S4). The N e hist posterior estimates were very broad (Fig. S9), suggesting that historical effective population sizes ranging between a few thousand and tens-of-thousands of individuals could have led to similar observed genetic diversities under the bottleneck model. expected loss of genetic diversity. The expected loss of neutral genetic diversity caused by a bottleneck with an effective size N e bot of 200 lasting for ten generations was evaluated using coalescent simulations in fastsimcoal 53 . As the long-term N e prior to the bottleneck will affect overall levels of diversity, three scenarios were analysed with N e hist simulated as 1,000, 10,000 and 50,000 respectively. Despite a severe and relatively long-lasting simulated bottleneck, the expected loss of microsatellite alleles was small, ranging from 6% with an N e hist of 1,000 to 18% with an N e hist of 50,000 (Fig. 4). Under the scenario that produced levels of genetic diversity most similar to the empirical values (N e hist = 10,000) only around 14% of alleles were lost. Furthermore, simulated bottlenecks of varying strength and duration, again with an N e hist of 10,000, suggest that only a combination of very strong (i.e. N e bot < 50) and long (i.e. 15-20 generations) bottlenecks would result in genetic diversity losses in excess of 50%. Consequently, the Antarctic fur seal is unlikely to have lost a substantial proportion of neutral genetic diversity despite an inferred steep decline in N e due to overhunting. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
We investigated the genetic consequences of historical commercial exploitation in the Antarctic fur seal. Information on the number of sealing vessels recorded as having visited each of the main sealing grounds allowed us both to formulate realistic priors for our ABC analysis and to interpret our results in the light of historical evidence. An innovation of this study was to use a comparative approach to shed light on the genetic diversity and recent demographic histories of multiple populations within a single species.
Consistent with previous studies of Antarctic fur seals, we found evidence for strong population structure [47][48][49][50]54 . However, our larger sample size of individuals and exhaustive geographical coverage yields a more comprehensive global picture. On the one hand, the discovery of four genetically distinct clusters builds upon recent studies suggesting that relict populations likely persisted in localities other than South Georgia 49,50 . This is in contrast to the previously held belief that Antarctic fur seals only survived in one relict population off South Georgia 28,38,45,46 and suggests that, at least in some regards, the species as a whole may have been more resilient to commercial exploitation than previously believed. Nevertheless, our data also suggest that at least one and possibly several populations may have gone locally extinct. In particular, Macquarie, Heard and Kerguelen Islands showed little in the way of population structure, consistent with the suggestion that this species was locally extirpated at Macquarie and Heard Islands, leaving these sites to be subsequently recolonised by individuals mainly from the Kerguelen Archipelago 47,55 . In addition, a strong signal of admixture at Marion and Crozet Islands implies that these two populations may have either gone extinct or been reduced to sufficiently small numbers for the local gene pool to be swamped by emigrants from nearby Bouvetøya and the Kerguelen Islands. This is supported by anecdotal evidence suggesting that Antarctic fur seals were not present after sealing ceased at Marion Island in the early 1930s 56 .
Two approaches were used to infer the strength of recent population declines. Heterozygosity excess provides a measure of the relative strength of population decline, whereas our ABC analysis specifically tested for www.nature.com/scientificreports www.nature.com/scientificreports/ an absolute reduction in N e (i.e. to below 600). By focusing on time priors that reflected the known sealing history, the ABC analysis also evaluated a clear hypothesis-that 18 th and 19 th century commercial sealing resulted in strong bottlenecks. Although in principle, similar patterns of genetic diversity might be produced by more ancient bottlenecks, these are unlikely to be detected reliably using microsatellite data when subsequent recovery occurred 57 . This is in keeping with a previous study in which a simple model incorporating a recent bottleneck was better supported than a more complex model that also included a small population size during the last glacial maximum followed by expansion 20 .
In practice, heterozygosity excess and ABC produced consistent results, with all of the populations identified by STRUCTURE showing a significant excess of heterozygosity as well as support for recent demographic declines in the ABC analysis. Similar results were also obtained after analyzing each of the sampling localities separately, although Heard Island failed to reach significance in the heterozygosity excess test, possibly due to our relatively small sample size of individuals from this locality. N e bot estimates from the ABC analysis were also quite similar across populations, although the mode of the posterior distribution was slightly lower for the easternmost population. This subtle difference in inferred bottleneck intensity is a reflection of the lower genetic diversity of the Kerguelen Islands, Heard and Macquarie Island populations. It is furthermore consistent with the historical sealing data (Fig. 2), which show that sealing effort was higher in these areas, with Kerguelen and Heard Islands in particular having been visited by more ships in total than any of the other main sealing grounds.
Overall, however, we were struck by the similarity of the N e bot estimates, whose maximal densities consistently fell within the range of around 150 to 200. This provides a good indication that census population sizes on all of the islands may have been similarly diminished, although the relationship between N e and census size has to be interpreted with caution 58 . Moreover, mammalian effective population sizes are usually around 2-3 times smaller than census population sizes 59 and this ratio may be even greater in fur seals due to their strongly polygynous mating system 24 . Our results therefore suggest that the global population of Antarctic fur seals that survived sealing was considerably larger than implied by historical accounts and may have numbered in the thousands.
Several factors may play a role in explaining our results. First, standard economic theory predicts that exploitation will not necessarily lead to extinction due to the escalating costs of harvesting increasingly depleted populations 1,60 . Consistent with this explanation, historical sealing data show a geographically concordant pattern of early intense exploitation followed by greatly diminished sealing effort during the mid to late 19 th century (Fig. 2). This reduction in sealing effort coincides with a crash in the price of fur seal pelts in America and China during the 1850s, which would have substantially decreased profitability 61 . Consequently, the available historical evidence suggests that sealing probably ceased to become economically viable around the time when the main sealing grounds had been largely depleted of animals.
Having said this, historical accounts suggest that the sealers were extremely thorough in the pursuit of their quarry. Not only did their high mobility allow them to quickly discover and exploit new hunting grounds 62 , but it was also commonplace for sealing gangs to be left ashore for extended periods to opportunistically harvest any animals that remained on land 28,63 . It therefore seems unlikely that appreciable numbers of individuals could have survived ashore at the main breeding grounds. This may explain the rarity of sightings of this species during the first few decades of the 20 th century 28,[33][34][35] and the fact that some breeding populations were not re-established until as late as the 1980s [39][40][41][42][43][44] .
It is conceivable that some juveniles, subadults and non-breeders could have escaped the sealers by not coming ashore. However, if this were the case, relict populations should have survived at all of the breeding grounds, whereas our data, together with recent studies and historical accounts, suggest that that Antarctic fur seals may have been locally extirpated at Macquarie, Heard, Crozet and Marion Islands 47,50,51 . We therefore speculate that tiny remnants of once much larger breeding populations may have persisted in a handful of particularly inaccessible locations that were rarely if ever visited by sealers. Candidates for sites that may have sheltered these relict populations include the Willis Islands off the coast of South Georgia 28 , the San Telmo Islets in the South Shetlands 49 Larsøya and Kapp Norvegia at Bouvetøya 36 and any one of a multitude of small islands in the Kerguelen Archipelago 48 .
This study contributes towards wider debate about the factors shaping genetic diversity and its importance for conservation. In particular, many studies have reported links between known historical bottlenecks and low levels of genetic diversity [64][65][66][67] . However, the implied loss of diversity is often greater than predictions from population genetic theory 68,69 , suggesting that other mechanisms such as selective sweeps could be involved 70 . Additionally, the general perception that bottlenecked populations have low genetic variability may also have been influenced by publication bias 68 . In the current study, it was not possible to directly quantify the amount of genetic diversity lost as a result of sealing due to a lack of historical samples. However, coalescent simulations provided an expectation of the likely magnitude of genetic erosion given the strength of the inferred demographic declines. These simulations suggest that a recent decline from an historical N e of 10,000 to an N e bot of 200 for ten generations is unlikely to have resulted in the loss of more than around 15% of alleles. Consequently, our study supports the argument of Amos and Balmford 68 that even extreme recent bottlenecks need not necessarily lead to major losses of genetic diversity. This may further help to explain why several other pinniped species including northern, Australian and Juan Fernandez fur seals also retained high levels of genetic diversity despite having been harvested in their tens of thousands to millions [71][72][73] .
To conclude, demographic reconstruction of an entire species provided insights into the resilience of a marine mammal that was hunted in very large numbers. Our study suggests that, although local extinctions took place, Antarctic fur seals persisted in sufficient numbers in several relict populations to retain a large proportion of the species᾿ pre-sealing genetic diversity. The extreme geographical isolation of some small breeding sites in addition to diminishing economic incentives for sealing likely saved the Antarctic fur seal from extinction. The subsequent recovery of this species is a testament both to its resilience and to the success of protection measures in the 20 th century.
www.nature.com/scientificreports www.nature.com/scientificreports/ Methods tissue sampling and microsatellite genotyping. Samples were collected from a total of 2,064 Antarctic fur seal individuals from eight different breeding locations (Fig. 1, for GPS coordinates see Supplementary  Table S1). Skin samples were taken from the interdigital margin of the foreflipper with piglet ear notching pliers 74 and stored at −20°C in 20% dimethyl sulphoxide saturated with sodium chloride. DNA was extracted using a chloroform/isoamylalcohol extraction protocol 75 and genotyped at 39 microsatellite loci as described in the Supplementary methods.
Because Antarctic fur seals are known to hybridise with Subantarctic fur seals (Arctocephalus tropicalis), we additionally sampled 91 pups from a reference population of A. tropicalis at Macquarie Island in order to check the full A. gazella dataset for individuals with a high proportion of A. tropicalis ancestry (see Supplementary methods and results). We found that 49 of the Antarctic fur seal individuals were inadvertently pure A. tropicalis, while five of the A. tropicalis reference individuals were hybrids (see Figs. S11 and S12). A further 15 Antarctic fur seal individuals were also identified as hybrids. After having removed all pure A. tropicalis individuals as well as the hybrids and any additional individuals that failed to genotype at more than four loci from the dataset, we were left with a total of 2,000 Antarctic fur seal individuals, comprising 197 from the South Shetland Islands, 1,042 from South Georgia, 396 from Bouvetøya, 166 from Marion Island, 18 from Crozet Islands, 51 from Kerguelen Islands, 22 from Heard Island and 108 from Macquarie Island. population structure and genetic diversity. For the full dataset of 2,000 Antarctic fur seals, we evaluated whether population structure could be detected without prior knowledge of sampling locations using STRUCTURE version 2.3.3 76 . This program distributes individual genotypes into k clusters by subdividing the dataset in a way that maximizes Hardy-Weinberg and linkage equilibrium within the resulting clusters. The membership of each individual to a given cluster is then estimated as q, which varies between 0 and 1, with the latter indicating full cluster membership. STRUCTURE was run using ParallelStructure 77 with the following parameters: 500,000 burn-in length, 1,000,000 MCMC replications, 1-10 assumed clusters (k) and 10 iterations for each k. The most likely number of clusters was evaluated using the maximal average value of Ln P(D) as well as the Evanno method 78 . Individual admixture proportions were integrated over all simulations using pophelper 79 . For comparison, we also implemented principle component analysis (PCA) using the R-package adegenet 80 . Genetic diversity was quantified as allelic richness (A r ) within hierfstat 81 . Observed heterozygosity (H o ) was calculated with adegenet 80 , the M-ratio was calculated with StrataG 82 , the inbreeding coefficient F is was calculated with diveRsity 83 and the number of private alleles was calculated with poppr 84 . Pairwise F st values together with their 95% confidence intervals were obtained using hierfstat 81 . Bottleneck analyses. We used two coalescent-based approaches to infer recent demographic histories.
First, we quantified the proportion of loci in heterozygosity-excess using the program BOTTLENECK v. 1.2.02 85 . This analysis was based on the two-phase model with 80% single-step mutations (TPM80), consistent with recent estimates from multiple pinniped species 20 . A total of 1,000 iterations were run for each population and all other parameters were left at default values. Statistical significance was evaluated using sign, standardized differences and Wilcoxon signed ranks tests.
Second, we tested for signatures of recent population bottlenecks using a coalescent-based approximate Bayesian computation (ABC) framework. Support was evaluated for two alternative demographic models ( Supplementary Fig. S10). Briefly, to evaluate whether 18 th to 20 th century commercial exploitation caused severe demographic reductions, we defined a "bottleneck model", which incorporated a strong decline in population size within strictly bound time priors that fit the time period hunting took place. The priors for preand post-bottleneck effective population size were drawn independently from the same distribution to allow the model not only to incorporate a bottleneck, but also longer-term reductions or expansions within realistic bounds. For comparison, we defined a "non-bottleneck model" which was identical to the bottleneck model with the exception of a bottleneck.
For both models, priors for current effective populations size (N e ) and historical pre-bottleneck effective population size (N e hist) were drawn from a log-normal distribution in order to encompass a wide range of values while mainly sampling from the most likely ones. We specified N e and N e hist ~ lognorm[logmean = 10.5, logsd = 1], which concentrated around 90% of sampling in the range between 3,000 and 15,000 while also occasionally simulating much larger N e values. This prior range encompassed the estimated current census population sizes across all islands, i.e. from 3,600 at Marion Island 86 to 5 × 10 6 at South Georgia 38 . For the bottleneck model, the effective population size during the bottleneck (N e bot) was drawn from a uniform distribution between 1 and 600 (N e bot ~ U [1,600]). Start and end time for the bottleneck were defined by t bot start and t bot end and drawn from uniform distributions between ten and 40 (t bot start ~U [10,40]) and one and 20 (t bot end~ U [1,20]) generations ago respectively. This places the bottleneck in the last one to four centuries, as the generation time of the Antarctic fur seal is estimated at around nine years 87 . For the neutral model, the time parameter corresponding to the historical population size (t hist ) was drawn from a uniform distribution ranging between 10 and 40 generations ago (t hist ~U [10,40]). The microsatellite mutation rate (µ) was drawn from a uniform prior with µ ~ U[10 −5 , 10 −4 ]. A stepwise mutation model was used with the geometric parameter GSMpar reflecting the proportion of multistep mutations, uniformly distributed from GSMpar ~U[0, 0.3].
Genetic data were simulated under each demographic model for 181 individuals, reflecting the smallest genetic population (as determined by STRUCTURE) in the empirical dataset, and 39 microsatellites using fast-simcoal2 17 within the R-package strataG 82 . This resulted in a total of 2 × 10 7 datasets. We then calculated three summary statistics (as means across all loci) for the simulated data: (1) allelic richness; (2) the proportion of low frequency (<5%) alleles; and (3) the M-ratio 88 , and retained the 5,000 simulations with summary statistics closest to the empirical data using a tolerance threshold of 5 × 10 −4 .
www.nature.com/scientificreports www.nature.com/scientificreports/ To evaluate model specification and fit, we implemented leave-one-out cross-validation using the cv4postpr function of the abc R-package 89 . Here, simulations are selected at random and their summary statistics are used as pseudo-observed data. The remaining simulations are then used to classify the data into the bottleneck or the neutral model. If ABC is capable of distinguishing between the two models, a large posterior probability will be assigned to the same model under which the pseudo-observed data were generated. This procedure was repeated 100 times and the posterior probabilities for a given model were averaged to derive the rate of misclassification. Furthermore, posterior probabilities were calculated for each model using the rejection method of the postpr function and the fit of the preferred model to the empirical data was evaluated using a formal hypothesis testing procedure implemented by the gfit function of the abc R-package.
Next, we tested the accuracy of the ABC parameter estimates using leave-one-out cross validation in the cv4abc function. For a randomly selected pseudo-observed dataset, parameters were estimated using all remaining simulations via ABC using the rejection method. This procedure was repeated 1,000 times and the mean prediction error was calculated. The smaller the prediction error, the better the fit to the true parameter value. After performing the four checks, we constructed posterior distributions of the parameters using a simple rejection algorithm with the abc function of the abc package.
To provide a measure of the fit of the preferred models to the empirical data, we conducted posterior predictive checks 23,52 . After generating posterior distributions for each parameter, 1,000 multivariate parameters were sampled from their respective posterior distributions and used to simulate summary statistics a posteriori based on the preferred model. Those summary statistics were then plotted as histograms and superimposed over the observed summary statistics 23 . expected loss of genetic diversity. Finally, we used neutral coalescent simulations in fastsimcoal2 53 to estimate the expected loss of microsatellite alleles for a bottleneck effective population size (N e Bot) of 200 lasting for ten generations. As the historical effective population size N e hist is unknown but will be an important determinant of the overall level of genetic diversity, we simulated three different scenarios with N e hist values of 1,000, 10,000 and 50,000. The post-bottleneck N e is unlikely to have a major effect on diversity, as only a few generations have passed since the bottleneck, and was thus simulated equal to the N e hist within each scenario. All other model parameters were the same as in the models above. We then compared the resulting allelic richness to the allelic richness generated in equivalent models with constant effective population sizes of 1,000, 10,000 and 50,000 respectively to estimate the proportional loss of alleles due to the bottleneck. To account for uncertainty in building coalescent trees, we simulated 100,000 datasets under both the bottleneck and constant size model and treated the difference in mean allelic richness between the distributions as the average loss of alleles. To further explore the effects of varying bottleneck strength and duration, we conducted a final round of simulations based on the scenario that provided the best fit to the empirical data (N e hist = 10,000, see Results). Specifically, we quantified the expected loss of alleles for N e bot ranging 25 to 600 in steps of 25 and for bottleneck duration ranging from one to 20 generations in steps of one. As hunting ceased by the turn of the 20 th century, we fixed t bot end at ten generations ago. 1,000 simulations were conducted for each combination of parameter values.

Data availability
All of the raw data 90 are available via the Zenodo repository, doi:10.5281/zenodo.3585717.

code availability
All data wrangling steps and statistical analyses except for the heterozygosity-excess tests were implemented in R. The complete documented analysis pipeline is available via our GitHub repository: https://github.com/apaijmans/ AFS_genetic_legacy. www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/