Risk for overexploiting a seemingly stable seal population: in ﬂ uence of multiple stressors and hunting

. Conservation efforts have mainly been focused on depleted species or populations, but many formerly reduced marine mammal populations have recovered to historical abundances. This calls for new management strategies and new models for ecological risk assessment that incorporate local density dependence and multiple environmental stressors. The harbor seal metapopulation in Swedish and Danish waters has increased from about 2500 to 25,000 over the past 40 yr. Trend analysis based on aerial survey data and somatic growth curves indicates that the population is close to carrying capacity. We performed a population viability analysis based on realistic life history parameters and investigated a range of potential scenarios caused by future stressors. If the population is able to resume its high intrinsic rate of increase at about 11% annually, when pushed down below carrying capacity, it can also sustain additional mortality such as modest hunting and infrequent epizootics. However, if xenobiotics will cause even a slight reduction in average fecundity, the population becomes signi ﬁ cantly more vulnerable. In the absence of epizootics, and given full reproductive capacity, hunting of a few hundred animals annually is not harmful to the long-term persistence of the population. Nevertheless, a slight decrease in growth potential, for exam-ple, caused by exposure to endocrine disruptors, makes even limited hunting risky. Our study shows how an apparently stable and abundant marine mammal population can be close to a point of rapid population decline. Thus, careful monitoring of population size, growth rate, health, and exposure to xenobiotics as well as recording of the age and sex structure of the hunt is required to avoid repeating the history of overexploitation and another population collapse. In the Anholt subpopulation, quasi-extinction is likely even under moderate levels of stress. In the Koster subpopulation, quasi-extinction is likely only under high levels of stress.


INTRODUCTION
The harbor seals (Phoca vitulina) studied here share a history of overexploitation with a majority of the other marine mammal populations worldwide. Many were hunted to near extinction and a few have faced true extinction (Magera et al. 2013, Lotze et al. 2017). However, after decades of active conservation measures and banned hunting, many marine mammal populations have recovered to historical numbers and a new era has begun with other challenges (Valdivia et al. 2019). Species that were once rare are seen as competitors and invaders when they resume natural abundances, demanding changes in the perception of the general public in order to avoid a second wave of overexploitation (Roman et al. 2015).
Historically, the hunting of marine mammals was motivated as a resource for food, blubber, and oil, and sometimes performed as a bounty hunt or culling campaign with the aim to reduce population size and thereby competition or interference with fisheries H€ ark€ onen 1999, Olsen et al. 2018). Hunting aiming at reducing populations is more efficient (and risky) in K-strategists such as the slowly reproducing seal populations as compared to many other animals with higher reproductive rates. The generation time for seals is 10-15 yr, and they give birth first at the age of 4-6 yr, and only bear one pup annually. These biological constraints limit the long-term maximum population growth rate to 10-12% in fully healthy seal populations (H€ ark€ onen et al. 2002). The low intrinsic rate of increase in combination with difficulties in estimating seal abundance in the field makes them vulnerable to overexploitation. Examples of marine mammal species that were strongly affected by overexploitation include the northern elephant seal (Mirounga angustirostris; Cooper and Stewart 1983), several species of monk seals (e.g., Monachus schauinslandi; Reijnders et al. 1993), and the Saimaa ringed seal (Pusa hispida saimensis; Kokko et al. 1999). The Caspian seal (Pusa caspica) population was reduced from 1.2 million at the end of the 19th century to 100,000 in 2005 (Harkonen et al. 2008). Populations of harbor seals and gray seals were extirpated from the southern Baltic and the continental North Sea coast (Reijnders et al. 1993(Reijnders et al. , H€ ark€ onen et al. 2007, and two species, the Caribbean monk seal (Monachus tropicalis) and the Japanese sea lion (Zalophus japonicus), are now extinct as a consequence of excessive hunting (McClenachan and Cooper 2008). History has shown that hunting or culling has been a very important stressor for most seal populations.
Several marine mammal populations are now recovering, and some have reached historical population sizes. The return from near extinction of humpback whales (Megaptera novaeangliae) is one such success story (Bejder et al. 2016); northern elephant seals (M. angustirostris) and Atlantic and Baltic gray seals (Halichoerus grypus) have also shown remarkable recoveries , Neuenhoff et al. 2018. When populations are devoid of pressures and resume their maximum growth potentials, numbers can increase exponentially. However, any environmental factor hampering the intrinsic growth rate quickly makes marine mammal populations vulnerable to increased mortality in the form of bycatches or hunting. This is reflected by the fate of the Baltic ringed and gray seal populations, where hunting caused a reduction in population size from 200,000 and 90,000 in the beginning of the 20th century to 25,000 and 20,000, respectively, in the 1940s (Harding and H€ ark€ onen 1999). Hunting continued after the 1940s and both populations crashed below 5000 seals when polychlorinated biphenyl (PCB) pollution led to sterility (Bergman 1999), which effectively decreased population growth rates. This combination of bycatches, hunting, sterility, and a pollution-induced disease complex (Bergman and Olsson 1985) almost eradicated the populations (Harding and H€ ark€ onen 1999).
The harbor seal is widespread across the Northern Hemisphere and is the most common seal species along the Swedish and Danish coasts of the Kattegat and the Skagerrak. Harbor seals show strong site fidelity to their breeding sites and form structured populations with many genetically differentiated subpopulations (Stanley et al. 1996, Goodman 1998, Olsen et al. 2014. The Swedish-Danish Skagerrak and Kattegat metapopulation is also genetically differentiated from the neighboring North Sea and Baltic Sea populations, and has been exposed to dramatic demographic fluctuations with historic culling, periods of exponential growth, and four events of increased mortality due to epidemic diseases within the last three decades (H€ ark€ onen et al. 2006(H€ ark€ onen et al. , Olsen et al. 2010). The population is now approaching its carrying capacity, and responding to a high density situation with reduced somatic growth and adult body size (Harding et al. 2018). Concurrently, interactions between seals and local fisheries have led to demands for increased hunting quotas (Olsen et al. 2010, Olsen et al. 2018). However, changes in seal populations caused by anthropogenic factors can result in a cascade of effects across trophic levels (Bowen 1997, Heithaus et al. 2008, Buren et al. 2014, potentially disturbing the ecosystem dynamics. This has been observed in the past when fish populations (e.g., cod, herring and sprat) changed following changes in seal populations and commercial fishing ( € Osterblom et al. 2007; T. H€ ark€ onen and K. C. Harding, unpublished manuscript). Because of that, modeling the dynamics of seal populations under different stressors is important for predicting possible impacts on the ecosystem as a whole.
Population viability analysis (PVA) is a quantitative approach for assessing causal factors affecting the population viability of a species (Beissinger and Westphal 1998) and is surprisingly accurate when retrospectively compared to empirical data (Brook et al. 2000). The construction of a carefully parameterized PVA model for pinnipeds is one step toward a better understanding of the sensitivity of the population to environmental variation and demographic changes. A PVA modeling approach can be used to evaluate the relative importance of different stressors and to investigate future scenarios of population management (Boyce 1992, Burgman et al. 1993, Gulve and Ebenhard 2000, Beissinger and McCullough 2002, Shaffer et al. 2002. Earlier risk analyses of the harbor seal population in the Kattegat-Skagerrak were carried out to evaluate the relative importance of the frequency and mortality in epizootics and host immunogenetics , 2003, Garnier et al. 2014). However, these models ignore density dependence and spatial structure, which are important factors regulating population growth.
Here we investigate the viability of the Swedish-Danish harbor seal population and elaborate several scenarios of single and combined effects of hunting, pollutants and outbreaks of infectious diseases. Data from mark-recapture studies, autopsies, satellite tagging, and longterm monitoring of population abundance are used as input values to produce the most detailed population model for harbor seals available. We also take geographic population structure into account by modeling the dynamics of each subpopulation within the study area separately. The findings can contribute to evaluations of the status of the marine mammal populations according to the EU Habitats Directive and International Union for Conservation of Nature (IUCN) criteria, and can be used to guide management decisions to avoid future overexploitation.

Study population
Aerial surveys of all haul-out sites in the Kattegat-Skagerrak harbor seal population ( Fig. 1) have been conducted annually during the peak haul-out in August, from 1979. From 1993, surveys were only conducted biannually of Danish sites. The proportion of the population hauled out on land during surveys is estimated at 65% (Heide-Jørgensen and H€ ark€ onen 1988, H€ ark€ onen et al. 2002). The entire area is surveyed three times annually, and the maximum count is used for trend analyses (H€ ark€ onen et al. 2006(H€ ark€ onen et al. , Galatius et al. 2014. A matrix of migrations to nearest neighbors was based on population genetic data, where 3-5% percent of the population migrate to the nearest colonies annually ( Fig. 1; Olsen et al. 2014).

Parameter values and initial conditions
Life history parameters are well known from extensive sampling and branding studies of the population, and we used parameter values given in H€ ark€ onen et al. (2002) for the baseline scenario. Survival of pups (at ages 0-1) was set at 75%, subadult survival (ages 1-4) at 89%, and adult survival (4+) was 95%. Age-specific fertility rates, that is, female pups per surviving female, were set at 17% for age 4, 33% for age 5, 47% for ages 6-27, and 35% for ages 28-38 (Heide-Jørgensen andH€ ark€ onen 1988, H€ ark€ onen et al. 2002). To account for stochasticity, a 5% random variation following a beta distribution (a = b = 5) was added every year. We parameterized a Leslie matrix and estimated the asymptotic population growth rate from the dominant eigenvalue, denoted k 1 and the stable age distribution, which is given by the corresponding eigenvector (Leslie 1948). The carrying capacity for the metapopulation was set at 35,000 individuals as a best guess motivated by several lines of evidence: historical hunting records (Heide-Jørgensen and H€ ark€ onen 1988) and recent signs of impaired somatic growth (Harding et al. 2018). Local carrying capacities for the subpopulations were assumed to be in proportion to current subpopulation size; thus, all subpopulations had the same relative capacity for increase from the current state (Table 1).

Population projections
An age-specific, stochastic metapopulation model for the harbor seal population was developed, and future population trajectories were simulated in the software MATLAB (see Appendix S1). Each scenario (100 yr each) was simulated 100 times, and the average population growth rate, final average population size, and quasi-extinction risk were recorded. Starting values for population size of each colony were taken from the most recent aerial survey (Table 1). The viability of the metapopulation was simulated for a range of future scenarios. A theoretical baseline scenario without any stressors was used as a control after which impacts from single stressors were investigated successively: hunting quotas, epizootic outbreaks, and reduction in fecundity due to xenobiotics. Subsequently, effects from combined stressors were systematically analyzed. The quasi-extinction threshold was set at one hundred females for each subpopulation. This sets the functional quasi-extinction to 1100 females, evenly spread among subpopulations, for the entire metapopulation in the Kattegat-Skagerrak. Quasi-extinction was defined to have occurred if the final population size was below the threshold at the end of the given timespan of the simulation (100 yr). Thus, single projections that reach below the threshold during the simulations are able to increase again v www.esajournals.org through immigration and/or chance events such as a sequence of favorable conditions. The proportion of all projections below the threshold at the end of the simulation for each subpopulation represents the local quasi-extinction risk.

Environmental stressors: hunting, epizootic outbreaks, and pollutants
We investigated population viability under four different hunting regimes with annual hunting quotas (Q) of 1%, 2.5%, 5%, and 10% of the population size. We first elaborate an age-unbiased hunting regime, where age classes are hunted in proportion to their sizes, and then investigate how a hunt biased toward pups, subadults, or adults can affect population viability. In the age-unbiased scenarios, individuals were first randomly selected from the total metapopulation and removed from their respective age classes and subpopulation, that is, hunting in each subpopulation is proportional to the subpopulation size (it also implies that hunting pressure increases on the remaining subpopulations when some go extinct). In the age-biased scenarios, individuals in the target age class are initially removed with a probability of 0.95, until the quota is filled, and if this is not sufficient to reach the hunting quota, randomly selected individuals of all age classes are harvested. Consequently, if subpopulations have different age structures, the subpopulations with most individuals in the target age class can be exposed to a proportionally higher hunting pressure.
The effect of a survey delay (D) of five or ten years was also analyzed, since the number of seals hunted any given year is calculated as a percentage of the population size a number of delay years ago due to variance in survey data (Svensson et al. 2011). In order to make a fair comparison of the effect of different survey delays, hunting was applied only after 10 yr in the simulated time in all simulations with hunting (Q > 0).
Repeated outbreaks of pathogens such as the phocine distemper virus (PDV;, influenza A (2014), and unknown pathogens (2007 and 2011) have plagued the Swedish-Danish harbor seal population (Osterhaus et al. 1988, Zohari et al. 2014). In the two PDV epizootics, mortality rates up to 66% were documented (H€ ark€ onen et al. 2002). We analyzed the impact of epizootic outbreaks by including a probabilistic function determining the probability of an epizootic outbreak to happen (E) and the mortality rate. Mortality caused by epizootics was based on the age distribution of dead individuals during the reported outbreaks in nature (H€ ark€ onen et al. 2007). Similarly to hunting mortality, individuals were randomly selected in the total metapopulation and removed from their respective subpopulations, with the maximum number of dead individuals not exceeding 65% of the total metapopulation size, mimicking the mortality caused by PDV. We tested effects of the following frequencies of epizootic outbreaks: 0.03, 0.07, and 0.10 (E; annual probability of occurrence).
We also tested effects of endocrine-disrupting pollutants on the fate of the metapopulation. Several types of pollutants such as PCBs and dichlorodiphenyl-trichloroethane (DDT) are known to bioaccumulate in seals and to impair reproductive health (Drescher et al. 1977, Bredhult et al. 2008, B€ acklin et al. 2011. In order to mimic the effect of an endocrine-disrupting substance , Bergman 1999, we reduced the average fecundity (F) by a factor 0.1-0.5 in the Notes: Column 2 shows the best estimate of true population size at each colony, that is, maximum numbers of observed seals corrected by the 65% haul-out fraction on land during surveys. Female population sizes are estimated to be 48% of the total population. Local carrying capacities are estimated by assuming a total carrying capacity is 35,000 for the entire Kattegat-Skagerrak metapopulation, divided by the relative size of each subpopulation at present. v www.esajournals.org simulations and assessed the impact of this reduction.
Since natural populations are exposed to many stressors simultaneously, we simulated several combinations of parameter values for hunting quota (Q) and delay (D), probability of epizootic outbreaks (E), and reduction in fecundity (F) due to pollutants (Table 2). Hunting quotas of Q = {0.01, 0.025, 0.05, 0.10} were tested in combination with hunting delays of D = {0, 5, 10}, probabilities of epizootic outbreaks of E = {0.03, 0.07, 0.10}, and fecundity of F = {0.5, 0.7, 0.9, 1.0} ( Table 2). The quasi-extinction risk and the longterm average population size after 100 yr for each subpopulation were recorded.

Multistressor sensitivity analysis
We investigated the sensitivity of quasi-extinction risk estimates to changes in the parameter values of stressors by elaborating the full range of possibilities (Naujokaitis-Lewis et al. 2009).
This analysis was done for all subpopulations, but we selected the most resilient (Koster), an intermediate (Hesselø) and the most vulnerable (Anholt) subpopulations to illustrate our results. Subsequently, the influence of age-biased hunting was investigated by targeting specific age classes, as described above. We tested hunting bias toward pups (age classes 0-1), subadults (age classes 1-4), and adults (age classes 5-38). Hunting quotas were always assumed to be reached. All simulations and analyses were conducted in MATLAB. In order to facilitate further elaboration of the PVA model and to facilitate applications to future management needs, the MATLAB code containing all the simulations and figures is provided in Appendix S1.

Growth trends in survey data
Growth functions were fitted to survey data in three different time periods in both sea regions (Fig. 2). In period 1 (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987), the Skagerrak Note: N 100 values give the population size after the 100-yr-long simulation period for the entire metapopulation (both females and males included). v www.esajournals.org population grew exponentially by 13.3% and the Kattegat grew by 12.7% in this early time period when the population was recovering from hunting. In period 2 (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), after the first PDV mass mortality in 1988, the population increased once more at an exponential rate, with an estimated growth of 15.2% in the Skagerrak and 9.6% in the Kattegat. Growth rates can exceed 12% only temporarily, due to skewed age structure after the PDV epizootics, and eventually reach stable age structure and the associated intrinsic rates of increase stabilize around 12% in the absence of environmental stressors (H€ ark€ onen et al. 2002).
If an exponential function is assumed also for the last period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), the average growth rates are 6.5% and 5% for the Skagerrak and the Kattegat, respectively. However, several lines of evidence suggest a recent decline in fecundity and survival (Harding et al. 2018), so logistic functions were also applied, with a predicted asymptote for the Skagerrak estimated at 8322 counted seals, while for the Kattegat, the asymptote is estimated to 10,140 counted seals. Accounting for the haul-out behavior, that is, only~65% of the population is hauled out during surveys (Heide-Jørgensen and H€ ark€ onen 1988), the logistic models suggest that the total carrying capacity for both areas can be around 30,000. However, given the short data series, the confidence limits are wide (95% CI 1977, 14,667 and 8230, 12,043 for Skagerrak and Kattegat, respectively) and the results can only be interpreted as an indication.

Simulation outcomes
In the theoretical best case scenario in the absence of environmental stressors, all local harbor seal populations reach their carrying capacities within 50 yr, when the asymptotic population size of the entire metapopulation approaches 34,935 seals (including males and v www.esajournals.org females; Fig. 3). Final average population sizes tended to decrease when multiple stressors were present, with more vulnerable subpopulations frequently going quasi-extinct (Appendix S1: Fig. S1).
Single stressor scenarios: hunting, epizootics, reproductive health All three single stressors, that is, hunting with or without a delay, mortality by epizootic outbreaks, and reduced fecundity, resulted in lowered asymptotic population sizes over the simulation period of 100 yr. However, the impact of these stressors differed in the magnitude of the reduction in population size and increase in quasi-extinction probability. While a low hunting quota (Q = 0.01) did not increase quasi-extinction probability, moderate (Q = 0.05) and high (Q = 0.10) quotas significantly increased quasiextinction probabilities of more vulnerable subpopulations in particular (Fig. 4), and a delay in quota calculation during moderate hunting did not affect quasi-extinction probability substantially for single stressor scenarios (Appendix S1: Fig. S2). In the most resilient of the subpopulations (Koster), the probability of quasi-extinction did not increase even under high hunting quotas during the 100-yr simulation period, but the population trends show steep declines in population sizes leading to the extinction of the Koster subpopulation if the hunting regime is maintained beyond the 100 yr of the simulations. Thus, a ten percent hunting quota was considered unrealistically high and was not further investigated.
When populations are exposed to pollutants decreasing female fecundity, population growth rates and asymptotic population sizes are reduced, but quasi-extinctions remain unlikely when fecundity is reduced by a factor up to F = 0.5. A decline of 50% in fecundity reduces the number of females in all subpopulations below one thousand animals (Appendix S1: Fig. S3, top panels). Similarly, epizootic outbreaks as single stressors result in temporary reductions in population sizes rapidly compensated by higher growth rates following mass mortality events. This is because we used a PDV specific mortality rate sparing young females (H€ ark€ onen et al. 2007), which does not result in increased quasiextinction probability (Appendix S1: Fig. S3, bottom panels). Nevertheless, it is important to note that the occurrence of an epizootic outbreak before complete recovery of the population from the previous event may reduce the population Fig. 3. Baseline scenario without environmental stressors. Each line represents a subpopulation growing to its carrying capacity. Lines represent average values across 100 replicates per scenario and therefore do not show the annual stochastic variation of 5% used in the model. Baseline parameter values: D = 0 (no delay in hunting quotas), Q = 0 (hunting quota is zero), F = 1 (maximum fecundity), and E = 0 (no epizootics). v www.esajournals.org size to quasi-extinction levels, in particular if mortality is not age-biased (not explored in this case).

Multiple-stressor scenarios
With the exception of moderate and high hunting quotas, which often led to quasi-extinctions of vulnerable subpopulations, none of the other stressors generated high risks of quasi-extinction to any of the subpopulations when introduced one by one. However, natural conditions are more complex, where populations are often affected by several types of stress simultaneously.
We created scenarios where multiple stressors are present and illustrated how hunting can affect the probability of quasi-extinction in different subpopulations. Under conditions where the annual epizootic probability is E = 0.07 and fecundity is reduced by 10% (F = 0.9) due to pollutants, many subpopulations face quasi-extinction under moderate hunting quotas (Q = 0.05). Introducing a delay of 10 yr in the calculation of the quota substantially increases the risk relative to a delay of 5 yr (Fig. 5). When hunting quotas are high (Q = 0.10), most subpopulations are likely to go quasi-extinct. When hunting quotas are low (Q = 0.01) or absent (Q = 0), none of the subpopulations face quasi-extinction. A further reduction in fecundity to 70% (F = 0.7) when hunting quotas are intermediate (Q = 0.05) will substantially increase the risk of quasi-extinction in many subpopulations (Appendix S1: Fig. S4).

Sensitivity analysis of quasi-extinction risk
To illustrate the sensitivity of quasi-extinction risk estimates to the full range of parameter values of multiple stressors, we chose three subpopulations based on their effective carrying capacities (Fig. 3): Anholt, the most vulnerable; Hesselø, with is intermediate; and Koster, the most resilient. We found that Anholt is at risk of quasi-extinction when the hunting quotas equal or exceed 0.025, but only when fertility is reduced and epizootic outbreaks are introduced (Fig. 6, first row). The risk of quasiextinction increases at Hesselø with hunting quotas at 0.05, especially when fecundity is reduced, even when epizootic outbreaks are rare (Fig. 6, second row). Koster is only likely to decline below the quasi-extinction threshold during the 100-yr simulation periods when hunting quotas are high (Fig. 6, third row). Our results show that the harbor seal population is resilient to multiple stressors during low hunting regimes (Fig 6, left column) while increasing hunting leads to extinction prone subpopulations especially under multiple stressors (Fig 6, right column).
We further explored the effect of hunting biased toward different age classes and found that a hunt directed toward pups poses a generally lower risk of quasi-extinction compared to a hunt directed toward older animals (Fig. 7). The subpopulation Anholt has high extinction risk also for a hunting quota of 0.025 and with an age-unbiased hunt (Fig. 6, upper right corner). Biased hunting toward older age classes gradually increases the extinction risk (Fig. 7), and the effects are magnified by introducing multiple stressors. An increase in the hunting quota to 0.05 makes quasi-extinction the most likely outcome even for subpopulations of intermediate size (e.g., Hesselø) when hunting is biased toward adults.

DISCUSSION
Harbor seals in Denmark and Sweden have recovered from low population sizes three times v www.esajournals.org during the last century and are now seemingly viable and abundant. However, signs of declining growth trends and the threats posed by xenobiotics and recurrent epizootics, combined with hunting and bycatches, warrant a careful risk assessment especially since demands for increased hunting are being raised. We constructed an age-structured metapopulation projection model with realistic carrying capacities and migrations rates, building on the extensive long-term studies conducted on this particular population. Quasi-extinction risks were estimated for a range of parameter values and single-and multiple-stressor scenarios. Moderate levels of single stressors (hunting, reduced fecundity, and epizootics) typically led to a substantial reduction in the long-term population size below the environmental carrying capacity, but the quasi-extinction risks of subpopulations seldom exceeded 1% (Table 2). However, hunting quotas above 5% of the metapopulation size led to high extinction risks for the smaller subpopulations and cannot be regarded as sustainable even in a single stressor scenario. Contrastingly, in a multistressor scenario, even the lower hunting quota of 2.5% and a slight reduction in fecundity (90% of baseline) together with an epizootic frequency of 7% reduced the metapopulation size down to two-thirds of carrying capacity and resulted in high quasi-extinction risks for Anholt and SW Kattegat, the most vulnerable subpopulations.
Epidemics of infectious diseases causing mass mortalities have been repeatedly observed in marine mammals, such as canine distemper Fig. 6. Quasi-extinction probability of the Anholt, Hesselø, and Koster subpopulations under a range of hunting quotas (Q), fecundity (F), and probabilities of epizootic outbreaks (E), keeping a 5-yr delay in the calculation of the hunting quotas. Anholt is the most vulnerable and Koster is the most resilient subpopulation. (A) Hunting is here assumed not to be age-biased. Note that hunting has the greatest impact on quasi-extinction risk. In the Anholt subpopulation, quasi-extinction is likely even under moderate levels of stress. In the Koster subpopulation, quasi-extinction is likely only under high levels of stress.
(CDV) among Baikal seals (Pusa sibirica) in Siberia in 1987-1988, when about 100,000 seals were found dead (Osterhaus et al. 1988, Grachev et al. 1989. In 2000, another epidemic occurred in the Caspian Sea killing more than 10,000 Caspian seals (P. caspica; Kennedy et al. 2000). The increase in disease outbreaks affecting aquatic mammals over the past decades (Gulland and Hall 2007) along with the increasing number of newly discovered strains of cetacean morbilliviruses (West et al. 2013, Groch et al. 2014, Stephens et al. 2014) has caused concern for marine mammal population health, demanding the inclusion of epizootics in models of population viability. One of the most well-documented epidemics in marine mammals occurred in the European harbor seal population when PDV spread across regions during two summers, in 1988summers, in and 2002summers, in (H€ ark€ onen et al. 2006summers, in , Stokholm et al. 2019). In the current study, we used an annual outbreak probability of 7% in several scenarios, which is equivalent to the observed time between the two PDV outbreaks (1 in 14 yr). Mortality rates and age structure of the dead seals of the modeled epizootic were also chosen to mimic PDV outbreaks (H€ ark€ onen et al. 2006(H€ ark€ onen et al. , Olsen et al. 2014). However, many other infectious diseases causing significant mortality might emerge (Siebert et al. 2007), and we elaborated a general scenario for repeated mass mortalities with frequencies of emergence varying from 1% to 10%. In order to evaluate threats from different pathogens, we did not include immunity of survivors to the disease. It is known that including immunity in modeling PDV epizootics in harbor seals causes a slight decline in extinction risk while new cohorts quickly dilute the proportion of immune survivors (Harding et al. , 2003(Harding et al. , 2005. We find that epizootics of the order of magnitude of PDV have the potential to dramatically reduce long-term carrying capacities in combination with low levels of other stressors (Table 2).
Polychlorinated substances such as PCBs, DDT, dioxins, and flame retardants have been shown to affect a number of endocrine processes in marine mammals. Effects were first noted in Baltic ringed seals, where fertility rates dropped to 15% in the 1970s (Helle 1980). A disease syndrome showed multiple disorders including occluded uterine horns (Bergman and Olsson 1985), osteoporosis and alveolar exostosis (Mortensen and Pedersen 1993). Furthermore, the immune function can also be affected , with observations showing a 25% reduction in the production of T cells in harbor seals feeding on fish with PCB levels mimicking Baltic conditions . Thus, seals exposed to these substances could be more vulnerable to chronic and epidemic diseases. These effects will result in reduced rate of population increase, either by reducing fecundity or increasing mortality rates or both , Sonne et al. 2020). In our model of harbor seals, pollutants such as PCBs were assumed to have a negative effect on fecundity. Decreased fecundity had severe effects on the metapopulation longterm abundance in combination with low levels of other stressors. A 30% reduction in fecundity reduced the long-term population size to Fig. 7. Age-biased hunting with a hunting quota of 0.025 in the Anholt subpopulation. Hunting directed toward adults can significantly increase the risk of quasi-extinction of vulnerable subpopulations relative to an age-unbiased hunting regime. Hunting directed toward pups had the weakest impact on quasi-extinction risk.
one-third of the carrying capacity (Table 2, Scenario 8).
Population abundance data based on surveys show large variation due to weather conditions and seal behavior and, consequently, it takes a minimum of 5 yr to identify a 5% change in abundance trends , Svensson et al. 2011. In order to mimic the uncertainty in population size estimates, we introduced a delay in calculating the hunting quota. We found the delay to have a small effect on long-term population extinction risks for most multistressor scenarios. However, an increase from a 5-yr delay to a 10-yr delay for a hunting regime of 2.5% slightly increases the quasi-extinction risk from 23% to 37% for the SW Kattegat population and from 97% to 98% for the Anholt population (Fig. 5). A higher hunting quota (5%) together with an increased delay in quota calculation (10 yr) led to a further increase in quasi-extinction risks from 0% to 9% for Lysekil, 1% to 20% for Marstrand, and 7% to 32% for Onsala. All these three subpopulations are intermediately resilient subpopulations. The age structure of the hunt is of critical importance, and an age structure skewed toward adult females led to a threefold increase in the extinction risk.

Multistressor scenarios
Previous studies of seal quasi-extinction risks have mostly focused on single stressors (Elmgren 1989, Harding et al. 2005, Olsen et al. 2014). However, natural populations are typically affected by multiple stressors (Munns et al. 2007, Chilvers 2012) that have to be taken into account to successfully design management efforts (Nacci et al. 2005, Munns 2006). Our multistressor analysis revealed that even if a seal population is healthy and close to carrying capacity, the combined effects of the several mild stressors can have severe consequences. Under the low hunting regimes (1-2.5% hunting quota), quasi-extinction risk is moderate (below 1%), while for higher hunting quotas of 5%, quasi-extinction risks quickly reach unsustainable levels (Mace and Lande 1991). Additionally, our study shows that a slight decrease in fecundity, together with realistic scenarios of random epizootic events, can lead to drastic population declines within the next 100 yr. Although concentrations of pollutants such as PCBs tend to decrease in the North Atlantic (Aguilar et al. 2002), new xenobiotics can contribute to combined effects under the multistressor scenario (Bjurlid et al. 2018).
In contrast to most other studies analyzing the entire metapopulation as a single population, we investigated the Swedish-Danish subpopulations of harbors seals as a metapopulation. The application of multiple-population viability analysis to assess the extinction risk of metapopulations is a useful tool that can guide management and conservation strategies (Neville et al. 2020).

Management implications
Hunting quotas should be coordinated among countries according to the 2006 HELCOM seal recommendation and the EU Habitats directive, where the main goals for the management of seal populations are defined as natural distribution, natural abundance, and a health status that ensures the persistence of the populations in the future. Consequently, quotas for hunting must, by international agreements and laws, ensure that these main goals are not jeopardized.
A hunting quota of 1% is equivalent to approximately 250 animals in the study area, and a 2.5% quota corresponds to 625 animals, based on the estimated metapopulation size at present. Currently, the annual hunting quota for harbor seals is 440 individuals in Sweden (Swedish EPA), although in 2019 only 292 were actually killed, and in Denmark about 10 are killed annually (Olsen et al. 2018), implying that current levels of hunting can already affect local carrying capacities, and that any additional future stress may cause severe population declines. Having experienced two outbreaks of PDV, influenza A and other infectious diseases (Olsen et al. 2010) over the past 40 yr, future epidemics are impending. Pollutants at only slightly higher levels than recorded at present can have detrimental effects on fecundity and immune system responses , and cumulative effects of known and unknown xenobiotics pose further concern (Sonne et al. 2020). Additional stressors not explored in the current study include disturbance at breeding sites caused by anthropogenic activities (Montgomery et al. 2007, poaching and bycatches. Furthermore, inaccuracies in the reporting of the hunt can cause errors in the calculation of the quotas and, therefore, it is also important to document the v www.esajournals.org age and sex of the hunted individuals, since a bias toward adult females has a threefold effect compared to hunting pups.

CONCLUSIONS
Our study illustrates how marine mammal populations can be abundant and close to carrying capacity, but nevertheless vulnerable to just slightly increased levels of anthropogenic stressors. The history of overexploitation of marine mammals worldwide provides the hard evidence for this. Currently, when populations of several marine mammals have recovered after decades on the red list, it is important to not repeat the historical mistakes. Our main conclusion is that abundant marine mammal populations, such as the harbor seals studied here, also require continued careful monitoring of population trends in abundance, reproductive health, and documentation of new stressors. In particular, increased exploitation from hunting must be carefully evaluated by population viability analyses before implementation.

ACKNOWLEDGMENTS
The study was financed by Viltforskningsanslaget, Swedish Environmental Protection Agency, and by the BONUS program BaltHealth which has received funding from BONUS (Art. 185), funded jointly by the EU, Innovation Fund Denmark (grants 6180-00001B and 6180-00002B), Forschungszentrum J€ ulich GmbH, German Federal Ministry of Education and Research (grant FKZ 03F0767A), Academy of Finland (grant 311966), and the Swedish Foundation for Strategic Environmental Research (MISTRA). The authors declare that they have no conflict of interest.