Demographic Reconstruction of Antarctic Fur Seals Supports the Krill Surplus Hypothesis.

Much debate surrounds the importance of top-down and bottom-up effects in the Southern Ocean, where the harvesting of over two million whales in the mid twentieth century is thought to have produced a massive surplus of Antarctic krill. This excess of krill may have allowed populations of other predators, such as seals and penguins, to increase, a top-down hypothesis known as the 'krill surplus hypothesis'. However, a lack of pre-whaling population baselines has made it challenging to investigate historical changes in the abundance of the major krill predators in relation to whaling. Therefore, we used reduced representation sequencing and a coalescent-based maximum composite likelihood approach to reconstruct the recent demographic history of the Antarctic fur seal, a pinniped that was hunted to the brink of extinction by 18th and 19th century sealers. In line with the known history of this species, we found support for a demographic model that included a substantial reduction in population size around the time period of sealing. Furthermore, maximum likelihood estimates from this model suggest that the recovered, post-sealing population at South Georgia may have been around two times larger than the pre-sealing population. Our findings lend support to the krill surplus hypothesis and illustrate the potential of genomic approaches to shed light on long-standing questions in population biology.

1. Introduction

Anthropogenic exploitation, particularly of ecologically important organisms, can have profound and often unexpected effects on natural ecosystems, influencing community structure, function, and productivity [1,2]. A prime example of this comes from the Southern Ocean, where uncontrolled exploitation of the great whales in the first half of the 20th century [3] resulted in significant ecosystem-level changes [4,5,6,7,8]. Large-scale commercial whaling commenced in the Southern Ocean in 1904 and peaked in 1930, when the annual catch reached almost 40,000 whales [3]. An estimated two million whales were taken between 1904 and 1990, including over 360,000 blue whales, 72,000 fin whales, 400,000 sperm whales, 200,000 humpback whales and 110,000 minke whales [9]. The blue whales alone accounted for around 40 million tonnes of biomass, equivalent to that of a billion people [10].

These hunted whales would have consumed vast amounts of Antarctic krill (Euphausia superba), a shrimp-like crustacean that forms the dominant prey of many mammals and birds of the Southern Ocean [8]. Consequently, whaling activities in the Southern Ocean are estimated to have produced a surplus of around 147–380 million tonnes of uneaten krill per year [6,11]. This competitive release of krill is thought to explain increases in populations of seals and penguins at South Georgia during the 1940s to 1970s, a top-down hypothesis referred to as the ‘krill surplus hypothesis’ [6]. However, krill abundance in the Scotia Sea has decreased substantially since the 1970s as the amount of sea ice has declined in response to rising global temperatures [12,13]. Consequently, many krill-dependent Antarctic predators are currently declining as a result of bottom-up changes to the Antarctic marine ecosystem [7,8,14,15].

While arguments have been made for both top-down and bottom-up effects playing important roles in the Southern Ocean ecosystem, their relative importance remains contentious [16,17]. In particular, ecosystem modelling studies have suggested that prey release is sufficient alone to explain 20th century trends in populations of some seals and penguins [4,5] and several empirical observations have lent support to the krill surplus hypothesis [8,18,19,20]. However, some authors have pointed out inconsistencies among species as well as incomplete niche overlap between, for example, resident penguin populations and migratory whale populations [21,22]. Arguably, the greatest barrier to understanding the importance of top-down and bottom-up effects in the Southern Ocean is a lack of pre-whaling population size baseline data [23]. This makes it challenging to evaluate whether the major consumers of krill increased as a consequence of the demise of the whales, which is important for understanding how these predators might respond to the subsequent recovery of certain whale species [16].

A compelling approach to circumvent this lack of historical baselines is to reconstruct recent population size changes of krill predators from molecular genetic data. This approach has so far only been attempted by a single study [23], which estimated the long-term coalescent effective population size (Ne) of the Antarctic minke whale from eleven nuclear markers and converted the resulting value into a census population size estimate. This fell within the range of several contemporary abundance estimates and was interpreted as meaning that the number of Antarctic minke whales was not unusually high after whaling. However, this study was criticised [24] because effective and census population sizes are not directly comparable, with Ne usually being much smaller [25], and because molecular diversity (θ) based estimates of the long-term coalescent Ne integrate information over long timeframes (in the order of 4Ne generations), so are unlikely to be strongly affected by recent anthropogenic perturbations [26].

Alternatively, changes in Ne can be inferred from the site frequency spectrum (SFS). The SFS is the distribution of allele frequencies of a given set of loci within a population [27] and ranges from rare ‘singletons’, in which an allele appears only once in the sample of individuals, to high frequency alleles that are carried by the majority of individuals. Importantly, the SFS is directly affected by a population’s demographic history; bottlenecks decrease the number of rare alleles [28], whereas population expansions lead to an excess [29,30]. Various methods have been developed for inferring demographic histories from the SFS [31,32,33,34] including the program fastsimcoal2, which uses a composite likelihood based framework to robustly infer demographic parameters under complex demographic scenarios [35]. Estimating the SFS requires high-resolution data from multiple individuals, so SFS-based demographic inference was historically limited to species for which large genomic datasets were available [36,37]. However, the emergence of reduced representation sequencing approaches such as restriction site associated DNA (RAD) sequencing [38] has meant that population genomic data can now be generated at reasonable cost for practically any species. Consequently, increasing numbers of studies are implementing SFS-based demographic reconstruction in wild populations [39,40,41,42].

The Antarctic fur seal (Arctocephalus gazella) is ideally suited to testing the krill surplus hypothesis using a population genomic approach. This polygynous and site-faithful pinniped [43,44,45,46] breeds on sub-Antarctic islands, with around 97% of the global population being concentrated around South Georgia [47]. Commercial sealing began shortly after the discovery of South Georgia by Captain James Cook in 1775 and continued unabated until the 1820s, by which point around 1.2 million Antarctic fur seals are believed to have been taken [48]. After the last commercial catch of 170 individuals at South Georgia in 1908, the species was considered all but extinct. For several decades, hardly any individuals were sighted ashore [49] until the Discovery expedition of 1936, when a small breeding population with 12 pups was observed at South Georgia [47]. This was followed by a period of explosive population growth during the 1960s and 1970s [50], which culminated in the South Georgia population reaching an estimated 2.7 million individuals in 1990 [47] and potentially as many as four to seven million individuals by the late 1990s [51]. However, it is difficult to gather reliable census data when population sizes are very large [51] and there is also some evidence to suggest that the 1990 census may have been conducted during an atypical breeding season [47]. Furthermore, the larger estimate [51] is based on data from an single breeding colony on Bird Island and it is unclear to what extent these data can be extrapolated to the whole of South Georgia.

Despite these uncertainties, the rapid mid-20th century growth of the Antarctic fur seal population at South Georgia has been described as ‘unprecedented in pinnipeds’ [52]. Temporal concordance between this explosive growth phase and the harvesting of the whales (Figure 1) is therefore considered a key piece of evidence supporting the krill surplus hypothesis [5]. Accordingly, it has been speculated that this release of krill increased the carrying capacity of the environment for Antarctic fur seals and allowed the post-sealing population at South Georgia to exceed its historical abundance [53,54]. However, census data from the pre-sealing population are lacking and the only available historical estimate of 2.5 million individuals, based on a reconstruction of the number of harvested pelts and a female-only, age-structured population model, has a very high level of associated uncertainty (95% CI = 1.5–3.5 million individuals) [55].

Figure 1

Graphical representation of temporal trends in Antarctic fur seal abundance and baleen whale harvesting. The Antarctic fur seal was heavily hunted in the late 18th and early 19th centuries and was considered virtually extinct by the early 1900s. Subsequently, explosive population growth during the second half of the 20th century coincided with the harvesting of the baleen whales [5]. Temporal changes in Antarctic fur seal abundance are reconstructed from empirical population size estimates from the scientific literature (purple points, [47,49,51,55,56,57]). Numbers of harvested baleen whales originate from [58]. Comparable whale census population size estimates from the same period are not available. Similarly, although data are available on the annual numbers of sealing vessels that visited South Georgia [59], the number of seals taken was often not recorded, meaning that it is not possible to depict temporal changes in the number of harvested seals. Original artwork by Elena Fissenewert.

Here, we used RAD sequencing and a coalescent-based maximum composite likelihood approach to reconstruct the recent demographic history of the Antarctic fur seal. Specifically, we used the empirical SFS together with coalescent simulations in fastsimcoal2 to estimate the likelihood of our data under two alternative demographic scenarios, one incorporating a severe reduction in Ne during the timeframe when sealing is known to have taken place (the ‘bottleneck model’), and the other assuming no recent changes in Ne (the ‘null model’). The best supported model was then used to estimate relevant parameters including Ne values before, during, and after the bottleneck. We hypothesised that Antarctic fur seals would show a genome-wide signature of a strong demographic reduction caused by sealing. In line with the krill surplus hypothesis, we further hypothesised that competitive prey release due to the harvesting of the whales may have facilitated the demographic recovery of the South Georgia population and allowed it to attain a larger size than was present before sealing.

2. Materials and Methods

2.1. RAD Sequencing Data

RAD sequencing data were generated for 70 Antarctic fur seal individuals sampled during 1996–2001 from Bird Island, South Georgia, as described by Humble et al. [60]. Our dataset includes previously published data for 57 individuals [60] plus unpublished data for a further 13 individuals. Briefly, skin plugs were collected in the field using piglet ear notching pliers and whole genomic DNA was extracted using a modified phenol-chloroform protocol [61]. RAD libraries were prepared following Etter et al. [62], with minor modifications as described in Humble et al. [60], and were 250 bp paired-end sequenced on an Illumina Hiseq 1500 (Hayward, CA, USA).

2.2. Bioinformatic Data Processing

Read quality was assessed using FastQC v0.11.9 and sequences were trimmed to 225 bp and demultiplexed using process_radtags in STACKS v1.41 [63]. The reads were then mapped to the Antarctic fur seal reference genome v1.4 [64] using BWA MEM v0.7.13 [65] with default parameters. The resulting SAM files were converted to BAM format, sorted by coordinates and indexed using SAMtools v1.11 [66]. The BAM files were then used as input for calculating genotype likelihoods and the SFS using -doSaf and -realSFS in ANGSD v0.935 [67]. Due to the absence of ancestral state information, we estimated the folded SFS. Individual genotype likelihoods were calculated assuming Hardy–Weinberg equilibrium using the algorithm provided by SAMtools (-doSaf 1, -GL 1). Only sites with a minimum mapping quality of 20 and a minimum qscore of 20 were included (-minMapQ 20, -minQ 20). The qscores around indels were adjusted and bad reads were discarded (-baq 1, -remove-bads 1). We only retained sites that were present in a minimum of 90% of all individuals with a minimum depth of coverage of five and maximum depth of coverage of 58 per individual (-minInd 63, -setMinDepth 315, -setMaxDepth 3600). To estimate the SFS using the -realSFS command, the maximum number of iterations of the EM algorithm was set to 1000.

2.3. Demographic Modelling

We first simulated SFS based on two alternative demographic scenarios. The bottleneck model simulated a reduction in Ne followed by exponential population growth. For this model, we estimated the pre-sealing effective population size (Nepre-sealing), the bottleneck effective population size (Nebot), the growth rate after the bottleneck (GR) and the post-sealing effective population size (Nepost-sealing). By contrast, the null model assumed no population growth or decline and estimated only Nepost-sealing. In both models, the defined initial search ranges for Nepost-sealing were log uniformly distributed between 5000 and 50,000 diploid individuals. For the bottleneck model, the defined initial search range for Nepre-sealing was uniformly distributed between 5000 and 50,000 individuals, while the defined initial search range for Nebot was uniformly distributed between 10 and 250 individuals. However, the composite maximum likelihood approach implemented in fastsimcoal2 uses these search ranges solely as starting values and parameter values can therefore exceed the upper limits [35]. Detailed sealing records [59] were used to simplify the bottleneck model by fixing the timing of the start and end of the bottleneck to 22 and 11 generations ago, respectively, assuming a generation time of ten years [68]. For both models, we simulated data under a fixed mutation rate of 2.5 × 10−8, in line with marine mammal mutation rate estimates from the literature [69]. The growth rate in the bottleneck model was defined as a complex parameter log (Nebot/Nepost-sealing)/11. A total of 100 replicate runs were performed for each model, including 100 estimation loops with 100,000 coalescent simulations. We did not include singletons in the simulations, as these have been found to be biased when sequence coverage is low [70]. Out of the 100 replicates for each model, the run with the highest maximum likelihood was retained. The best model was then determined based on the delta likelihood values (difference between the estimated and observed likelihoods) and Akaike’s information criterion (AIC) = 2 * ln(likelihood) + 2 * K where ln is the natural logarithm, likelihood is the maximum likelihood, and K is the number of parameters in the model. Moreover, to evaluate the fit of the best model to our data, we simulated 600 SFS based on the parameter estimates from the best model and visually compared the resulting SFS with the observed SFS.

Finally, we investigated the uncertainty of our parameter estimates using a non-parametric bootstrap approach. Specifically, we used ANGSD to bootstrap the SFS and to generate 600 SFS estimates based on data that were subsampled with replacement. This was implemented using the option-bootstrap within the ‘realSFS’ module of ANGSD. Then, for each of these 600 SFS, the parameters for the model were re-estimated based on 100 replicate runs, each including 100 estimation loops with 100,000 coalescent simulations. For each SFS, the run with the top maximum likelihood was retained and used for the bootstrap distribution. For each parameter, 95% confidence intervals (CIs) were then calculated based on the resulting 600 bootstrap estimates.

3. Results

RAD sequencing produced an average of 1,409,505 (range = 328,537–4,488,652) 250 bp paired-end reads per individual. For each individual, a mean of 99.4% of the reads were aligned to the reference genome (range = 95.9–99.8%). A total of 31,328,941 positions with estimated genotype likelihoods were used to estimate the SFS, of which 1.2% were polymorphic. In line with the known demographic history of the Antarctic fur seal, the bottleneck model achieved the highest support, having both a higher maximum likelihood and a lower AIC than the null model (Table 1). Furthermore, simulated SFS based on the best supported model showed a similar distribution to the observed SFS (Supplementary Figure S1), indicating that this model provides a good fit to the empirical data. Based on this model, Nepre-sealing, Nebot, and Nepost-sealing were estimated at 12,506, 534, and 29,319, respectively (Table 1). To investigate model uncertainty, we used a non-parametric bootstrapping approach to produce 600 site frequency spectra by subsampling the empirical SFS with replacement. Although there was some variation in the uncertainty of the estimated parameter values (95% CIs: Nepre-sealing = 10,589–12,760, Nebot = 452–615 and Nepost-sealing = 14,667–33,222), none of the 95% CIs overlapped (Figure 2a). Furthermore, the ratio of Nepost-sealing to Nepre-sealing averaged 2.0 (95% CI: 1.37–2.64) across all 600 simulations and was always greater than one (Figure 2b).

Table 1

Relative likelihoods of the two alternative demographic models together with AIC values and parameter estimates for the pre-sealing effective population size (Nepre-sealing), the bottleneck effective population size (Nebot), and the post-sealing effective population size (Nepost-sealing). See the Materials and methods for details.

Model	Max(log10(likelihood)) a	Number of Parameters	AIC	Nepre-Sealing	Nebot	Nepost-Sealing
Bottleneck model	−882,361.6	3	4,063,431	12,506	534	29,319
Null model	−882,378.8	1	4,063,506	–	–	14,789

a Based on the best likelihood among the 100 independent runs for each model.

Figure 2

Reconstruction of the recent demographic history of the Antarctic fur seal based on population genomic data. (a) Schematic of the best supported demographic model showing the estimated parameter values for Nepre-sealing, Nebot, and Nepost-sealing. Dashed grey lines represent the fixed values for the start and end of the bottleneck. The priors for the estimated parameters are described in the Materials and methods. For all three parameters, distributions of 600 non-parametric bootstrap estimates are plotted together with box plots showing median values, and 25th and 75th percentiles, with the whiskers representing 95% confidence intervals. The Ne estimates are plotted on the same scale along the x-axis for comparability, while the distributions of Nepre-sealing and Nebot can be seen in the zoom-ins. (b) Distribution of the ratio of Nepost-sealing to Nepre-sealing over all 600 bootstrapped datasets. The box plot above the bar chart shows the median value, and 25th and 75th percentiles, with the whiskers representing the 95% confidence interval.

4. Discussion

We used RAD sequencing to infer recent changes in the size of the Antarctic fur seal population at South Georgia. In line with expectations based on known history of this species, our best supported model incorporated a recent bottleneck followed by exponential population growth. Furthermore, maximum likelihood estimates from this model suggest that the post-sealing Ne may have been around twice as large as the pre-sealing Ne. Our study builds upon ecosystem modelling [4,5] and empirical studies of multiple krill consumers [8,18,19,20] that lend support to the krill surplus hypothesis.

4.1. Bottleneck Inference

Our best supported model included a substantial demographic reduction during the known time period of commercial sealing. However, the maximum likelihood estimate of Nebot was both larger (534 versus ~150–300) and more precise than previous estimates based on microsatellites [59,71,72]. This probably reflects differences in both the genetic markers themselves and the analytical approach. One key methodological difference is that previous studies used approximate Bayesian computation, which produces estimates that are constrained by pre-defined priors, whereas fastsimcoal uses search ranges as starting values but the resulting parameter estimates can exceed the upper limits [35]. Accordingly, an exploratory analysis by Hoffman et al. [59] produced a more comparable Nebot estimate of around 700 when the priors for this parameter were relaxed.

An Nebot of 534 equates to a census population size substantially larger than Larsen’s [49] estimate of 30 individuals at South Georgia in 1911. Furthermore, mammalian effective population sizes are usually several times smaller than census population sizes [25], implying that thousands of animals may have escaped sealing at South Georgia. This is consistent with previous observations of high genetic diversity [71] and rapidly decaying linkage disequilibrium across the genome [60] in the South Georgia population. Furthermore, recent studies have uncovered a strong population genetic structure across the species’ circumpolar range [60,72,73,74], implying that at least four refugial populations survived commercial sealing. Thus, it appears that the species as a whole may have been more resilient to commercial exploitation than was previously believed, probably because small remnants of historically much larger populations were able to persist in a handful of remote and inaccessible locations [75].

4.2. Pre- and Post-Sealing Effective Population Sizes

A lack of pre-whaling baseline data on seal and penguin populations has made it challenging to evaluate the ecosystem-level consequences of historical whaling [23]. We circumvented this issue by reconstructing temporal changes in Antarctic fur seal abundance from RAD sequencing data. This approach has the advantage of producing Ne estimates that are directly comparable between different time points. In contrast to previous attempts at demographic reconstruction using microsatellites [59,71,72], we were able to precisely estimate Nepre-sealing, as indicated by the relatively narrow range of parameter estimates resulting from the non-parametric bootstrapping (95% CI = 10,589–12,760). Our estimate of the post-sealing Ne was somewhat less precise, with the 95% CI ranging from 14,667 to 33,222. This is probably because the population boom was short lived, lasting for only a few generations, which would not have been long enough to fully erase the imprint of the long-term Ne on the SFS. Nevertheless, the 95% CIs of Nepre-sealing and Nepost-sealing did not overlap and the ratio of Nepost-sealing to Nepre-sealing was consistently greater than one across all 600 bootstrapped datasets. Consequently, although there is some uncertainty associated with our parameter estimates, statistically speaking, there is a low probability that Nepost-sealing and Nepre-sealing are the same.

We also recognize that our Ne estimates are small in comparison to available census size estimates, which run into the millions of individuals (e.g., [47,51]). However, this is to be expected given that Antarctic fur seals exhibit strong polygyny [43,46], natal site fidelity [45], and population structure [60,72,73], all traits that are expected to reduce the ratio of the effective to census population size. Accordingly, our Ne estimates are comparable to other estimates from the literature [76] including genomic estimates of the long-term coalescent Ne from several pinniped species [26].

Some authors have argued that competitive prey release may have allowed the Antarctic fur seal population at South Georgia to increase beyond its pre-exploitation size [53,54]. However, these arguments were largely motivated by anecdotal observations, such as heavy damage to the tussock grass in the mid-1980s, which Bonner [53] described as ‘a new phenomenon associated with high density occupation by the recovering population’. While these kinds of observations can be prone to confirmation bias, our results suggest that the post-sealing population at South Georgia may indeed have been larger than the pre-sealing population. A similar conclusion was also reached by the authors of a paleolimnological study at Signy Island, which found over four times more Antarctic fur seal hairs in lake sediments from the 1990s than at any point during the past 6000 years [18]. While both of these studies ostensibly support the krill hypothesis, it is important to recognize that the density of Antarctic krill had already begun to decline before the fur seal population reached its peak in the 1990s [12]. One possible explanation for this lag [5] could be that the biomass of krill at South Georgia was sufficient to support the growing fur seal population until the 1990s, by which point the krill surplus would have come to an end (see below).

4.3. The Role of Bottom-Up Effects

Although it has been argued that the krill surplus hypothesis provides a valid explanation for increases in populations of penguins and seals during the mid-20th century [4], it cannot account for more recent declines of many Antarctic predators [4,8,10,14,22]. In practice, the excess krill biomass appears to have been eroded during the last quarter of the 20th century by a decline in the primary productivity of the Southern Ocean [4]. This has been linked to a reduction in the bioavailability of whale-recycled iron [4,10,11], as well as to climate-driven declines in sea ice extent and krill density in the Western Antarctic Peninsula and Scotia Sea [8,12,13,16]. In the case of the Antarctic fur seal, the number of females breeding at a long-term study population on Bird Island shows tight linkage to local krill availability [14] and there has also been a recent switch from positive to negative dependent pup mortality [77]. Consequently, large-scale changes in krill biomass may help to explain both the rapid post-sealing recovery of the fur seal population and its subsequent, ongoing decline.

4.4. Caveats

Although our study illustrates the utility of coalescent based simulation for inferring recent demographic histories, several important limitations are associated with this approach. First, computation of the SFS from genomic data can be challenging, particularly for low coverage data where inaccurate genotype calls can lead to biased estimation of the SFS [78,79,80]. To compensate for this, we implemented a genotype likelihood approach that incorporates uncertainty due to sequencing errors, variation in the depth of coverage and alignment quality [67]. Additionally, we excluded singletons from our analysis as these can be especially error prone in low coverage datasets [81]. Finally, non-parametric bootstrapping of the SFS allowed us to quantify the uncertainty associated with our maximum likelihood estimates. This was generally quite low, suggesting that our results are reasonably robust.

Second, alternative demographic trajectories can potentially produce the same SFS in a single population [82]. This is a limitation inherent to all demographic inference approaches in which demographic scenarios must be pre-specified. However, the recent history of the Antarctic fur seal is extraordinarily well documented [48,59,83,84] so we believe that our bottleneck model is a good representation of the true demographic history both in terms of the time parameters and population size priors. However, we cannot exclude the possibility that an alternative (although arguably less parsimonious) scenario might fit the data equally well.

Third, our study focused on a single species out of a suite of Southern Ocean predators. Thus, although our results are consistent with the krill surplus hypothesis, they do not allow us to infer causality. In the future, it would be interesting to conduct parallel analyses of SFS from multiple Southern Ocean predators including seals, penguins and smaller whale species. Although this approach would entail significant challenges in terms of data collection, it would allow the investigation of competitive release via the analysis of concurrent demographic patterns across multiple coexisting predators.

5. Conclusions

We used RAD sequencing to reconstruct the recent demographic history of the Antarctic fur seal population at South Georgia. We found evidence of a severe demographic reduction due to commercial sealing. Furthermore, maximum likelihood estimates of the pre- and post-sealing Ne from our best supported model were consistent with independent observations suggesting that prey release may have facilitated a rapid demographic increase in the recovering Antarctic fur seal population at South Georgia. Our results build upon previous modelling and empirical studies supporting the krill surplus hypothesis, although ongoing climate change has been linked to more recent declines of many Southern Ocean predators including Antarctic fur seals.