## Abstract

Genetic diversity provides the raw material for populations to respond to changing environmental conditions. The evolution of diversity within populations is based on the accumulation of mutations and their retention or loss through selection and genetic drift, while migration can also introduce new variation. However, the extent to which population growth and sustained large population size can lead to rapid and significant increases in diversity has not been widely investigated. Here, we assess this empirically by applying approximate Bayesian computation to a novel ancient DNA dataset that spans the life of a southern elephant seal (*Mirounga leonina*) population, from initial founding approximately 7000 years ago to eventual extinction within the past millennium. We find that rapid population growth and sustained large population size can explain substantial increases in population genetic diversity over a period of several hundred generations, subsequently lost when the population went to extinction. Results suggest that the impact of diversity introduced through migration was relatively minor. We thus demonstrate, by examining genetic diversity across the life of a population, that environmental change could generate the raw material for adaptive evolution over a very short evolutionary time scale through rapid establishment of a large, stable population.

## 1. Introduction

Understanding how micro-evolutionary processes lead to changes in evolutionary potential in the form of genetic diversity is a fundamental goal of population genetics. Population genetic diversity can be gained either *in situ* through mutation, through immigration or through some combination of the two [1,2]. Variation in the form of new mutations is expected to enter a population at a rate of 2*N*_{e}*u*, where *N*_{e} is the effective population size and *u* the neutral mutation rate per generation, while genetic drift decays variation at the rate of 1/2*N*_{e} [3]. New mutations within a population can either be fixed or (more commonly) lost [4,5]. Fixation time of new unlinked neutral mutations within a randomly mating population is strongly reliant on *N*_{e}, taking some 4*N*_{e} generations when *N*_{e} is constant [4,6]. However, change in population size, a fundamental characteristic of many wild populations, can significantly impact on this time frame [4,6]. Changes in genetic diversity can also result from the movement of individuals (migration) and their genes (gene flow) into or out of a population. Thus, rapid changes in population genetic diversity may result from demographic processes (e.g. immigration) or *in situ* processes (e.g. mutation and genetic drift). These processes will also be influenced by changes in *N*_{e}, which may result from varying physical conditions, resource availability, habitat availability, predator density, human impacts [7] or skewed effective breeding ratios.

While theoretical models predict changes in levels of genetic diversity within a population over time, it is uncommon to be able to track such changes in vertebrate populations empirically, particularly over time frames longer than those of pedigree studies. However, molecular analyses of ancient DNA allow these processes to be assessed more fully, especially when ancient DNA samples span the life of the population. Here, we investigate such a dataset, provided by the founding and subsequent extinction of a colony of southern elephant seals. This population existed along the Victoria Land Coast (VLC), Antarctica for much of the Holocene, resulting from the expansion of suitable habitat owing to altered environmental conditions [8,9]. Our study is inspired by the observation that this short-lived population had apparently accumulated extensive novel variation, most of which was lost when the population went to extinction [9], but the relative contribution of migration and *in situ* mutation to this increase in genetic diversity remains unknown.

Radiocarbon dates from 223 individual seal remains provided upper- and lower-bound age estimates on the population of 7200–270 years before present (YBP) [8,9]. Subsequent analyses using a southern elephant seal-specific substitution rate estimated from these dated samples supported the notion of an accelerated short-term time-dependent rate [10], and agreed with estimates obtained from various other aDNA datasets (e.g. [11,12]), and the mean human mitochondrial control region pedigree rate estimate [13]. Importantly, population genetic analyses of VLC founder dynamics conducted using this substitution rate fitted well the time frame recovered from geological dating of Holocene raised beaches along the VLC, and the fate of the Ross Sea Ice Sheet, which indicated that this coastline would have been inhospitable to southern elephant seals prior to around 8000 YBP [8,14]. Elephant seals require breeding habitat with access to open water, and the VLC coastline was beneath the Ross Sea Ice Sheet prior to this time [8]. The retreat of the ice shelf was followed by a period of relative warmth when the seal colony apparently thrived and expanded, before more recent cooling around 1000 YBP led to an increase in land-fast sea-ice—a likely cause of the demise and eventual extinction of this population [8,9].

Molecular analyses of these ancient DNA samples from the VLC and from contemporary samples from seven putative populations, representing all four major extant southern elephant seal breeding stocks [15], have provided considerable insight into the origin and demographic history of the VLC seals [9]. Our earlier research demonstrated that the VLC seals were an independent breeding colony—most likely founded by seals from Macquarie Island (MQ)—which grew in size through the Mid-Holocene, before declining drastically around 1000 YBP, resulting in eventual extinction and the possible return of some VLC animals to the source population at MQ [9]. Inferred timings of these demographic events derived from molecular data are consistent with inferences based on fossil radiocarbon dates and geological data. However, it was not clear whether this scenario was sufficient to explain the very high levels of diversity of the VLC population or the relative importance of migration compared with *in situ* mutation.

We therefore applied approximate Bayesian computation (ABC) methods to the VLC and MQ dataset to test hypotheses about the timing and relative importance of different processes, thereby testing hypotheses about life history and behaviour (such as the likelihood of the VLC being a breeding as opposed to a moulting population, and the degree of insularity of the newly founded population). ABC enables the joint assessment of a large number of increasingly complex demographic models and allowed us to examine in detail processes leading to significant changes in genetic diversity, from initial founding to the eventual extinction of a population.

## 2. Material and methods

### (a) Data

The VLC and MQ dataset analysed here comprised 223 ancient VLC samples, ranging in age from ±7200 to 270 YBP, and 49 contemporary samples from MQ (see electronic supplementary material, table S1; GenBank accession numbers: FJ168073–FJ168343). Details of laboratory DNA methodology (including extractions, PCR, sequencing of a 325 bp fragment of the HVR1 section of the mitochondrial control region, cloning, ‘independent replication’ and other strict aDNA protocols), as well as calculation of aDNA calendar dates from radiocarbon dating that incorporated a time-dependent Southern Ocean marine reservoir effect, are presented elsewhere [9].

### (b) Coalescent models and summary statistics

Based on earlier population genetic results, including isolation-with-migration analyses [9], a series of increasingly complex models were examined within an ABC framework. These ranged from a simplified single population (model 1) or model with two populations that diverged 12 000 YBP (model 2), to those incorporating exponential population growth and bi-directional migration (model 7; figure 1). Models of intermediate complexity examined colonization from VLC to MQ (model 3, the opposite of what we believe occurred), colonization from MQ to VLC (model 4), ongoing migration from MQ to VLC after colonization (model 5) and a single migration event from VLC to MQ when the VLC population was extirpated (model 6). We also added alternative versions of models 4–7 that included exponential growth (consistent with known demographic potential) in VLC for the first 500 years after colonization, followed by a stable population (models 4b–7b). For models with colonization, we specified that the timing of colonization was 8000 YBP, consistent with the retreat of ice from VLC and the first available breeding habitat. This was a conservative choice given that a more recent colonization would require more rapid population growth.

We used serial coalescent theory, as implemented in Bayesian Serial SimCoal [16], to simulate 1 million datasets from each of our 11 models. We set the mutation rate to 9.8 × 10^{−7} mutations site^{−1} yr^{−1}, as estimated from dated aDNA samples, with a transition/transversion ratio of 0.8564 and a mutation rate gamma distribution shape parameter of 0.2003 [9]. The sample sizes, sample dates and length of the locus (325 bp) matched our collected data (see electronic supplementary material, table S1). We used uniform priors on population size and migration rates (table 1; electronic supplementary material S2) so as to specify little prior knowledge about these quantities. Uniform priors for population size and migration rates have been commonly recommended or used in ABC [17–19]. The prior on back-migration in models 6 and 7 was also uniform 0–100% (see electronic supplementary material, table S2). For both on-going migration from MQ to VLC and the back-migration event from VLC to MQ, 100% migration corresponded to full replacement of the sink population with genes from the source population, while 50% migration corresponded to replacement of only half the sink population. In addition, we evaluated the sensitivity of our conclusions to alternative priors by generating 333 333 simulations from each of the 11 models with (i) log uniform migration rates from *e*^{−9} to *e*^{0}, or (ii) wider priors on the size of the colonizing population (uniform 0–10 000). Colonizing populations up to 10 000 in a single generation are highly unlikely given that MQ only has 10 000 females.

From the real data and from the simulations, we then calculated a series of summary statistics on each of three statistical groups: MQ (*n* = 49), VLC 3000–7100 YBP (*n* = 64) and VLC 0–3000 YBP (*n* = 159). Statistical groups were chosen to examine temporal dynamics in VLC while maintaining a sufficiently large sample size for the ancient group. However, we also conducted a sensitivity analysis after adding an additional statistical group for the earliest VLC years (6000–7100 YBP) and trimming the other ancient VLC statistical group to 3000–6000 YBP. The sensitivity analysis was conducted with 333 333 simulations from each of our 11 models.

We used the number of haplotypes (*H*), number of segregating sites (*S*) and average number of pairwise differences within each of the statistical groups (*π*) as summary statistics (table 2; electronic supplementary material S3). We also used *F*_{ST} values calculated between MQ and each of the two VLC statistical groups. We chose these summary statistics because they reflect the demographic changes in which we were interested [17], namely changes in population size (*H*, *S* and *π* at two time points) [20,21] and the exchange of migrants through time between MQ and VLC (*F*_{ST}s) [22] (see electronic supplementary material, table S3). Preliminary simulations also suggested that they were sufficient to distinguish among our competing hypotheses.

### (c) Approximate bayesian computation

After simulating the models and calculating the summary statistics, we conducted our ABC analysis in three steps: model selection, parameter estimation and quality control. For model selection, we used a weighted multinomial logistic regression to calculate the relative support for each model [23] from the 5% of all 11 million simulations for which *δ* (the Euclidean distance between the observed and simulated summary statistics) was smallest. Parameter estimation followed a similar process. From the 0.5% of simulations from each model for which *δ* was smallest, we used weighted local linear regression to estimate the posterior distributions of each parameter [24]. Before the parameter calculations, we log-transformed MQ and final VLC population sizes, and logit-transformed migration rates and the size of the VLC colonizing population. Both transformations were applied to facilitate parameter estimation and projection back onto untransformed axes. We report the median, mode and 95% highest probability density (HPD) for each parameter. Calculations were done in R v. 2.15.3 with the abc package [25].

For quality control, we followed the recommendations of Bertorelle *et al.* [17] and generated pseudo-observed datasets (PODs). We sampled 1000 parameter sets from the posterior distributions of each model and used these to simulate PODs. By sampling from the posteriors, we account for our full knowledge and uncertainty about the true model parameters. Previous authors have used only the best parameter estimates [17] or the full prior distributions [26]. We re-ran our model selection calculation on each of 100 PODs from each model (1100 total PODs) to evaluate our model selection procedure when the true model was known. This process allowed us to calculate model choice error rates. In addition, these results also allowed us to calculate a *corrected* relative support for each model following the procedure of Fagundes *et al.* [26]:
for each potential model *X*, the originally identified best model A, and the originally observed relative support *a* for model A. The procedure calculates the probability that *X* is the correct model given our initial (uncorrected) observation that model A was chosen as the best model with relative support *a*. The procedure accounts for difficulties differentiating between the models while also accounting for the fact that model A might have been chosen with strong support (*a* close to 1). We implemented the procedure by calculating the relative support for model A from each of the 100 PODs produced by each of our 11 models. For each model, we then calculated a probability density across the range of relative supports (0–1). We expected that models that were rarely confused with model A would have densities centred near zero relative support, while models similar to model A (including model A itself) would have densities centred closer to 1. We then evaluated these densities at the observed value (*a*) of relative support for model A (e.g. to evaluate the above equation).

We also assessed model fit against our observed data using tail-area probability, or *p*-value, tests for each summary statistic [27]. We used false discovery rate corrections on the statistics within each model [28], as well as the method of Ghirotto *et al.* [29] to combine *p*-values across statistics. When combining *p*-values, we assigned 0.001 to those *p*-values at the edge of the empirical distribution in order to avoid infinite values in the calculation. Low *p*-values indicated that the observed data were unlikely given a particular model. Finally, we performed a principal components analysis (PCA) on the summary statistics and plotted the observed data along with simulated datasets from the prior and from the posterior [27]. Ideally, the observed data would fall within the posterior for the selected model.

To test the accuracy of our parameter estimation procedure, we estimated model parameter for each POD *i* and compared the estimates against the known, true parameters (*θ _{i}*) used to simulate each POD. We then calculated the bias and root mean square error (RMSE) for each parameter across all

*n*PODs: and

We also calculated the factor 2 statistic, which is the proportion of estimated values that fall between 50 and 200% of the true value [30]. We used coverage to assess the accuracy of our confidence intervals, defined as the proportion of simulations that fall within the 50 or 90% HPD around the estimated parameter. Finally, we calculated the coefficient of determination (*r ^{2}*) for each parameter, which is the proportion of parameter variance that can be explained by the summary statistics [17].

## 3. Results

### (a) Model selection

The raw posterior probabilities from our ABC analysis suggested the highest support (0.70) for model 4b (table 3), the scenario that involved rapid growth of the VLC population after colonization from MQ, but that did not include ongoing migration from MQ or back-migration to MQ. We found lower support for model 6b (0.24), which did include migration back to MQ near the collapse of VLC. There was very little support for scenarios in which VLC and MQ were panmictic, where they were completely independent, where MQ was colonized by VLC or where VLC grew slowly.

However, when we simulated datasets to validate our model choice results, we found that distinguishing between certain models could be difficult, at least after such models had been fitted to the observed data (a more stringent test than used in many previous studies; electronic supplementary material, table S4). In particular, models 1, 3, 5, 6, 7, 5b and 7b were misidentified as a different model more than 50% of the time (see electronic supplementary material, table S4). Model 4b was the most accurately identified scenario, with 90% of simulations correctly chosen as model 4b. Most of the remaining model 4b simulations were mistaken for model 6b, a model that included 4b as a limiting case when back-migration was low. Encouragingly, however, and of particular relevance to our study, few scenarios were misidentified as model 4b with high support (e.g. greater than or equal to 0.70; electronic supplementary material, figure S1). The exception to this general rule was model 5b, which, like 6b, included 4b as a limiting case when the migration rate was low. We corrected for the possibility of misidentification by calculating for each model (table 3; electronic supplementary material, figure S1). We found that most of the posterior support, once corrected for errors in model choice, was split between models 4b and 5b, with lower support for models 6b and 7b. This result suggested that rapid growth soon after VLC colonization was highly likely, but determining whether or not migration continued after colonization could be more difficult. Back-migration to MQ near the collapse of VLC was less well supported (probability of 0.27 that the true scenario was models 6b or 7b).

Further model checking with tail probabilities indicated that at least one summary statistic failed the test (*p* < 0.05) for each of models 1–7, but none of the statistics failed for models 4b–7b (those with rapid growth; electronic supplementary material, table S5). The combined *p*-value statistic indicated particularly poor fits for models 1 and 2 (*p* < 0.001), but did not reject any of the other models. The PCA analysis on the summary statistics indicated poor fits for models 1–4, with the observed data clearly falling outside the simulations from the posterior (see electronic supplementary material, figure S3). The other models provided better fits to the data, according to this test.

Under alternative model priors or alternative statistical group definitions, model 4b continued to be selected as the best model (table 3). However, under lognormal priors for migration rate, the model choice procedure indicated lower support for model 4b than before, with the remaining support distributed primarily among models 5b, 6b and 7b. This result is to be expected because very low values of migration rates (close to 0) were commonly selected under lognormal priors, and so the functional differences among models 4b, 5b, 6b and 7b were smaller. Under wider priors for the size of the colonizing population, support for model 4b remained high (68%). Support also remained substantially higher for model 4b (51%) than for alternative models (less than or equal to 18%) with four statistical groups.

### (b) Parameter estimation

We then estimated the posterior distributions for both of the well-supported models (4b and 5b) and found similar parameter estimates from both (table 1 and figure 2). In particular, we found that MQ was about two orders of magnitude smaller than VLC, consistent with the lower genetic diversity in this population (table 1 and figure 2). Posterior parameter estimates for the migration rate between MQ and VLC also clarified why model 5b received similar support to model 4b: the migration rate posterior for model 5b was nearly 0 (95% HPD of 0–0.003; table 1 and figure 2). With a very low to zero migration rate, model 5b was functionally similar to model 4b, and the two were therefore difficult to distinguish. From the width of the posterior, it was clear that we did not have much power to detect the effective number of females that colonized VLC, though a founding size of less than 200 seals was unlikely (figure 2).

In testing the accuracy of our parameter estimates, we found that the relative bias and RMSE were quite low for the MQ and ancient VLC effective size (less than 50%), but substantially higher for the final VLC effective size with a positive bias of up to 100% (table 4). Our estimates of final VLC size, however, were in the right order of magnitude. This accuracy was also reflected in the factor 2 calculations, which were highest for MQ size and lowest for final VLC size. The coverage of our 50% and 90% HPDs were quite accurate, with mean values across parameters of 44 and 88%. The high values for our coefficients of determination suggested that our choices of summary statistics were informative. We concluded that our estimates of MQ and ancient VLC population sizes were most accurate, but we could only approximate the final VLC size to within an order of magnitude.

When we compared the fit of simulated datasets based on our posterior distributions, we found that they were a substantially better fit to the observed data than were the prior distributions (figure 3). However, the best model continued to predict haplotype diversity in the ancient VLC population (*H*_{VLC1}) that was somewhat too low, while the predicted number of segregating sites in the recent VLC population remained somewhat too high (*S*_{VLC2}). On the other hand, the key prediction of final haplotype number in the VLC was a strong match. The fit for model 4b was also substantially better than for a similar model that did not include rapid initial population growth in VLC (model 4). The latter model substantially under-predicted haplotype diversity in VLC (*H*_{VLC1} and *H*_{VLC2}), suggesting that rapid growth was important for the generation and maintenance of genetic diversity (see electronic supplementary material, figure S2).

## 4. Discussion

Our analysis indicated that the high levels of genetic diversity attained within the VLC southern elephant seal colony (table 2) probably arose largely by rapid mitochondrial control region *in situ* mutation as a function of early rapid population growth to a large effective population size, which effectively reduces the probability of extinction for new (single-copy) alleles [5]. This followed the release of a large expanse of suitable breeding habitat along the VLC, after the retreat of the Ross Sea Ice Sheet in the Early- to Mid-Holocene [8]. ABC provides strong support for the idea that MQ was the source of founders of the VLC colony, consistent with previous findings [9] (though the mtDNA marker provides long-term inference only about the movement of females). Further, our analysis did not support the idea that seals from VLC were purely a non-breeding moulting colony, that VLC was the source of the MQ population or that seals from the two locations were panmictic. We could also rule out any of the other major extant southern elephant seal breeding colonies as sources of either founders or observable levels of gene flow [9]. Thus, the significant levels of genetic diversity observed in the VLC population probably resulted from a combination of rapid substitution rate at the mitochondrial control region and rapid population growth to a large effective population size soon after colonization of newly available habitat resulting from Holocene climate change.

The model with the highest posterior probability is one of initial founding of VLC by a cohort of MQ seals (95% HPDs of female effective size: 157–953), followed by rapid population growth and a lack of gene flow between the two locations until final extinction of the VLC colony, some 7000 years later [9]. Interestingly, inclusion of an initial rapid growth phase would appear critical to fitting the models to our data (table 3 and figure 3; electronic supplementary material, figure S2). The next best-fit model is otherwise the same, also including a rapid growth phase, but distinguished by very low levels of ongoing migration from MQ to VLC (95% HPDs of migration rate 0–0.003). Although immigration from an unsampled ‘ghost’ population cannot be fully excluded, there are no apparent candidate populations not already directly or indirectly sampled, and data from de Bruyn *et al.* [9], which sampled all major extant populations, illustrate a clear connection between VLC and MQ in the mtDNA networks. One other possible interpretation is that migrants to VLC were from other now-extinct colonies, not from Macquarie. The most likely candidate, based on geographical proximity, is the population from the northwestern coast of Tasmania, extinct since prehistoric times [31], though 1500 km further from VLC than MQ. In this context, it seems unlikely that the VLC population would be founded by the Tasmanian (or the even more distant historical colony in New Zealand) rather than by the MQ population. The colonies on Juan Fernandez and St Helena were probably small and exterminated by hunting in the nineteenth and twentieth centuries. Both were very distant from VLC: over 4000 km to Juan Fernandez and over 6000 km to St Helena (the northernmost record for this species [32]). The chance that either of these would form a useful/relevant reference sample for this study seems very remote.

Simulated datasets based on posterior distributions fitted the observed data far better than the prior distributions (see electronic supplementary material, figure S2), including a strong fit between observed recent VLC haplotypic diversity and that predicted (figure 3). There were two exceptions, with the model prediction for haplotype diversity being low for the older VLC sample, and high for the number of segregating sites in the more recent VLC sample. The lower than expected observed number of segregating sites may simply reflect higher than expected variance of substitution rate among sites. While the best-supported model was one of no post-founder gene flow, one possible explanation for higher than expected haplotypic diversity early on is a post-founder ‘connection event(s)’ (see [33]). Although this may be consistent with our second (model 5b) and third (model 6b) best-fitting models (table 3), the level of migration indicated in those models is very small. Other possible explanations include stochastic sampling effects and sequence error in the older samples due to post-mortem damage. The latter is unlikely given the controls and replication undertaken to ensure accurate genotyping, including multiple replicate extractions through to sequencing, and cloning (see [9] for further details).

The key result remains the match between the very high haplotype diversity at VLC, including extensive novel diversity compared with the source population MQ [9], and model expectations derived using ABC and ancient DNA. While the generation of genetic diversity in this case is presumed to be neutral and the mutation rate relatively high (mitochondrial control region), our results suggest that given the right conditions, environmental change could generate the raw material for adaptive evolution (affecting multigenic phenotypic traits) through rapid population growth over a very short evolutionary time scale.

## Funding statement

M.L.P. was supported by National Science Foundation Graduate Research and David H. Smith Conservation Research fellowships. Research funding was provided by National Science Foundation grant no. PLR-0807292 and support from the Italian Antarctic Programme (PNRA).

## Acknowledgements

Two anonymous reviewers provided very helpful comments that improved the manuscript. We thank M. Blum for assistance with calculating the coefficient of determination and Oscar Gaggiotti for helpful comments that improved the manuscript.

- Received November 25, 2013.
- Accepted January 3, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.