Recent demographic histories and genetic diversity across pinnipeds are shaped by anthropogenic interactions and mediated by ecology and life-history

A central paradigm in conservation biology is that population bottlenecks reduce genetic diversity and negatively impact population viability and adaptive potential. In an era of unprecedented biodiversity loss and climate change, understanding both the determinants and consequences of bottlenecks in wild populations is therefore an increasingly important challenge. However, as most studies have focused on single species, the multitude of potential drivers and the consequences of bottlenecks remain elusive. Here, we used a comparative approach by integrating genetic data from over 11,000 individuals of 30 pinniped species with demographic, ecological and life history data to elucidate the consequences of large-scale commercial exploitation by 18th and 19th century sealers. We show that around one third of these species exhibit strong genetic signatures of recent population declines, with estimated bottleneck effective population sizes reflecting just a few tens of surviving individuals in the most extreme cases. Bottleneck strength was strongly associated with both breeding habitat and mating system variation, and together with global abundance explained a large proportion of the variation in genetic diversity across species. Overall, there was no relationship between bottleneck intensity and IUCN status, although three of the four most heavily bottlenecked species are currently endangered. Our study reveals an unforeseen interplay between anthropogenic exploitation, ecology, life history and demographic declines, sheds new light on the determinants of genetic diversity, and is consistent with the notion that both genetic and demographic factors influence population viability.


Introduction
Unravelling the demographic histories of species is a fundamental goal of population biology and has tremendous implications for understanding the genetic variability observed today 1,2 . Of particular interest are sharp reductions in the effective population size (Ne) known as population bottlenecks 3,4 , which may negatively impact the viability and adaptive evolutionary potential of species through a variety of stochastic demographic processes and the loss of genetic diversity [5][6][7][8] . Specifically, small bottlenecked populations have elevated levels of inbreeding and genetic drift, which decrease genetic variability and can lead to the fixation of mildly deleterious alleles and ultimately drive a vortex of extinction 6,[8][9][10] . Hence, investigating the bottleneck histories of wild populations and their determinants and consequences is more critical than ever before, as we live in an era where global anthropogenic alteration and destruction of natural habitats are driving species declines on an unprecedented scale 11,12 . Unfortunately, detailed information about past population declines across species is sparse because historical population size estimates are often either non-existent or highly uncertain 13,14 . A versatile solution for inferring population bottlenecks from a single sample of individuals is to compare levels of observed and expected genetic diversity, the latter of which can be simulated under virtually any demographic scenario based on the coalescent [15][16][17] . A variety of approaches based on this principle have been developed, one of the most widely used being the heterozygosity-excess test, which compares the heterozygosity of a panel of neutral genetic markers to the expectation in a stable population under mutation-drift equilibrium 18 . Although theoretically well grounded, these methods are highly sensitive to the assumed mutation model, which is seldom known 19 .
A more sophisticated framework for inferring demographic histories is coalescent-based Approximate Bayesian Computation (ABC) 20 . ABC has the compelling advantages of making it possible to (i) compare virtually any demographic scenario as long as it can be simulated, (ii) estimate key parameters of the model such as the bottleneck effective population size and (iii) incorporate uncertainty in the specification of models by using priors. Due to this flexibility, ABC has become a state of the art approach for inferring population bottlenecks as well as demographic histories in general [20][21][22][23][24][25][26] .
Although the widespread availability of neutral molecular markers such as microsatellites has facilitated numerous genetic studies of bottlenecks in wild populations, the vast majority of studies focused exclusively on single species and were confined to testing for the presence or absence of bottlenecks. We therefore know very little about the intensity of demographic declines and how these are influenced by anthropogenic impacts as well as by factors intrinsic to a given species. For example, species occupying breeding habitats that are more accessible to humans would be expected to be at higher risk of declines, while species with highly skewed mating systems tend to have lower effective population sizes 27 and might also experience stronger demographic declines as only a fraction of individuals contribute towards the genetic makeup of subsequent generations.
Consequently, in order to disentangle the forces shaping population bottlenecks, we need comparative studies incorporating genetic, ecological and life history data from multiple closely related species within a consistent analytical framework.
Another question that remains elusive due to a lack of comparative studies is to what extent recent bottlenecks have impacted the genetic diversity of wild populations. While a number of influential studies of heavily bottlenecked species have indeed found very low levels of genetic variability [28][29][30][31] others have reported unexpectedly high genetic variation after supposedly strong population declines 24,[32][33][34] . Hence, it is not yet clear how population size changes contribute towards one of the most fundamental questions in evolutionary genetics -how and why genetic diversity varies across species 2,[35][36][37] . To tackle this question, we need to compare closely related species because deeply divergent taxa vary so profoundly in their genetic diversity due to differences in their life-history strategies that any effects caused by variation in Ne will be hard to detect and decipher 36,37 .
Finally, the relative contributions of genetic diversity and demographic factors towards extinction risk remain unclear. While historically there has been a debate about the immediate importance of genetic factors towards species viability 5,7 there is now growing evidence that low genetic diversity increases extinction risk 8,38 and on a broader scale that threatened species show reduced diversity 7 . Nevertheless, due to a lack of studies measuring bottlenecks consistently across species, it remains an open question as to how the loss of genetic diversity caused by demographic declines ultimately translates into a species extinction risk, which is assessed by its International Union for Conservation of Nature (IUCN) status.
An outstanding opportunity to address these questions is provided by the pinnipeds, a clade of marine carnivores inhabiting nearly all marine environments ranging from the poles to the tropics and showing remarkable variation in their ecological and life-history adaptations 39 . Pinnipeds include some of the most extreme examples of commercial exploitation known to man, with several species including the northern elephant seal having been driven to the brink of extinction for their fur and blubber by 18 th to early 20 th century sealers 13 . By contrast, other pinniped species inhabiting pristine environments such as Antarctica have probably had very little contact with humans 13 . Hence, pinnipeds show large differences in their demographic histories within the highly constrained time window of commercial sealing and thereby represent a unique 'natural experiment' for exploring the causes and consequences of recent bottlenecks.
Here, we conducted a broad-scale comparative analysis of population bottlenecks using a combination of genetic, ecological and life-history data for 30 pinniped species. We inferred the strength of historical declines across species from the genetic data using two complimentary coalescent-based approaches, heterozygosityexcess and ABC. Heterozygosity-excess was used as a measure of the relative strength of recent population declines, while a consistent ABC framework was used to evaluate the probability of each species having experienced a severe bottleneck during the known timeframe of commercial exploitation, as well as to estimate relevant model parameters. Finally, we used Bayesian phylogenetic mixed models to investigate the potential causes and consequences of past bottlenecks while controlling for phylogenetic relatedness among species. We hypothesised that (i) extreme variation in the extent to which species were exploited by man should be reflected in their genetic bottleneck signatures; (ii) both breeding habitat and mating system should have an impact on the strength of bottleneck signatures across species; (iii) past bottlenecks should reduce contemporary genetic diversity; and (iv) heavily bottlenecked species with reduced genetic diversity will be more likely to be of conservation concern.

Results
Genetic data. We analysed a combination of published and newly generated microsatellite data from 30 pinniped species, with a median of 253 individuals and 14 loci per species (see Methods and Supplementary Table 1 for details). Measures of genetic diversity, standardised across datasets as the average per ten individuals, varied considerably across the pinniped phylogeny, with observed heterozygosity (Ho) and allelic richness (Ar) varying by over two and almost five-fold respectively across species (Supplementary Table 2).
Both of these measures were highly correlated (r = 0.92) and tended to be higher in ice breeding seals, intermediate in fur seals and sea lions, and substantially lower in a handful of species including northern elephant seals and monk seals (Fig 1A).
Bottleneck inference. We used two different coalescent-based approaches to infer the extent of recent population bottlenecks. First, the amount of heterozygosity-excess at selectively neutral loci such as microsatellites is an indicator of recent bottlenecks because during a population decline the number of alleles decreases faster than heterozygosity 3 . Recent bottlenecks therefore generate a transient excess of heterozygosity relative to a population at equilibrium with an equivalent number of alleles 18 . Here, we quantified the proportion of loci in heterozygosity-excess (prophet-exc) for each species, which was highly repeatable across a range of mutation models (see Methods and Supplementary Table 3). Consequently, we focused on a two-phase model with 80% single-step mutations (TPM80), which is broadly in line with mammalian mutation model estimates from the literature 40 as well as posterior estimates from our ABC analysis (Supplementary table 4). Fig. 1B shows a heatmap of prophet-exc across species, which is bounded between zero (all loci show heterozygosity-deficiency, an indicator of recent expansion) and one (all loci show heterozygosity-excess, an indicator of recent decline) whereby 0.5 is the expectation for a stable population. Considerable heterogeneity was found across species, with northern and southern elephant seals, grey seals, Guadalupe fur seals and Antarctic fur seals showing the strongest bottlenecks signals. By contrast, the majority of ice-breeding seals showed heterozygosity-deficiency, consistent with historical population expansions.
Second, we used ABC to select between a bottleneck and a neutral model as well as to estimate posterior distributions of relevant parameters. To optimally capture recent population size changes across species, we allowed Ne to vary over time in both models within realistic priors (see Methods) while the bottleneck model also included a severe decrease in Ne to below 800 during the time of peak sealing. ABC was clearly able to distinguish between the two models, with the posterior probability of correct model classification being 75% for the bottleneck model and 71% for the neutral model, and for every species the preferred model showed a    although the pre-bottleneck effective population size (Nehist) also had a prediction error below one in both models, visual inspection of the cross-validation results revealed high variation in the estimates and a systematic underestimation of larger Nehist values, so this parameter was not considered further.  Table 4B). Therefore, although studies of individual species are usually limited by uncertainty over the underlying mutation characteristics, our ABC analyses converged on similar estimates of mutation model and rate across species, allowing us to appropriately parameterise our bottleneck analyses.
To explore whether our results could be affected by population structure, we used STRUCTURE 41 to infer the most likely number of genetic clusters (K) across all datasets . For all of the species for which the best supported value of K was more than one (n = 12), we recalculated genetic summary statistics and repeated the bottleneck analyses based on individuals comprising the largest cluster. Using the largest genetic clusters did not appreciably affect our results, with repeatabilities for the genetic summary statistics and bottleneck signatures all being greater than 0.9 (see Supplementary table 6 for repeatabilities1 and Supplementary Fig. 3, which is virtually identical to Fig. 1).   prophet-exc and pbot (Fig. 4). A substantial 77% of the total variation in Ar was explained (Fig. 4C Fig 4A), decreased nearly five-fold from the species with the lowest pbot to the species with the highest pbot (β = -2.00, CI [-3.35, -0.74] Fig 4A) and was on average 27% higher in ice than in land-breeding seals (β = 1.70, CI [0.31, 3.29], Fig. 4A). Due to multicollinearity among the five predictor variables (Supplementary Table 7), standardized β estimates (Fig 4D) can be hard to interpret because of potential suppression effects 42 . This is reflected in the low unique R 2 values of the predictors relative to the marginal R 2 of the full model (Fig. 4C). However, the structure coefficients (Fig. 4E) Fig. 1 and Supplementary table S1.

Conservation status, bottleneck signatures and genetic diversity.
To investigate whether population bottlenecks and low genetic diversity are detrimental to species viability, we asked whether contemporary conservation status is related to both the strength of past bottlenecks and to Ar. Based on data from the IUCN red list 43 , we classified species into two groups; the first comprised species listed as 'least concern' while the second combined species listed as 'near threatened', 'vulnerable' or 'endangered' into a 'concern' category. Using a phylogenetic mixed model, we did not find any clear differences in either heterozygosity-excess or pbot with respect to conservation status (Fig. 6A, B). By contrast, average Ar was around 1.3 alleles lower in the concern category, although there was some uncertainty with the 95% credible interval of β ranging from -0.07 to 2.52. into either a concern or least concern category depending on their current IUCN status as described in the main text. Shown are the raw data for each category together with standard Tukey box plots for (A) Ar, (B) prophet-exc and (C) pbot. Marginal R 2 and standardised β estimates are shown for Bayesian phylogenetic mixed models with standardized predictors (see Materials and methods for details).

Discussion
To explore the interplay between historical demography, ecological and life-history variation, genetic diversity and conservation status, we used a comparative approach based on genetic data from over 80% of all extant pinniped species. To model bottleneck strength, we used two approaches that capture different but complementary facets of genetic diversity resulting from population bottlenecks. Using ABC, we contrasted a bottleneck model incorporating a severe decrease in Ne during the time of peak sealing in the 18 th and 19 th Centuries with a neutral model. The resulting bottleneck measure, pbot is the probability (relative to the neutral model) that a species' observed genetic diversity is similar to the diversity of a population that experienced a severe reduction in Ne below 800, and therefore provides an absolute bottleneck measure. By contrast, heterozygosity excess (prophet-exc) theoretically captures sudden recent reductions in Ne even in fairly large populations 18 , and therefore provides a relative bottleneck measure. Concretely, given the average sample size of individuals and loci used in this study, we would expect to detect an excess of heterozygosity at the majority of loci (i.e. prophet-exc > 0.5) when a 100-to 1000-fold reduction in Ne occurred within the last 4Ne generations, regardless of the magnitude of Ne (see simulations in 18 ).
ABC analysis supported the bottleneck model for more than a third of the species. The strongest bottlenecks (Nebot <50) were inferred for the northern elephant seal, a textbook example of a species that bounced back from the brink of extinction 44 , as well as for the two monk seals and the Saimaa ringed seal, species with very small geographic ranges and a long history of anthropogenic interaction 13 . Slightly weaker bottlenecks were estimated for seven further species including Antarctic and Guadalupe fur seals, both of which share a known history of commercial exploitation for their fur 13 . At the other end of the continuum, several Antarctic species that have not been commercially hunted such as crabeater and Weddell seals showed unequivocal support for the neutral model in line with expectations. Surprisingly, several otariid species known to have been hunted in the hundreds of thousands (e.g. South American sea lions) to millions (e.g. northern fur seals) did not show clear support for a bottleneck as strong as simulated in our analyses. This suggests that sufficiently large numbers of individuals must have survived despite extensive sealing, possibly on inaccessible shores or remote islands 45 .
We hypothesised that not all pinniped species were equally affected by commercial exploitation partly due to intrinsic differences relating to a species' ecology and life-history. In line with this, we found a strong influence of breeding habitat on bottleneck signatures, with both prophet-exc and pbot being higher in species that breed on land relative to those breeding on ice. A likely reason for this is that terrestrially breeding pinniped species were more profitable due to their generally higher population densities and accessibility, and therefore probably experienced more intense hunting. We also found that heterozygosity-excess was strongly linked to sexual size dimorphism (SSD), with highly polygynous species like elephant seals and some fur seals showing the strongest footprints of recent decline. While this could reflect the increased ease of exploitation and thus higher commercial value of species that predictably aggregate in very large numbers to breed, species with higher SSD also have highly skewed mating systems making them potentially more vulnerable to severe decreases in Ne when key males are taken out of the system. By contrast, we did not find an effect of SSD on the ABC bottleneck probability pbot, suggesting that although sexually dimorphic species experienced the greatest declines, these were not necessarily as severe as simulated in the ABC analysis (Ne < 800). This is probably because many species reached 'economic extinction' well above this threshold, when populations became too small to sustain the sealing industry.
Although vast numbers of species are declining globally at unprecedented rates 12 we still lack a clear understanding of how recent declines in Ne affect contemporary genetic diversity in wild populations 2,36 . Here, we explained a large proportion of the five-fold variation in allelic richness (Ar) observed from the most to the least diverse pinniped species. First, Ar was strongly associated with pbot but not with prophet-exc, in agreement with the theoretical expectation that populations have to decline to a very small Ne 3 , as was simulated in our ABC analysis, to lose a substantial proportion of their diversity. Second, we showed that global abundance across species was tightly linked to Ar, despite the likely impact of bottlenecks and the limited time-window for the recovery of genetic diversity. As differences in genetic diversity across species are largely determined by long-term Ne 2 , this implies that contemporary population sizes across pinnipeds must to some extent resemble patterns of historical abundance, and hence that many bottlenecked species have to a large extent rebounded to occupy their original niches. Third, Ar was higher in ice-breeding relative to land-breeding seals. However, a low unique R 2 of breeding habitat in our model suggests that this probably reflects the more intense bottleneck histories of land-breeding seals rather than a true ecological effect.
Finally, we compared genetic diversity and bottleneck strength between species that are currently classified by the IUCN as being of conservation concern versus those that are not. We found that Ar was on average around 22% lower in species within the concern category, consistent with previous evidence from a broad range of species 7 . While three out of the four pinniped species with the strongest estimated bottlenecks are currently listed as endangered, species from both categories did not overall differ in their bottleneck signatures. Our comparative study of population bottlenecks is therefore encouraging: population bottlenecks do not necessarily result in reduced genetic diversity, population viability and adaptive potential. As shown here, global bans on commercial sealing at the beginning of the 20 th Century allowed many surviving pinniped populations to recover in abundance. Those that have not sufficiently rebounded illustrate the two fundamental conservation challenges, especially as biodiversity loss and climate change continue at unprecedented rates: halting population declines and promoting population recovery.

Methods
Genetic data. Microsatellite data were obtained from a total of 30 pinniped species including three subspecies of ringed seal (summarised in Supplementary Table 1). We generated new data for five of these species as Phylogenetic, demographic, life history and conservation status data. Phylogenetic data were downloaded from the 10k trees website 46 and plotted using ggtree 47 . The three ringed seal subspecies were added according to their separation after the last ice age 48 . Demographic and life-history data for each species were obtained from 49 . While most data stayed untransformed, we calculated sexual size dimorphism (SSD) as the ratio of male to female body mass, and log-transformed abundance across species to account for the several orders of magnitude differences in population sizes. Data on conservation status were retrieved from the IUCN website (http://www.iucnredlist.org/, 2017) 43 .

Data cleaning and preliminary population genetic analyses.
In order to maximise data quality, we checked all datasets by eye and generated summary statistics and tables of allele counts to identify potentially erroneous genotypes including typographical or formatting errors. In ambiguous cases, we contacted the authors to verify the correct genotypes. As several of the datasets included samples from more than one geographical location, we used a Bayesian approach implemented in STRUCTURE version 2.3.4 41 to infer the most likely number of genetic clusters (K) across all datasets. For computational and practical reasons, we used the ParallelStructure package in R 50 to run these analyses on a computer cluster. For all of the species for which the best supported value of K was more than one, we recalculated genetic summary statistics and repeated the bottleneck analyses based on individuals comprising the largest cluster and calculated repeatabilities including confidence intervals for all variables using the rptGaussian function in the rptR package 51 . We also tested all loci from each dataset for deviations from Hardy-Weinberg equilibrium (HWE) using χ 2 and exact tests implemented in pegas 52 and applied Bonferroni correction to the resulting p-values. Overall, 6% of loci were found to deviate from HWE in both tests and as these are unlikely to affect our comparative analyses, we focused subsequently on the full datasets.
Genetic diversity statistics. In order to examine patterns of genetic diversity across species, we calculated observed heterozygosity (Ho) and allelic richness (Ar) with strataG {Archer:2017vj} as well as the proportion of low frequency alleles (LFA), defined as alleles with a frequency of <5%, using self-written code. For maximal comparability across species with different sample sizes, we randomly sampled ten individuals from each dataset 1000 times with replacement and calculated the corresponding mean and 95% confidence interval (CI) for each summary statistic.

Heterozygosity-excess.
We quantified heterozygosity-excess using the approach in 18 Table   1). The microsatellite mutation rate µ was refined after initial exploration and drawn from a uniform prior with µ ~ U [10 -5 , 10 -4 ] which lies within the range of current empirical estimates 40,58 . The mutation model was defined as a generalized stepwise mutation model with the geometric parameter GSMpar reflecting the proportion of multistep mutations, uniformly distributed from GSMpar ~U[0, 0.3]. The neutral model was defined with five parameters (Fig. 7b). Ne, Nehist, µ and GSMpar were specified with the same priors as previously defined for the bottleneck model and the time parameter corresponding to the historical population size thist was drawn from a uniform distribution ranging between 10 and 70 generations ago (thist ~U [10,70]). where k is the number of alleles. All statistics were calculated using a combination of functions from the strataG package and self-written code. For the ABC analysis, we used a tolerance threshold of 5 x 10 -4 , thereby retaining 5000 simulations with summary statistics closest to those of each empirical dataset. For estimating the posterior probability for each scenario and each species, we used the multinomial regression method 20,63 as implemented in the function postpr in the abc package 26 where the model indicator is the response variable of a polychotomous regression and the accepted summary statistics are the predictors. To construct posterior distributions from the accepted summary statistics for the model parameters, we a local linear regression approach 20 as implemented in the abc function of the abc package.

Evaluation of model specification and model fit via cross-validation.
We evaluated whether ABC can distinguish between our two models by performing a leave-one-out cross validation implemented by the cv4postpr function of the abc package. Here, the summary statistics of one of the existing 2 x 10 . simulations were considered as pseudo-observed data and classified into either the bottleneck or the neutral model using all of the remaining simulations. If the summary statistics can discriminate between the models, a large posterior probability should be assigned to the model that generated the pseudo-observed dataset. This was repeated 100 times and the resulting posterior probabilities for a given model were averaged to derive the rate of misclassification. We furthermore checked for each species that the preferred model for each species provided a good fit to the empirical data by conducting a formal hypothesis test using the gfit function in abc.
Specifically, we used the median distance between the accepted and observed summary statistics as a test statistic, whereby the null distribution was generated using summary statistics from the pseudo-observed datasets. Hence, a non-significant p-value indicates that the distance between the observed summary statistics and the accepted summary statistics is not larger than the expectation based on pseudo-observed data sets, i.e.
the assigned model provides a good fit to the observed data.
Evaluation of the accuracy of parameter estimates via cross-validation. In order to determine which parameters (i.e. population sizes, times and mutation rates and models) could be reliably estimated, we used leave-one-out cross validation implemented in the cv4abc function from the abc package to determine the accuracy of our ABC parameter estimates. For a randomly selected pseudo-observed dataset, parameters were estimated via ABC based on the remaining simulations using the rejection algorithm and a prediction error was calculated. This is possible because we know the "true" parameter values from which a given pseudoobserved dataset was simulated. This procedure was repeated 1000 times and a mean prediction error ranging between 0 and 1 was calculated, where 0 reflects perfect estimation and 1 means that the posterior estimate does not contain any information about the true parameter value 26 .
Bayesian phylogenetic mixed models. Finally, we used Bayesian phylogenetic mixed models in MCMCglmm 64 to evaluate the ecological and life-history variables affecting bottleneck strength and genetic diversity, and to test whether bottleneck history and genetic diversity are predictive of contemporary conservation status.
Details of all the models are given in the supplementary material. All of the response variables were modelled with Gaussian distributions, while the predictors were fitted as fixed effects and the phylogenetic covariance matrix as a random effect. Predictors in models containing binary fixed effects were standardised by two standard deviations to allow a direct comparison between the effect sizes 65 . In models without binary fixed effects, the predictor variables were standardised by one standard deviation. For all models, we report the marginal R 2 as in 66 . Some of the predictors in our models were inter-correlated and multicollinearity might lead to suppression effects and make the interpretation of regression coefficients difficult 42 . We therefore reported standardized β estimates, structure coefficients, ( % , ) and unique R 2 values for all variables in all models. The structure coefficients represent the correlation between a predictor and the fitted response of a model independent of the other predictors, and therefore reflect the direct contribution of a variable to that model. On the other hand, the unique R 2 is the difference between the marginal R 2 of a model including and a model excluding a predictor, which will be small when another predictor explains much of the same variation in the response 42 . All model estimates were presented as the posterior median and 95% credible intervals (CIs).
We used uninformative priors with a belief (shape) parameter v = 1 for the variance-covariance matrices of the random effects and inverse-Wishart priors with v = 0.002 for residual variances. For each model, three independent MCMC chains were run for 1,100,000 iterations, with a burn-in of 100,000 iterations and a thinning interval of 1000 iterations. Convergence was checked visually and by applying the Gelman-Rubin criterion to three independent chains. All of the upper 95% confidence limits of the potential scale inflation factors were below 1.05.
Data and code availability. All data wrangling steps and statistical analyses except for the heterozygosity-excess tests 53 were implemented in R 67 . The documented analysis pipeline along with the raw data can be accessed via GitHub (https://github.com/mastoffel/pinniped_bottlenecks) and is fully reproducible.