Population genomics through time provides insights into the consequences of decline and rapid demographic recovery through head‐starting in a Galapagos giant tortoise

Abstract Population genetic theory related to the consequences of rapid population decline is well‐developed, but there are very few empirical studies where sampling was conducted before and after a known bottleneck event. Such knowledge is of particular importance for species restoration, given links between genetic diversity and the probability of long‐term persistence. To directly evaluate the relationship between current genetic diversity and past demographic events, we collected genome‐wide single nucleotide polymorphism data from prebottleneck historical (c.1906) and postbottleneck contemporary (c.2014) samples of Pinzón giant tortoises (Chelonoidis duncanensis; n = 25 and 149 individuals, respectively) endemic to a single island in the Galapagos. Pinzón giant tortoises had a historically large population size that was reduced to just 150–200 individuals in the mid 20th century. Since then, Pinzón's tortoise population has recovered through an ex situ head‐start programme in which eggs or pre‐emergent individuals were collected from natural nests on the island, reared ex situ in captivity until they were 4–5 years old and subsequently repatriated. We found that the extent and distribution of genetic variation in the historical and contemporary samples were very similar, with the latter group not exhibiting the characteristic genetic patterns of recent population decline. No population structure was detected either spatially or temporally. We estimated an effective population size (N e) of 58 (95% CI = 50–69) for the postbottleneck population; no prebottleneck N e point estimate was attainable (95% CI = 39–infinity) likely due to the sample size being lower than the true N e. Overall, the historical sample provided a valuable benchmark for evaluating the head‐start captive breeding programme, revealing high retention of genetic variation and no skew in representation despite the documented bottleneck event. Moreover, this work demonstrates the effectiveness of head‐starting in rescuing the Pinzón giant tortoise from almost certain extinction.

postbottleneck contemporary (c.2014) samples of Pinzón giant tortoises (Chelonoidis duncanensis; n = 25 and 149 individuals, respectively) endemic to a single island in the Galapagos. Pinzón giant tortoises had a historically large population size that was reduced to just 150-200 individuals in the mid 20th century. Since then, Pinzón's tortoise population has recovered through an ex situ head-start programme in which eggs or pre-emergent individuals were collected from natural nests on the island, reared ex situ in captivity until they were 4-5 years old and subsequently repatriated.
We found that the extent and distribution of genetic variation in the historical and contemporary samples were very similar, with the latter group not exhibiting the characteristic genetic patterns of recent population decline. No population structure was detected either spatially or temporally. We estimated an effective population size (N e ) of 58 (95% CI = 50-69) for the postbottleneck population; no prebottleneck N e point estimate was attainable (95% CI = 39-infinity) likely due to the sample size being lower than the true N e . Overall, the historical sample provided a valuable benchmark for evaluating the head-start captive breeding programme, revealing high retention of genetic variation and no skew in representation despite the documented bottleneck event. Moreover, this work demonstrates the effectiveness of headstarting in rescuing the Pinzón giant tortoise from almost certain extinction.
The Pinzón giant tortoise, Chelonoidis duncanensis (previously Chelonoidis ephippium; Turtle Taxonomy Working Group, 2017), is endemic to Pinzón Island (18 km 2 area) in the Galapagos ( Figure 1) and consists of a single population. Historically, Pinzón giant tortoises numbered in the thousands, but mass harvesting for food by humans in the early to mid-1800s dramatically reduced the population size (Townsend, 1931 were introduced to the island and depredated all hatchling tortoises, resulting in no young Pinzón giant tortoises surviving from that time onward (Pritchard, 1996). Surveys conducted in the 1960s located 100 mature C. duncanensis on the island and estimated the census population size to be 150-200 individuals, all born before the introduction of rats (MacFarland, Villa, & Toro, 1974). At that point, the Pinzón giant tortoise was a species of "living dead" due to the absence of any recruitment, seemingly destined to become extinct when the last of the remaining ageing adults died.
Faced with the potential loss of the Pinzón giant tortoise, the Charles Darwin Research Station initiated a head-start programme in 1965, collaboratively managed with the Galapagos National Park Directorate. Eggs or pre-emergent individuals were collected from natural nests on Pinzón Island and reared in captivity until they were 4-5 years old. Juveniles at this age can avoid predation by rats and were repatriated (Cayot, 2008). The programme successfully repatriated more than 800 juvenile tortoises over the past 50 years.
However, unless the cause of hatchling mortality in the wild was addressed, the programme would have to operate in perpetuity (Cayot, 2008). In recognition of this, the Galapagos National Park carried out a rat poisoning campaign in December 2012. In 2014, the first wild-born hatchlings of C. duncanensis were observed in over a century, and the rat eradication project was considered a success (Tapia, Malaga, & Gibbs, 2015). Taken together, these efforts have helped rescue the Pinzón giant tortoise from extinction.
The known history of decline and recovery of C. duncanensis provides a rare opportunity to perform a direct evaluation of the relationship between current genetic diversity and past demographic events, including the impact of the head-start programme on prospects for long-term survival of the species. Importantly, a large number of museum specimens (n = 86 adults, Figure 1) were collected from Pinzón Island during a 1905-1906 California Academy of Sciences expedition to Galapagos (Van Denburgh, 1914). Given the long lifespan (>100 years) and generation time of Galapagos giant tortoises (25 years, Throp, 1975), coupled with the well-documented history of decline in the Pinzón species, these museum specimens are likely representative of the prebottleneck population. Thus, they provide a valuable reference point for investigating changes in both the extent and distribution of genetic variation over the past century. Notably, despite the historical specimens being collected more than 100 years ago, only a single generation likely elapsed since their collection due to the "recruitment pause" between the late 1890s and the 1960s. As the historical specimens were collected from a naturally reproducing population, the distribution of genetic variation in that sample can be used to evaluate whether human-assisted survival of individuals through the head-start programme has caused a skew in genetic variation.
In this study, we assessed the genetic impacts of a known bottleneck, using pre-and postbottleneck samples of Pinzón giant tortoises-a species that has now demographically recovered. We paired restriction-site-associated DNA sequencing (RAD-seq) with targeted capture techniques to collect genome-wide single nucleotide polymorphism (SNP) genotypic data from historical (prebottleneck) and contemporary (postbottleneck) samples. We then used this SNP data set to directly investigate changes in genome-wide diversity over time. Additionally, given that the majority of the tortoises in the contemporary sample are the product of the head-start programme, we evaluated the degree to which genetic variation has been impacted by this conservation intervention. recorded. Blood (0.1-1.0 ml) was collected from the brachial artery.

| MATERIAL S AND ME
Additional details regarding sample collection and DNA extraction can be found in the Supporting Information Supplemental Methods in Appendix S1. Field notes suggest that specimens were collected from throughout the island. Femurs attached to carapaces were sampled from 78 of these specimens accessioned at the California Academy of Sciences (Supporting Information Data S1). All individuals were adults, 57 females and 21 males, with carapace lengths ranging from 53 to 87 cm (Van Denburgh, 1914). DNA was extracted from wedge cuts of bone in a dedicated ancient DNA laboratory at The University of British Columbia Okanagan using a modified version of extraction protocol Y from Gamba et al. (2016)

described in the Supporting Information
Supplemental Methods in Appendix S1.

| Molecular and bioinformatic methods
The detailed information regarding data collection and bioinformatic processing can be found in Supporting  The reference "genome" used for alignment of nuclear capture loci consisted of 140-bp-long target sequences, as well as 100 bp of flanking sequence on either end obtained from a draft genome of Chelonoidis abingdonii (Quesada et al., in preparation). Sequences were processed using the BAM pipeline in PALEOMIX (version 1.2.6, Schubert et al., 2014), which employs other standard bioinformatic tools alongside native scripts to support the pipeline. To allow the sequences to be compared between historical captures and contemporary RAD-seq data, the fastq files retained following the clone_ filter step in the STACKS workflow were run through PALEOMIX, using the same procedure as for the historical individuals, excluding the DNA damage correction, starting at the mapping stage.
Genotype calling was performed on the combined BAM files generated from the historical and contemporary individuals and from the contemporary individuals alone to produce two data sets: SNPs genotyped in both temporal samples, and a larger pool of SNPs genotyped in the contemporary sample only. We tested several different combinations of mapping and genotype calling approaches on the contemporary individuals to evaluate the impact of different workflows on the outcomes. The details and results are provided in the Supporting Information Supplemental Methods in Appendix S1.
Here, we present the analyses used with the data set from the best combination. We implemented stringent filtering of the variable sites and ultimately retained only the first SNP per target region to produce a data set of loci for population genetic analyses. Only individuals with >50% of SNPs genotyped were retained.
We assessed error associated with genotyping of SNP loci by calculating the number of genotype mismatches between two pairs of replicate contemporary individuals (tortoise ID #'s A025 and G154, Supporting Information Data S1) that had been processed independently from DNA extraction onwards.

| Within-population diversity
We used the genotypic data to calculate standard measures of within-population genetic diversity for the contemporary and historical samples separately, including heterozygosity and G IS using GENODIVE V 2.0b27 (Meirmans & Van Tienderen, 2004). Individual inbreeding coefficients were calculated using VCFTOOLS (Danecek et al., 2011). Pairwise relatedness (Queller & Goodnight, 1989) was calculated within the contemporary and historical samples separately using the Related package (Pew, Muir, Wang, & Frasier, 2015) in the R statistical package, version 3.2.2 (R Development Core Team, 2015). A genotyping error rate of 4.5% was used (empirically determined, see Supporting Information Table 1 in Appendix S1).
Heterozygosity, inbreeding coefficients and pairwise relatedness were calculated using the small (n = 2,218 SNPs, including loci in common among temporal samples) and large (n = 7,730 SNPs, including only loci in the contemporary sample) data sets to assess the impact of the number of loci on the results.
Effective population sizes (N e ) of the contemporary and historical sample groups were calculated using the bias-corrected measure of linkage disequilibrium (Hill, 1981;Waples, 2006;Waples & Do, 2010), as implemented in NeESTIMATOR V2.1 (Do et al., 2014). We explored the effects of the number of individuals on estimates of N e by creating 50 random subsets in increments of 10 individuals, from 10 to 100 individuals, pulled from the contemporary sample group and using both the 2,218 and 7,730 SNP data sets. We estimated N e for each subset using a minor allele frequency (MAF) cut-off of 0.05 and viewed the results in R.

| Population substructure analyses
The evidence for substructure within the combined sample (contemporary plus historical individuals) was assessed using Bayesian clustering analysis, as implemented in STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). Run length was set to 500,000 Markov chain Monte Carlo replicates after a burn-in period of 100,000 using correlated allele frequencies under an admixture model with alpha set to 0.5. We varied the number of clusters (K) from one to four, with ten iterations of each. The most likely number of clusters was determined by plotting the log probability of the data (ln Pr(X|K)) across the range of K values tested and selecting the K where the value of ln Pr(X|K) plateaued, as suggested in the STRUCTURE manual.
We also used the model-free discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010) implemented in adegenet (Jombart, 2008) in R. The best-fit value of K was selected using the find.clusters function and Bayesian information criterion (BIC). The chosen value of K was based on the minimum number of clusters after which the BIC decreased by a negligible amount.
To assess whether there were differences in the genetic diversity captured in the head-start programme over time, we used carapace size as a proxy for age of sampled individuals in order to divide samples from the contemporary sample group into two age classes and repeated the calculation of diversity metrics (see Within-population diversity, above). The age classes were "adult" samples, those with a curved carapace length >65 cm (n = 82), and "young" samples, those with carapace lengths <35 cm and >15 cm (n = 29). Individuals <15 cm were excluded because they are wildborn hatchlings that were not part of the head-start programme.
Admittedly, these age classes are somewhat arbitrary, but were chosen to represent nonoverlapping segments of the contemporary sample group that was head-started either in the early years of the programme ("adults") or very recently ("young"). This grouping allowed comparisons of levels of diversity within each time point and of the patterns of pairwise relatedness within and among the two temporal sample groups.
The Phi ST metric of differentiation was calculated between temporal sample groups in GENODIVE with significance assessed using 999 permutations. Exact tests for differences in allele frequencies between the historical and contemporary sample groups, and between the historical sample and two age classes within the contemporary sample group, were performed in GENEPOP v4.5 (Raymond & Rousset, 1995), with significance assessed using an adjusted pvalue based on the correction for false discovery rate described by Benjamini and Yekutieli (2001).
To evaluate evidence for spatial genetic structure within the contemporary sample group, we compared the straight line geographic distance between pairs of individuals at the time of sampling, calculated using the GEOGRAPHIC DISTANCE MATRIX GENERATOR (Ersts, 2012), and their pairwise relatedness. Only individuals with a curved carapace length >50 cm (i.e., mature adults, n = 99) were used in this analysis to exclude recently repatriated cohorts that have not yet had time to disperse away from the release sites ( Figure 1b). No spatial data were available for the museum specimens, precluding a similar analysis for the historical sample group.

| RE SULTS
We  Table 2 in Appendix S1, see Supporting Information Supplementary Results in Appendix S1 for genotype quality results). Levels of heterozygosity and inbreeding were very similar in the historical and contemporary sample groups, as was the mean relatedness among individuals (  Table 3 in Appendix S1, see Supporting Information Supplementary Results in Appendix S1). Thus, the 2,218 SNP data set that overlapped with the historical sample group was used for all downstream analyses.
We found no geographical clustering of closely related contemporary adult individuals (R 2 = 0.0006, Supporting Information Figure 1 in Appendix S1). Population structure tests (i.e., STRUCTURE and DAPC) did not split the contemporary and historical samples into separate clusters (Supporting Information Figure 2 in Appendix S1).
Similarly, levels of differentiation were low (but still statistically significant) between temporal samples as measured by Phi ST , and only a small proportion of loci exhibited significantly different allele frequencies (Table 2).
Adult and young subsets of the contemporary sample were selected for comparison as they are at the two ends of the age continuum, representing either the early periods of the head-start programme ("adults") or more recent head-start cohorts ("young").
Heterozygosity was slightly higher, and measures of inbreeding, G IS and F, were more negative in the young than the adult group (  Queller and Goodnight (1989) relatedness; SNP, single nucleotide polymorphism.
Measures are for the historical and contemporary samples and subset "age" classes within the two contemporary samples. The analyses were carried out on two SNP data sets: a data set including the 2,218 SNPs common to both temporal samples, and a 7,730 SNP data set from the contemporary sample only. a "Adult" refers to the subset of individuals in the contemporary sample group with a curved carapace length >65 cm; "Young" refers to the subset of individuals in the contemporary sample with a curved carapace length <35 and >15 cm. distributions ( Figure 2b). All Phi ST values were significant, but very low (Table 2). There was slightly greater differentiation between the young contemporary samples and the historical samples than between the adult contemporary samples and the historical samples (Table 2). These results were not influenced by the inclusion of three highly related young individuals (Supporting Information Table 4 in Appendix S1).  (Cornuet & Luikart, 1996) and the fact that the size of our data set precluded the application of site frequency spectrum-based F I G U R E 2 Frequency distributions of pairwise relatedness. Relatedness values are as follows: (a) among individuals within the contemporary and historical sample groups; and (b) among the two size classes "adults" and "young" from the contemporary sample group, and the relatedness values between pairs of "adult" and "young" contemporary individuals. Relatedness estimates were calculated following Queller and Goodnight (1989) based on 2,218 SNP loci. SNP, single nucleotide polymorphism F I G U R E 3 The effective population size estimates from subsets of different sample sizes taken from the pool of contemporary Pinzón giant tortoises, using a minor allele frequency of 0.05. The left and right plots are separated to allow for different scales on the y-axes. For each sampling size, there were 50 subsets drawn based on the 2,218 (blue) and 7,730 (black) SNP data sets. The points are jittered, and the bars indicate 95% jackknife confidence intervals. Point estimates that were negative or infinity, or had confidence intervals that included infinity are not presented (29% of the runs, exclusively in the subsets with sample sizes of n = 10-30 individuals). Thus, for example, for a sample size of 10, only one estimate is presented based on the 7,730 SNPs, and no values are presented for the 2,218 SNP data sets. SNP, single nucleotide polymorphism  (Beheregaray, Ciofi, Geist et al., 2003). In that case, a major volcanic eruption approximately 74,000-120,000 years ago was inferred as the likely cause for this demographic change (Beheregaray, Ciofi, Geist et al., 2003).

| Genomic consequences of population decline and rapid demographic recovery
Based on the results of the present study, it seems likely that the bottleneck affecting Pinzón giant tortoises was not as dramatic both in terms of magnitude of size reduction and duration as the event experienced by the Volcano Alcedo tortoises.
Population genetic theory predicts a decrease in N e due to a bottleneck (reviewed in Charlesworth, 2009;Wright, 1940). We had intended to compare N e between the pre-and postbottlenecked samples, but the small sample size for the historical group (25 successfully genotyped of 78) precluded accurately estimating N e (Table 1). Our simulations that explored the impact of sample size via downsampling of the contemporary group indicated that 50 or more individuals are required for accurate estimates of N e (Figure 3). This conclusion follows those of England, Cornuet, Berthier, Tallmon, and Luikart (2006), who found that sample sizes equal to or greater than the true N e are required for the linkage disequilibrium method to produce reliable estimates. Our estimated N e in the contemporary sample was roughly the same when calculated based on the 2,218 or 7,730 SNP data sets (N e = 58 and 59, respectively, Table 1). This value is higher than a previous estimate for Pinzón giant tortoises of N e = 26 (95% CI = 17, 45) based on microsatellite genotypic data for 24 individuals (Garrick et al., 2015). Interestingly, conversion of the SNP-based N e translates to an estimated census population size (N c ) of 536 based on the Frankham (1995) ratio (N e = 0.11 N c ), which more closely approximates the number of different individuals that were encountered (n = 420) during a population survey conducted concurrently with our sample collection in 2014.
One persistent pattern found in all previous studies of Pinzón giant tortoises using microsatellites has been heterozygote deficit relative to Hardy-Weinberg expectations and significantly positive inbreeding coefficients Garrick et al., 2015;Jensen et al., 2015). Similar results were found in a recent study employing double-digest RAD (ddRAD) sequencing to quantify genetic diversity and reconstruct population structure across all 12 extant species of Galapagos giant tortoises .
At >26,000 SNPs, Miller et al. (2018) found significant inbreeding within the 10 Pinzón individuals included in that study. In contrast, we detected heterozygote excess in both the historical and contemporary samples, and the two measures of inbreeding we calculated are slightly negative and, in the case of G IS , not significantly different from zero (Table 1). These patterns could be due to assembly parameters allowing promiscuous mapping of reads, but we rigorously tested our assembly approach and the same pattern of slightly negative inbreeding coefficients was repeatedly found (see Supporting Information Table 1 in Appendix S1). At least one other study has found that the methodological differences between ddRAD-seq and single-enzyme RAD-seq can produce different heterozygosity patterns in the same samples (Flanagan & Jones, 2018). Although the extent to which technical artefacts associated with polymorphisms in restriction enzyme cut sites may bias estimation of population genetic parameters is a matter of some debate (see Andrews et al., 2014;Puritz et al., 2014), moving forward, genome-wide sequencing, even at low coverage (e.g., Therkildsen & Palumbi, 2017), should ameliorate such concerns.

| Genetic legacy of the head-start programme
The demographic recovery of the Pinzón giant tortoise species was achieved via a head-start programme. However, this action carried the possibility of skewing genetic contributions to subsequent generations due to unequal representation of the surviving individuals.
Even if all 150-200 adults on the island immediately following the bottleneck (MacFarland et al., 1974) contributed offspring to the head-start programme, genetic diversity could have become skewed due to the over-representation of certain families in the head-start generations. This scenario may have occurred if some mate pairs naturally produced more offspring than others, and/or if eggs or pre-emergent individuals were collected nonrandomly. A previous microsatellite study evaluated the genetic representativeness of the cohorts of the head-start programme hatched in 2007 and 2009, relative to a sample of 72 adults from the wild population (Jensen et al., 2015). The conclusions were that the cohorts were genetically diverse, but that there was genetic variation in the wild adults not represented in the head-start cohorts. Surprisingly, there were a number of alleles present in the hatchlings that were not found in the sample of adults, indicating that despite constituting a large proportion of the reproducing population, that adult sample was not fully representative of the breadth of genetic variation in the population (Jensen et al., 2015).
TA B L E 2 Measures of genetic differentiation calculated between the historical and contemporary sample groups, and between "age" class subsets within the contemporary sample group Notes. "Adult" refers to the subset of individuals in the contemporary sample group with a curved carapace length >65 cm; "Young" refers to the subset of individuals in the contemporary sample with a curved carapace length >15 and <35 cm. *Denoting significance at p < 0.001 for the Phi ST , DAF, the proportion of loci with significantly different allele frequencies (adjusted p-value 0.0059).
In the current study, we were able to more fully assess how successful the head-start programme has been in maintaining the extent and distribution of genetic diversity by calculating pairwise relatedness among individuals in the contemporary and historical samples and comparing the distributions of relatedness values. As individuals in the historical sample were naturally reproducing, they serve as a benchmark for the typical distribution of relatedness prior to the bottleneck and establishment of the head-start programme.
We found that the distribution of pairwise relatedness in the contemporary sample matched the distribution for the historical sample (Figure 2a), indicating that overall, the head-start programme collected eggs/individuals in a way that was not biased towards certain families. We further compared the relatedness distributions of the adult and young subsets within the contemporary sample to see whether this finding applies to both the early and recent periods of the head-start programme and not just overall. We found largely overlapping distributions of relatedness within the two age subsets, with a slight excess of higher relatedness in the subset of young individuals ( Figure 2b). Although difficult to test directly, this increase in relatedness in the young cohort may be due to an unusually large head-start cohort collected in 2009 and reintroduced 4 years later (Jensen et al., 2015). There were some indications that the young group is more differentiated than the adults from the historical group, as there were a larger proportion of loci with significantly different allele frequencies and a slightly higher Phi ST value (Table 2).
These patterns are likely due to background levels of genetic drift operating in a small population over time. More broadly, levels of heterozygosity and inbreeding were similar between the two age class subsets (Table 1), suggesting that genetic diversity appears to have been captured consistently by the head-start programme over a ~50 year period.

| Insights from temporal sampling
In this study, we had the benefit of historical samples to provide context for interpreting levels of diversity observed in the contemporary population. Prebottleneck sampling of multiple individuals is unavailable for most species of conservation concern, and so indirect estimates of the severity and genetic impacts of bottlenecks must be relied upon. In some cases, the bottlenecked population can be compared to a stable population of the same species (Whitehouse & Harley, 2001) or related species (Akst, Boersma, & Fleischer, 2002;Waldick, Kraus, Brown, & White, 2002) to indirectly assess the genetic impacts of population decline. The previous archipelago-wide studies of Galapagos giant tortoises have taken this approach and compared levels of variation in each species to gain insights into population history (e.g., Beheregaray, Ciofi, Geist et al., 2003;Ciofi, Milinkovitch, Gibbs, Caccone, & Powell, 2002;Garrick et al., 2015). However, interpreting the baseline provided by comparing Pinzón giant tortoises to other Galapagos tortoises has been complicated as C. duncanensis maintains higher levels of genetic variation, particularly in the mitochondrial genome, than most of the other extant species, despite having gone through a substantial bottleneck (MacFarland et al., 1974;Pritchard, 1996). In a parallel study of mitochondrial genetic variation in the temporal samples of Pinzón giant tortoises, Jensen et al. (2018) found haplotypic diversity was equal over the mitogenome, although insights from particular genes or subsets of genes were inconsistent. Here, through direct comparison of pre-and postbottleneck samples, it is apparent that the extent and distribution of nuclear genetic variation have been maintained through time (Table 1, Figure 2).
Although the study design used here provides unique insights, methods for collecting and analysing population genomic data from historical specimens, in the context of a temporal study, are still maturing. In Appendix 1, we detail important considerations for study design and implementation, including those related to initial SNP discovery, trade-offs associated with balancing sample size and SNP number, and quality control.

| CON CLUS IONS
Over the past several decades, considerable attention has been given to documenting the genetic consequences of population declines (Frankham, 2005;Frankham et al., 2017), yet rarely have there been opportunities to test theoretical population genetic predictions using temporal pre-and postbottleneck sampling or to evaluate the impacts of conservation programmes using samples from before the intervention. In the case of Pinzón giant tortoises, the harvesting of many individuals in the early 20th century for museum collections has provided an opportunity to directly assess genetic patterns associated with population decline and recovery, in this case facilitated through a head-start programme. Given the increased capacity to mine the genome even from highly degraded sources of DNA, empirical studies using temporally spaced samples will continue to enrich our understanding of evolutionary processes and help inform conservation action.

ACK N OWLED G EM ENTS
We acknowledge the skill and commitment of members of the

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
The data sets generated in this manuscript are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.7tp3sg0.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section at the end of the article.

Insights from temporal sampling
As methods for collecting and analysing population genomic data from historical specimens in the context of a temporal study are still maturing, we provide here some important considerations for study design and implementation. We chose to use restrictionsite-associated DNA sequencing (RAD-seq) for the contemporary sample to discover variable regions of the genome that could be targeted using hybridization baits in the historical samples. This procedure may have introduced some ascertainment bias, in that only RAD loci known to have variable sites in the contemporary samples were targeted. Our filtering of variable sites was designed to minimize this potential bias by selecting the first occurring single nucleotide polymorphism (SNP) in the RAD-seq read when both contemporary and historical individuals were considered, which was not necessarily the SNP identified in the initial analysis.
In our case, there were 48 SNPs retained in the final data set that were not polymorphic in the contemporary sample group, suggesting that SNPs other than those identified during the initial analysis were ultimately retained. Additionally, in the data set used for temporal analyses, we applied the minor allele frequency filter to the combined sample, to ensure that the loci retained were informative in both temporal groups.
During the filtering process for the SNP data set, a trade-off was necessary between the number of historical individuals retained and the number of SNP loci that met filtering criteria. Here, the final data set consisting of 25 historical individuals genotyped at 2,218 SNP loci is lower than the data set initially targeted (baits were constructed for 8,918 SNPs), but is still a very encouraging outcome that demonstrates the possibilities afforded by targeted capture.
Furthermore, the larger number of loci genotyped in the contemporary samples (7,730 SNPs) provided almost exactly the same diversity estimates as the subset of 2,218 loci (Table 1), lending confidence that our findings would remain unchanged with additional loci from the historical samples. The mean depths for the historical and contemporary SNP data sets were about equal, with genotyping error rates similar to what would be expected based on the coverage (Fountain, Pauli, Reid, Palsbøll, & Peery, 2016).
We were able to directly assess genotyping error rates using two pairs of contemporary samples that had been run in duplicate.
Duplicate historical individuals were not included in the study design, so a similar calculation specific to that data set is not possible.
Errors in genotypes can arise during data collection steps, including PCR errors during amplification, sequencing errors and postmortem DNA degradation in the case of historical samples; or during data processing, including alignment errors and filtering steps (e.g., read depths contributing to genotype calling). We expect the levels of genotyping error due to sequencing errors and data processing to be similar between the historical and contemporary data sets, as the sequencing was run using the same chemistry and instrument, and using the same data quality filters. Postmortem DNA degradation was accounted for during the data processing using MAPDAMAGE (Jonsson, Ginolhac, Schubert, Johnson, & Orlando, 2013) to rescale nucleotide quality scores suspected to be impacted by DNA degradation. Thus, in our data set, we expect the genotyping error rates calculated from the contemporary samples to be reasonably representative of error rates in the historical samples. This may be especially the case given that the historical individuals had on average greater read depth for genotypes (18.3×) than contemporary individuals (14.0×). Directly assessing genotyping error rates using duplicate samples would be preferable, however, and should be factored into the design of future studies.