Ninety years of change, from commercial extinction to recovery, range expansion and decline for Antarctic fur seals at South Georgia

With environmental change, understanding how species recover from overharvesting and maintain viable populations is central to ecosystem restoration. Here, we reconstruct 90 years of recovery trajectory of the Antarctic fur seal at South Georgia (S.W. Atlantic), a key indicator species in the krill‐based food webs of the Southern Ocean. After being harvested to commercial extinction by 1907, this population rebounded and now constitutes the most abundant otariid in the World. However, its status remains uncertain due to insufficient and conflicting data, and anthropogenic pressures affecting Antarctic krill, an essential staple for millions of fur seals and other predators. Using integrated population models, we estimated simultaneously the long‐term abundance for Bird Island, northwest South Georgia, epicentre of recovery of the species after sealing, and population adjustments for survey counts with spatiotemporal applicability. Applied to the latest comprehensive survey data, we estimated the population at South Georgia in 2007–2009 as 3,510,283 fur seals [95% CI: 3,140,548–3,919,604] (ca. 98% of global population), after 40 years of maximum growth and range expansion owing to an abundant krill supply. At Bird Island, after 50 years of exponential growth followed by 25 years of slow stable growth, the population collapsed in 2009 and has thereafter declined by −7.2% [−5.2, −9.1] per annum, to levels of the 1970s. For the instrumental record, this trajectory correlates with a time‐varying relationship between coupled climate and sea surface temperature cycles associated with low regional krill availability, although the effects of increasing krill extraction by commercial fishing and natural competitors remain uncertain. Since 2015, fur seal longevity and recruitment have dropped, sexual maturation has retarded, and population growth is expected to remain mostly negative and highly variable. Our analysis documents the rise and fall of a key Southern Ocean predator over a century of profound environmental and ecosystem change.


| INTRODUC TI ON
The unregulated harvesting of marine predators causes defaunation and trophic downgrading with negative repercussions for ecosystem processes (Dulvy et al., 2021;Estes et al., 2011;Luypaert et al., 2020;McCauley et al., 2015).When harvesting ceases, the targeted species are often heavily depleted and their path to recovery is thereafter determined by novel ecosystem states, predator-prey interactions, harvesting regimes and environmental stressors (Hobbs et al., 2009;Jachowski et al., 2015;McNellie et al., 2020;Stier et al., 2016).Because these effects create transient population states (Hastings, 2004), only long-term analyses can identify which novel conditions enable recovery (Hastings et al., 2018), and inform best management practices (e.g.Larkin, 1996;Lindegren et al., 2009).However, as the supporting population data are often unavailable, the status of many affected species remains uncertain.
After near extirpation by 1907and rediscovery in 1933-1936(Mackintosh, 1967)), the total population grew to an estimated 369,000 animals in 1976 (Payne, 1979), and from there to 1.55 million animals in 1991 (Boyd, 1993).This represents 95%-98% of the species' global abundance (Boyd, 1993;Forcada & Staniland, 2018), making the Antarctic fur seal at South Georgia the second most abundant pinniped in the Southern Hemisphere and the most abundant eared seal in the World (Costa et al., 2006).This abundant predator is considered a prime indicator species for resource management in the Southern Ocean ecosystem (Agnew, 1997;Croxall & Prince, 1979).
However, the status of this population remains uncertain because the latest abundance estimate did 'not reflect the pup production had 1990/91 been an average year' (Boyd, 1993, p. 22).Furthermore, the adjustment made to estimate the total population size was also derived from unrepresentative data, because 'nothing is [was] known of the structure of the male herd other than the mean age of a small sample of breeding bulls' (Payne, 1979, p. 9).Subsequently, and without any additional population surveys and supporting data, South Georgia was reported to carry between 4.5 and 6.2 million fur seals in 1999/2000.This is of concern because accurate population estimates are essential to evaluate the conservation status of a species and to grant adequate levels of protection.
Without a real and credible baseline, the notion of the Antarctic fur seal as super-abundant species contributed to (i) unwarrantedly depicting it as an ever-increasing environmental nuisance (Bonner, 1985;Convey & Hughes, 2022;Foley & Lynch, 2019), despite this being a naturally abundant species (Cook, 1777, v. 2, p. 213;Weddell, 1825, p. 54); and (ii) delisting it as a 'Specially Protected Species' in 2006 by the Antarctic Treaty Consultative Parties (SCAR, 2006), without a specific action plan for monitoring and reassessment (Jabour, 2008).Almost coincident with this decision, the population of Bird Island, the likely epicentre of historical recovery at South Georgia, collapsed in 2009 and has continued to decline ever since (Forcada & Hoffman, 2014).Antarctic fur seals have increasingly become more sensitive to external stressors (Forcada et al., 2008;Forcada & Hoffman, 2014;Krause et al., 2022), and their numbers are also declining elsewhere in the Southwest Atlantic sector of the Southern Ocean (Hofmeyr, 2016;Krause et al., 2022), where most of the global population currently resides.
To produce an improved population trajectory and to identify the drivers of population size changes, we develop new methods for longterm population assessment.Conventional methods to estimate pinniped abundance are based on counts of animals at breeding locations on land or ice taken at peak breeding dates (Buckland & York, 2018;Hammond et al., 2021).These counts are subsequently adjusted for detectability bias, which occurs when animals are within detection range but cannot be perceived by observers; and when they are unavailable while within or outside detection range.This availability bias is typically modelled from haul out studies (Southwell et al., 2008;Thompson et al., 1997;Womble et al., 2020), which may not account for the absences of individuals due to delayed accession to breeding, and sexual segregation in haul out behaviour due to variable rates of fecundity, breeding propensity, immigration and temporary emigration.This leads to additional adjustments typically based on life tables (e.g.Hammill et al., 2013;LaRue et al., 2021;Pitcher et al., 2007).However, these require demographic parameters that are often unavailable (e.g.Payne, 1979), which may not be specific to the focal population and may not account for spatiotemporal variation as required in serial abundance estimates.
Instead, we use integrated population models (IPMs), which jointly analyse individual and population-level data in a single statistical model, allowing the use of available information, and to achieve comprehensive error propagation throughout the model (Besbeas et al., 2002;Newman et al., 2014;Schaub & Abadi, 2011;Schaub & Kéry, 2021).IPMs allow for imperfect detection and observation error and carry this uncertainty over datasets, producing better and often more precise inference than other methods.They effectively combine in situ population counts and demographic data to directly infer demographic rates and availability bias.Here, we jointly analyse age-at-count and reproduction data, as well as capture-mark-recapture (CMR) data from the tagging and genetic sampling of individuals to obtain key demographic rates.Using flexible CMR model structures, we account for observed (counted) and unobserved states, estimate recruitment and breeding probabilities in both females and males, and temporary emigration in territorial males, allowing us to infer total population numbers while deriving adjustments for survey counts.Moreover, we produce the first robust life history and demographic assessment for the male of this species, which is the key to accurately evaluating the total Antarctic fur seal population size.
We apply these methods to photographic counts from the latest dedicated aerial surveys of South Georgia (2006/07 and 2008/09), which we combine with genetic, demographic, reproductive and count data from an extensively studied representative population at Bird Island (Doidge et al., 1984;Forcada & Hoffman, 2014).We thus provide an assessment of total population and range expansion at South Georgia, and a population trajectory spanning the 90 years since commercial extinction, subsequent increase, range expansion and local decline at Bird Island, in relation with external drivers.

| South Georgia survey data
Two helicopter surveys of the entire South Georgia archipelago, conducted from 11 to 14 December 2006 and from 3 to 9 December 2008, provided aerial photographs from digital SLR cameras of the seal breeding range near the peak pupping days.Counts were extracted from images by at least two observers working independently using consistent training and methodology, and according to four categories: pups, females, territorial males and non-territorial males (Appendix S1).The pupping season (mid-November to mid-January) typically peaks between 6 and 10 December with high spatial synchrony and around the peak pupping day, numbers of newborn pups, breeding females and territorial males at each location reach their maximum values and the counts are less variable than at later dates (Duck, 1990;Forcada et al., 2005;Hofmeyr et al., 2007).
Antarctic fur seals at Bird Island (54°00′ S, 38°02′ W), South Georgia, were surveyed through complete island counts on seven more occasions.Visual counts were conducted from elevated vantage points by trained fur seal biologists in December 1988, 1989, 1990, 1994and 1997. In December 2003and 2016, additional
Genetic samples were collected during daily surveys of nearly all the territorial individuals observed at SSB since 1995, and between-year genetic recaptures are available for 13 seasons (1995-2007 inclusive).
For each sample, total genomic DNA was extracted using a modified phenol-chloroform protocol and genotyped at nine highly polymorphic microsatellite loci as described by Hoffman et al. (2003).None of the loci showed significant deviations from Hardy-Weinberg or linkage equilibrium and the genotyping error rate was very low at <0.005 per reaction.The probability of identity was also very low (1.354 × 10 −12 ), indicating that identical genotypes most likely represent resampled individuals.To identify genetic recaptures, the microsatellite genotypes were checked for duplicate entries using the program identity (Allen et al., 1995) as described in detail by Hoffman et al. (2006).

| Females
We used a different IPM for each sex with an age-stage-structured population projection model (Caswell, 2001) as main building block, assuming a birth pulse and a pre-breeding census.The female model describes survival and accession to breeding (recruitment) for immature females (pre-breeders; P) of ages 3-7, and the transition of individuals through breeding states breeder (B) and mature non-breeder (N).The expected number of females in year t + 1 is where N j,t is the number of females for ages j = 1, …, 7, and N N,t and N B,t are number of mature non-breeders and breeders respectively at year t.Apparent survival parameters are: 0,t for pups from birth to tagging at 1.5 months of age, 1,t from tagging to the end of the first year (census year), j,t for ages j = 2, …, 7, and N,t and B,t for mature non-breeders and breeders, between years t and t + 1. Transition parameters between t and t + 1 are a j,t , defined as the probability that an immature female (pre-breeder) of age j starts to breed at that age (Pradel & Lebreton, 1999); parameters B−,t and N−,t are the probabilities of being in states B or N in year t and in year B or N at t + 1, conditional on surviving between years.F is the proportion of females at birth, and is the immigration rate, defined as number of breeding immigrants in year t + 1 per mature female in the population in year t (Abadi, Gimenez, Ullrich, et al., 2010).
Here, CMR data contribute to estimate survival, recruitment and breeding probability.A multistate Cormack-Jolly-Seber markrecapture model (MSCJS) was formulated as a Bayesian state-space model (Gimenez et al., 2007), including 11 states: P j , for j = 0-7, where P 0 is pups, states P 1 to P 7+ are pre-breeders of ages 1-7 and above, B is breeders, N is mature non-breeders or skippers and D is an absorbing state for dead or permanently emigrated females.Prebreeders at ages j = 3-7 access breeding (state B) with probability a j , or stay as pre-breeders with probability 1 − a j .Some pre-breeders, typically of ages seven 7 or more, are never observed breeding., seen as breeder (3), with probability p B t , seen as mature non-breeder (4) and with probability p N t , and not seen (5).This is formulated as where rows are states P 0 to P 7+ , B, N and Dead, and columns are observation types (1-5) at year t.We denote the CMR likelihood as  MSCJS (h| , a, , p), where h is the female recapture history as defined by observation types.Further details on Bayesian multistate model implementations can be found in Gimenez et al. (2007) and Kéry and Schaub (2012).
Data on pup production, equivalent to the breeding female population, contribute to estimate population sizes and demographic rates, including immigration.The state-space formulation is thus further composed of a state process defined by the matrix population model, with likelihood  sys (N| , a, , ), and uses Poisson distributions to introduce demographic stochasticity as follows.and where N B includes new recruits and returning mature females as breeders, and N BI includes the same demographic groups plus immigrants.The distinction between N B and N BI allows for estimates of number of immigrants N I , and the fecundity (or pregnancy) rate, F t , which is the number of mature females (N B + N N ) in year t which survive and breed in year t + 1 (excluding recruits and immigrants in year t + 1), divided by the number of mature females surviving to year t + 1.
The observation process conditional on this state process is based only on the count of breeding animals (n BI ), equivalent to the sum of daily pups produced in each season, y 2,t , which is observed with little error (Forcada et al., 2005).We model these counts as y 2,t ∼ Poisson N BI,t , where N BI,t is the number of breeders including immigrants, with likelihood  obs y 2 | N B .The product of state and observed process likelihoods is the likelihood of the state-space To evaluate the proportion of the estimated study population, , that represents a female count at the SSB (denoted with superscript S) obtained between 1 November and 31 December (e.g.n s j for day j), we model the observed daily total female counts with a generalised additive mixed model (GAMM).We select a Poisson-log-normal structure (Gelman & Hill, 2006), because the daily count series data tend to be overdispersed, with where d j is the jth day of the survey season, b k is the basis of a penalized regression spline and j are temporal random effects, with j ∼ N 0, 2 .The hyperparameter measures overdispersion, with = 0 corresponding to a simple Poisson generalised additive model.

Alternatively, we consider
where is an autocorrelation parameter with values between −1 and 1, and j,t−1 is the residual associated with observation at day t − 1.
The likelihood is  c (n| , b, , , ), where are unknown parameters (penalties), b is the vector of k basis functions and λ are smoothing parameters.
With the daily female count model, the likelihood of the IPM is the joint likelihood of the different parts; The combined likelihood produces two derived quantities, where M is for mature females, B for breeders and E n s j is the predicted count at SSB on the jth day.These are used as multipliers of observed counts at other locations to obtain estimates of total number of mature , where n c j is the female count at colony c on day j, and the second quantity is equivalent to pup production.These equations assume that the dynamics of daily counts at SSB and other locations are highly correlated, which is accurate within Bird Island and across other locations (see Section 3).
To account for pre-breeders and obtain an estimate of total number of females at South Georgia, we use where is the estimated number of pre-breeders of age k, for k = 1 to 7+, at SSB.

| Males
The male IPM has a combined age-stage matrix model with preterritorial males (P) of ages 1-11+ (ages 11 or older), and mature states including territorial (T) and mature non-territorial (N).The expected number of males in year t + 1 is where N j,t is the number of males for ages j = 1-11, and N N,t and N T,t are numbers of mature non-territorial and territorial in year t.
Apparent survival parameters are 0,t for male pups from birth to tagging, and between year t and year t + 1 are 1,t for weanlings, from tagging to the end of the first year; j,t for males of ages j = 2-11; and N,t and T,t for mature males.Transition parameters between years t and t + 1 are a j,t , defined as the probability that an immature male of age j starts to hold territory at that age; and T−,t and N−,t are the probabilities of being in states T or N in year t and in T or N at t + 1, conditional on surviving between years.The proportion of males at birth is M , and is the proportion of immigrants among territorial males.
CMR data contribute to estimate apparent survival, first-time accession to territory and probability of being territorial having held a territory before.Pup production data contribute to estimate population sizes and all demographic rates.For implementation, the CMR analysis was split into separate likelihoods given the different type, sample size and short temporal overlap between datasets.Thus, pups with implanted PIT-tags at birth (since 2001) represent the pre-territorial component, and are used to model j,t and a j,t ; while genetic samples from territorial males (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) contribute to estimate their survival, territory tenure and temporary emigration rates.Because the elapsed time between birth and accession to territory can be 10 years or more, observations of pre-breeders were limited and sparse (over 22 years, only 4.6% of 1773 PIT-tagged male pups were observed as territorial males), making IPM fitting and selection inefficient.As a valid alternative (Besbeas et al., 2002), we therefore use estimates from a simplified MSCJS recruitment model (Pradel & Lebreton, 1999) of the PIT-tag data as informative priors for j,t and a j,t in the IPMs.
Nearly all the territorial males are genetically sampled every season, allowing for robust estimates of abundance using CMR models.
We use a multistate Jolly-Seber model, selecting a hierarchical statespace formulation (Kéry & Schaub, 2012) with parameter-expanded data augmentation (Royle & Dorazio, 2009).Here, the number of individuals alive during the study is N s , and b t is the probability that a new territorial male enters the population in year t, so that the number of new territorial males in a year is B t = N s b t .The augmented dataset contains M individuals, of which N s is a subset, allowing for derivation of b t in a restricted dynamic occupancy model framework (Royle & Dorazio, 2009).This relies on t , defined as the probability that a male from M holds territory for the first time in year t, so that The state process for sampled male i and year t, is then formulated as where the states in rows and columns are P (not yet territorial in the population), T, N and Dead (or permanently emigrated), between years t and t + 1.The corresponding observation process is where p T t is the recapture probability of territorial males.As nonterritorials went unobserved, the multistate formulation part is analogous to a temporary emigration model (Schaub et al., 2004).We denote the likelihood as  MSJS ( , , , p , N|h).Because the number of new territorial males each year (B t ) confounds locally born males starting to hold territory with immigrants that never held a territory at SSB before, these two components were separated with the demographic part of the IPM and the combined datasets (see below).
Data on pup production contribute to estimate population sizes of pre-territorial males and all demographic rates, including immigration.The state-space formulation is thus further composed of a state or system process, defined by the matrix population model with likelihood  sys (N| , a, , ) and uses Poisson distributions to introduce demographic stochasticity: The observation process conditional on this state process is based on the count of male pups produced each year, N 0,t , and demographic rates from the PIT-tag recaptures for the pre-territorial component.Estimates of N N,t+1 , N T,t+1 and B t are directly obtained from the MSJS model component, which allows for immigration, and an estimate of immigrants is obtained as The proportion of the estimated study population, , that represents a territorial male count at SSB, obtained between 1 November and 31 December, is obtained from daily total male counts n s j and a Poisson GAMM (see female IPM).The combined IPM likelihood produces the derived quantities N S M,t ∕ E n s j and N S T,t ∕ E n s j , where M and T are for mature and territorial males respectively, and E n s j is the predicted count at SSB on the jth day.Used as multipliers of observed counts at other locations, we obtain estimates of total number of mature males and territorial males , where n c j is the male count at colony c on day j.Adjusting for pre-territorial males, an estimate of total number of males at South Georgia is where NS is the estimated number of pre-territorials of age k, for k = 1-11+, at SSB.

| Model implementation and selection
The total likelihood in IPMs is the product of likelihoods for the different datasets, which are assumed to be independent (Besbeas et al., 2002), although our counts include seals that were also marked.
As shown elsewhere, violation of the independence assumption has little impact on model estimation and performance (Abadi, Gimenez, Arlettaz, et al., 2010;Weegman et al., 2021), and thus a significant bias in our results is highly unlikely.
We use Markov Chain Montel Carlo methods in BUGS language and the program JAGS (Plummer, 2003), run from R (v4.2.2; R Core Team, 2022) with the packages jagsUI (Kellner, 2021) and rjags (Plummer, 2022) to fit the IPMs.To specify priors, we use normal distributions with numbers of pre-breeders, breeders and mature non-breeders observed in the first year for the initial population sizes.
For the capture-recapture part, model fit was initially assessed with frequentist multistate goodness-of-fit testing (Pradel et al., 2003) using R2ucare (Gimenez et al., 2018).Selection of best models was subsequently based on WAIC (Gelman et al., 2014;Watanabe, 2010), and the IPMs were fitted with the best mark-recapture modelling options.We selected temporal variation in parameters to better reflect the population state at the time of the surveys.For the female model, we used independent temporal random effects for j and a j for ages 1-7 and 3-7 respectively, ψ NB , log(ι), p P1 , p P2 , p B and p N (e.g. , where NB i,t ∼ Normal 0, 2 NB ), with diffuse Uniform priors for 2 * .We used correlated temporal random effects for B , N and, BN , where i,t ∼ MVN 0, Σ ,t (multivariate normal distribution), where Σ ,t is the temporal variance-covariance matrix describing the temporal variance of each parameter ( ) and the temporal covariance among multiple parameters (Link & Barker, 2005).
For this matrix, we chose priors from an inverse Wishart distribution (Gelman & Hill, 2006) with scale matrix Ω as the identity matrix and four degrees of freedom.
In the male IPM, we chose the best combinations of , and with temporal random effects or constant values over time depending on WAIC-based model selection.As the available data for the Jolly-Seber likelihood did not support estimation with complex model structures, to test full time variation in T , N , TN and NT , we used correlated temporal random effects.For juveniles and pre-territorials, j and a j were selected as informative Beta priors with shape parameters obtained from an independent best fit of a multistate mark-recapture model of the PIT-tag data set.These priors incorporate temporal variation as observed in the data record, which reflects uncertainty due to environmental stochasticity at the time of the survey counts.
In both IPMs, the GAMMs of daily counts were implemented using the methods of Wood (2016), using multivariate normal priors on the smooth coefficients in the JAGS models; and the smoothing parameters and smoothing penalty matrices directly specified the prior multivariate normal precision matrix.We implemented preliminary model selection based on WAIC to evaluate the need to incorporate overdispersion and autocorrelation before incorporating the Poisson models in the IPMs.
For inference, we produced 200,000 iterations of three Markov chains using dispersed parameter values as starting values and discarded the first 50,000 samples of each chain as burn-in, thinning the remainder to every 10th sample.We assessed chain convergence visually using trace plots, through the mixing of the chains and sample autocorrelation plots, and using the R (potential scale reduction factor) statistic (Gelman & Shirley, 2011).
We used posterior simulations for inference on generation time and longevity, based on the core matrix model of the female IPM.
Longevity was the age at which survivorship for a cohort fell below a critical value; and the generation time was the time required for a population to increase by a factor equivalent to the net reproductive rate (R0; average number of offspring produced by an individual in its lifetime taking mortality into account ;Caswell, 2001).

| Range expansion, trends and trajectory
We produced a total abundance estimate for South Georgia from the latest surveys, and evaluated range expansion comparing the distribution of breeding females (and total pup production) for three surveys conducted in 1971-1975, 1991 and 2007-2009 respectively.
To obtain a population trajectory at Bird Island, we combined all of the pup production data available since the cessation of sealing.This includes data from 1933, 1936 and 1950 with estimated pup counts of 29, 56 and 1140 respectively (Bonner, 1968); field counts from 1957/58 to 1963/64 (Bonner, 1968); corrected field counts from 1971/72 to 1973/74 (Payne, 1977) and from 1975/76 to 1976/77 (Croxall & Prince, 1979); and eight estimates of pup productivity for 1989-1991, 1995, 1998, 2004, 2009 and 2017 based on comprehensive whole island counts.We produced estimates of mature female abundance from 1989 to 2017, and productivity from 1933 to 2017, extending the SSB capture-recapture analysis back to the 1985 season and using a simplified IPM for mature females (Appendix S2).
We identified phases of change in external pressures on the fur seal population including climate, the physical and biological environment, and commercial extraction of Antarctic krill (Euphausia superba).For climate effects, we used the Southern Annular Mode index (Marshall, 2003), available at http://www.nerc-bas.ac.uk/ icd/gjma/sam.html,which is the leading mode of variability in the Southern Hemisphere atmospheric circulation, and contributes to sea surface temperature (SST) variability at South Georgia (Meredith et al., 2008).We used records of satellite SST Daily Optimum Interpolation (OI), AVHRR Only, Version 2.1 (Huang et al., 2020), averaged over a grid of 0.25° latitude (−54 to −53.5°) and longitude (−39 to −38°) resolution, from 1981 to the current day.We obtained commercial krill extraction for the CCAMLR management Area 48.3, which includes South Georgia, from 1977 to 2020 via https://www.ccamlr.org/en/fisheries/krill -fishe ries.And we used seasonal mean krill size in the fur seal diet at Bird Island as indicator of krill availability (Reid & Arnould, 1996), from 1990 to date; a high mean size in a given year reflects a reduced recruitment of juvenile krill and correlates with low krill biomass, with repercussions in the current and following fur seal breeding seasons (Forcada et al., 2008;Reid et al., 1999).Even though krill density estimates obtained with fishery-independent acoustic surveys are available (Fielding et al., 2022), the temporal resolution, interannual sampling consistency and the extreme variation in the annual estimates were deemed inadequate for this analysis due to a lack of precision and statistical power in any analysis based on this covariate.
We used time-frequency maps of coherence (correlation) between Southern Annular Model (SAM), SST, krill size and estimated log-population growth rate [ln(λ)] using wavelet transform multivariate analyses to investigate levels of synchrony in cycles among variables based on their wavelet modulus ratio (Keitt, 2008).This uncovered the non-stationary and temporal relationship between associated fluctuations in climate, physical and biological environments, and fur seal prey and fur seal population dynamics.

| Animal ethics
All animal handling procedures were evaluated and approved by the Animal Welfare and Ethics Review Body at the British Antarctic Survey, and this research was conducted under annual Regulated Activity Permits issued by the Government of South Georgia and the South Sandwich Islands (GSGSSI).Samples originating from South Georgia Islands were imported under permits from the Convention on International Trade in Endangered Species, and exported under permits issued by the GSGSSI and the UK Department for Environment, Food and Rural Affairs, under European Communities Act 1972.

| RE SULTS
The assessment of IPMs, including mark-recapture and GAMM model choices is detailed in Appendix S3.The following results describe estimates of vital and life history rates, demographic parameters, population structure and population size changes at SSB, Bird Island and across South Georgia.

| Female vital rates
The mean recapture probability was 0.128 [95% credible interval: 0.112, 0.147] during the first 3 years of life p P 1 when most individuals were unobserved.During the next 4 years, the mean recapture probability increased to 0.242 [0.242, 0.341] in pre-breeders p P 2 .For breeders p B it was 0.878 [0.856, 0.898] and it was best modelled as constant over time, like p P 1 and p P 2 .Recaptures of mature non-breeders p N were best modelled with temporal random effects and a population mean of 0.791 [0.671, 0.925].
Mean apparent survival rates were 0.35 [0.03, 0.82] for preweaned and weaned pups (first year females, P 1 ); they were high for pre-breeders aged 1 and 2, and were highly variable but lower for pre-breeders of ages 3-8 and above (Figure 1a).and for 2016-2022, shortly after the most severe El Niño on record, it was reduced to 0.150 [0.114, 0.193] (Figure 2b).
For females, the probability of first pupping at age (α 3-7 ) increased from 0.06 [0.01, 0.13] at age 3 to 0.43 [0.38, 0.53] at age 5, and decreased subsequently (Figure 1c).Estimates during the F I G U R E 1 Medians and 95% credible intervals of (a) female apparent survival for each breeding stage, (b) male apparent survival, (c) recruitment and breeding probability in females and (d) male accession to territory tenure for pre-territorials and experienced territorials.Black points and solid black bars are temporal means, and green and red points with dashed bars are estimates for the South Georgia surveys in 2007 and 2009.(Figure 2c).The observed mean age of breeding females declined significantly until 2015, and then increased during a subsequent period of persistently high juvenile mortality that started after 2015 as the surviving population aged.Similarly, the mean observed age at maturation increased from 4.4 (SE = 0.2) to 5.3 (SE = 0.1), contributing to a higher mean age of breeders, and therefore increasing the estimated mean generation time, from 10. 6 [10.1-11.3] to 11.9 [11-13] (Figure 2c).

| Male vital rates
The recapture probability of pre-territorial males p P , which was best modelled as constant over time and across ages 1-11, was 0.120 [0.101, 0.141].The recapture probability of territorial males p T , also best modelled as constant, was 0.925 [0.847, 0.999].Apparent survival rates for pre-territorial males were higher at ages 2-4 and declined in older individuals.After age 8, when most pre-territorials started to be competitive, survival estimates were highly variable, with a similar mean to territorial males (Figure 1b).For these animals, the temporal mean was 0.651 [0.533, 0.787], and the mean estimate prior to the surveys was slightly lower (0.588 [0.520, 0.679]).For temporary emigrants, survival was lower but highly variable (0.534 [0.090, 0.952]).
The probability of defending a territory for the first time (α  In contrast, the interannual variation of commercial krill extraction was uncorrelated with ln(λ), and the annual mean extraction represented a small percentage of the extant krill stock (see Section 4).

| Daily counts and adjustment factors
In 2007, the peak of the pupping season was on 8 December and the survey of South Georgia was conducted from 11 to 14 December.In 2009, the peak day, 6 December, was the central day of the survey, from 3 to 9 December (Figure 3).The best GAMM fit for daily numbers of seals at SSB in 2007 and 2009 indicated a higher overdispersion in female counts ( F,07 = 0.141 0.106, 0.184 and F,09 = 0.101 0.069, 0.140 ), than in territorial male counts ( M,07 = 0.022 0.001, 0.061 and M,09 = 0.026 0.001, 0.069 ) (Figure 3).The results of our analyses of 8 years with complete surveys of Bird Island are shown in Appendix S4 (Table S1), and the best models for the corresponding SSB female counts are shown in Appendix S4 (Figure S1).Four of eight surveys were in years affected by environmental anomalies and low pup production, which provided a good representation of long-term variation in peak female counts.Because most of the surveys were conducted when female numbers were close to the highest daily SSB counts (Fig- ure S1), variation in mean counts over the survey dates was relatively low (%CV = 1.6%-12%), and was well reflected in the GAMM predictions E n s j,t .
The Pearson's product moment correlation coefficient for paired whole Bird Island counts and predicted SSB counts (E n s j,t ; Table S1) was 0.992 [95% CI: 0.959, 0.999] and the correlation between the Bird Island counts and mean SSB counts was 0.982   and declined by 42. 1% [27.0, 54.5] from 2010 to 2017, which represents a mean annual rate of decline of 6.7% [3.9, 9.4] (Figure 4).
High levels of variation between 1989 and 2009 precluded the effective detection of significant trends in pup production, but after 2010 it fell by 44.8% [34.8, 53.4] over 8 years.represents a mean annual rate of decline of 7. 2% [5.2, 9.1].Given the high synchrony between SSB and Bird Island, we expect this trend to continue after 2017.Extending the time series of available productivity estimates back to the first half of the 20th century, Figure 4b shows a rapid exponential increase until 1975 (Payne, 1977).The mean productivity reached a maximum in the early 1980s, which was followed by a period of relative stability, and then the most recent drop after the abrupt decline of 2009 to values comparable to or below those of the early 1970s.
The geographical distribution of breeding females is shown in Figure 5c, along with the distribution range from the nearly complete surveys in 1971-1975and in 1991(Boyd, 1993;;Payne, 1977) to illustrate range expansion.The fur seal population of Bird Island comprised approximately 37% of the South Georgia total in 1975, but thereafter declined to 17.1% in 1991 and 5. 2% [4.7, 5.6] in 2009.
The breeding distribution in 2007-2009 extends eastwards along most of the suitable breeding habitat on the north coast, excluding unsuitable glacier and fjord ends, mostly at the northeast end.In the south, the breeding range only expands along the southwest.

| DISCUSS ION
This is the first assessment of the Antarctic fur seal population at South Georgia over a century of global change.We document its recovery from near extinction to become the most abundant otariid in the world, followed by a local decline.Such information provides the basis for an improved assessment of the conservation of the species worldwide, which is timely given the profound ecosystem changes occurring in the SW Atlantic region of the Southern Ocean, where the species is currently declining.
Using CMR and reproduction data from a representative study population, and visual and photographic count data on various geographical scales, we developed IPMs to address availability bias, one of the most common problems of evaluating marine predator populations.Methods to account for availability bias that rely on life tables have limitations when demographic data are inaccurate and unrepresentative (Hammond et al., 2021).To address this, Bayesian state-space population models can provide a flexible framework for data integration including detailed demographic parameters and other effects (e.g.Thomas et al., 2019).They can be extended to evaluate population viability and conservation management strategies (Saunders et al., 2018), assess population collapses (Maunder, 2004) and gain detailed insights into the drivers of population dynamics with limited data on marked individuals and study periods (Nater et al., 2021).We extended this framework to estimate availability bias, based on the proportion of the abundance of a representative population (SSB, Bird Island) that constituted a mature individual count, obtained in a given survey area, day and year.We incorporated uncertainty from different data sources, supported repeatability over past survey years and accounted for temporal specificity (Appendix S5), without making additional assumptions other than spatial synchrony in counts.
Spatial synchrony is often assumed, but untested, when applying abundance adjustments from data obtained at a single location to other locations (e.g.LaRue et al., 2021;Russell et al., 2019).In contrast, our results support the hypothesis of maximum temporal synchrony in numbers of breeders among locations, the high temporal correlation between peak female counts and total productivity, and the assumption that SSB count-based population adjustments have a wider spatial application, including at the opposite geographical end of the breeding distribution within South Georgia.
The IPMs considered separate demography and population dynamics by sex, given the large differences in life history between females and males.Sex structure in IPMs is preferable, although uncommon (but see Cleasby et al., 2017;Gamelon et al., 2021;Péron & Koons, 2012;Plard et al., 2019), and in our case it was justified by between-and within-sex differences in data types and data record periods, in a species with a high degree of sexual dimorphism.
This allowed for the most comprehensive demographic analysis for the male of the species to date, including first estimates of survival, recruitment, breeding propensity and temporary emigration; which highlighted major differences in population components by sex.
The estimated percentage of females in the breeding population was 68, whereas in males it was only 18, and the majority of these held active territories.This difference is typically observed in highly polygynous pinnipeds (e.g.Condit et al., 2014;Oosthuizen et al., 2019), and reflects strong survival-reproduction trade-offs.
Once becoming territorial at a mean age of 10 years, survival was limited to 3 years on average, reflecting the high survival cost of territoriality.If a territorial male emigrated, the probability to return as a territory holder was very low.As a result, up to 80% of males were pre-territorial, whereas the proportion of female pre-breeders was much lower, at around 32%.This indicates that most males were unobserved during counts of breeding beaches, and were hence unavailable to the observers during the survey period (see also Boyd et al., 1998;Staniland & Robinson, 2008).This is commonly found in pinniped surveys (e.g.LaRue et al., 2021;Pitcher et al., 2007), but requires representative and accurate adjustment estimates from non-count data.
The female IPM showed changes in vital rates and associated changes in female population structure, which since 2015 have been exacerbated by the continuing population decline at Bird Island (Forcada & Hoffman, 2014).The same pattern most likely occurs in nearby areas of South Georgia, where the species shares common foraging grounds and is affected by common climate-environment drivers (Meredith et al., 2008).This includes an unprecedented declining trend in first-year survival, which subsequently affects age at maturation, recruitment, mean breeding age and longevity.As a result, the population of Bird Island, which comprises approximately 5% of the global population, has declined by over 42% in less than one fur seal generation (9-11 years).This decline has been even more precipitous at the South Shetland Islands (Krause et al., 2022), South Georgia is in the Scotia Sea, which holds over half of the estimated global Antarctic krill biomass (Siegel & Watkins, 2016), is also the area with the highest regulated commercial krill extraction (Nicol & Foster, 2016), and where greenhouse warming is driving increases in ocean temperatures, changes in sea-ice seasonality and extent, increases in UV radiation and ocean acidification, affecting the marine biota (Constable et al., 2014).Here, the krill distribution is likely to be contracting Poleward, while increasing migration towards the seabed (Atkinson et al., 2019;McBride et al., 2021;Schmidt et al., 2011), both of which decrease its availability to millions of highly dependent predators such as Antarctic fur seals.
In the wider krill fishery management area encompassing South Georgia, CCAMLR Area 48.3, the commercial extraction of Antarctic krill peaked in the late 1980s, declined in the early 1990s (Figure 4a) with the collapse of the Soviet Union fishing fleet, and significantly increased since then (Nicol & Foster, 2016), although to much lower levels.Thus, the cumulative extraction between 1985 and 1991 was 1,108,473 tonnes, or between 6.4% and 7.1% of the available krill biomass estimates of 17,211,300 and 15,563,986 for that area for 2000 and 2019 (Hill et al., 2016;Krafft et al., 2021).Whereas subsequently, the extracted percentage has reduced to a mean of around 0.67%-0.74%per annum.This extraction is temporally uncorrelated with fur seal population growth, and is probably marginal compared to the fur seal peak demand on krill.This suggests that commercial krill extraction in the area encompassing a large part of the female Antarctic fur seal foraging summer and winter ranges (Arthur et al., 2017;Bamford et al., 2021;Boyd et al., 1998;Staniland & Robinson, 2008) could be sustainable.However, the effects of krill extraction on fur seal populations should be re-evaluated in the light of increasing environmental pressures on both the krill and the fur seals.Besides climate warming, these pressures include the effects of regional humpback whale recovery (Zerbini et al., 2019) and its direct competition for food with fur seals and other abundant krilldependent predators such as macaroni penguins.The present humpback whale annual krill consumption could be between 7% and 10% of the total krill stock (data from Baines et al., 2021), but data on the effects on fur seals and the wider ecosystem are largely unavailable.
We have shown that accurate and representative abundance estimates are the key to evaluating the recovery of depleted species, not only retrospectively, but also from the contemporary perspective of greenhouse warming, habitat degradation and shifting harvesting practices, which increasingly affect marine food webs (Daskalov et al., 2007;du Pontavice et al., 2020;Hocevar & Kuparinen, 2021).
We show that some of these effects have resulted in transient population states in Antarctic fur seals at South Georgia via changes in population structure and fertility, which carry over as negative population momentum (Forcada et al., 2008) and limit the ability of this population to recover from perturbations while facing new external pressures.As the fur seal recovery trajectory is expected to be fundamentally different from the long-term population at equilibrium (Hastings et al., 2018), previous projections of observed fur seal populations based on previous transient states, particularly without corroborating survey data, are very inaccurate, with an estimated upward bias of between 0.6 and 2.7 million seals (e.g.Bonner, 1985;Boyd, 2002;Foley & Lynch, 2019;SCAR, 2006).Therefore, these should not replace current abundance data in support of policy, and to achieve management objectives.

ACK N O WLE D G E M ENTS
We thank the British Antarctic Survey (BAS) field assistants and counts were obtained based on digital SLR camera images.Hereafter, an Antarctic fur seal breeding season is referred to by the second year of the split-year austral summer (e.g.2007-2008 is referred to as 2008).
First year female survival was lowest in 2006 and 2010 after strong climatic anomalies.For 2002-2015, the temporal mean was 0.441 [0.388, 0.493],

II
P1 I P2 I P3 I P4 I P5 I P6 I P7 I P8 I B I N P1 I P2 I P3 I P4 I P5 I P6 I P7 I P8 I P9 I P0 I P11 I T I N Annual medians and 95% credible intervals of vital and life history rates, and abundance at Special Study Beach (representative population) estimated from the integrated population models.(a) Apparent survival of experienced breeders (magenta) and skippers (cyan) and the immigration rate of breeders (black); (b) apparent survival of first year females, including a moving average filtering of order three (points, lines and polygon in magenta), and the fecundity rate (black); (c) seven-year rolling window estimates of longevity, the generation time and the observed mean age at breeding; and (d) abundance by demographic group at the representative female population.

From
2010 to 2015 the percentage of breeders declined significantly by 12.6% with respect to 2001-2009, and the proportions of skippers and pre-breeders increased.After 2015, breeders increased by 11.4% and 22.6% in relation to 2010-2015 and 2001-2009 respectively, while pre-breeders declined by 39.3% and 34.7% for the same periods, and the proportions of skippers remained stable throughout.The SSB female population declined from 2002 to 2022 (Figure 2d), with an estimated mean λ of 0.948 [0.938, 0.958].The fastest decline (λ = 0.927 [0.904, 0.953]) occurred from 2016 to 2022 (Figure 2d), when the survival of breeders and first year females decreased significantly (Figure 2b).This resulted in very low numbers of older juveniles (ages 1-7) in subsequent years.Extending the data series back to 1984, the long-term mean λ for mature females was 0.972 [0.966, 0.977].Most of the decline occurred during the last 13 years (2010-2022) of the 39year series (Figures 2d and 4b), during which the mean λ was 0.930 [0.917, 0.944], corresponding to an annual decline of 7%.Previously, the time series was characterised by two phases of near stability, from 1984 to 1990 (mean λ = 0.983 [0.957, 1.012]), and from 1991 to 2009 (mean λ = 0.997 [0.991, 1.004]), with an abrupt decline in 1991.Declines in ln(λ) were best explained by significant increases in SAM, summer SST and mean krill size (Figure4c-e).The SAM has been mostly in its positive phase since the early 1990s, with a steady decline in variation while increasing its mean; and summer SST (December-February) increased by 0.88°C over 40 years (Figure4a).The time-frequency map of coherence (correlation) between these indices showed a time-varying relationship between coupled fluctuations in climate and the physical environment, which propagated to mean krill size and ln(λ) (Figure 4f-h), especially during the decade preceding the abrupt 2009 decline.Under similar climatetemperature coupling, with current values of mean summer SST(3.34°C[3.03, 3.65]  since 2010) and highly positive SAM values, the mean ln(λ) is expected to be mostly negative and highly variable (Figure4d,e).
[0.902, 0.997].Additionally, the correlation between the Bird Island counts and total SSB productivity was 0.875 [0.444, 0.977].The correlation for paired predicted SSB counts and counts at Maiviken (−54.2435°,−36.5033°), another long-term monitoring site situated TA B L E 1 Female population structure at the Special Study Beach (SSB) of Bird Island, South Georgia.In each season, the total number of breeders is the sum of new breeders born at SSB (Recruits), returning breeders (Breeders) and immigrants from other areas.Skippers or non-breeders are mature breeders which defer breeding in each season.In parentheses are 95% credible intervals.approximately 105 km to the east, on the South Georgia mainland (data from Hollyman, Trathan, & Collins, 2021), at peak pupping over 11 years was 0.899 [95% CI: 0.650, 0.974].Thus, there is a high level of synchrony in the counts and productivity across South Georgia, supporting the applicability of SSB correction factors throughout the entire range.

3. 5 |
Mature female abundance and pup production at Bird Island Bird Island estimates of the absolute abundance of mature females and pup production by year are shown in Figure 4b (also Appendix S4; Table S1).Abundance fluctuated between 1989 and 2009, F I G U R E 3 Seasonal daily counts of mature females, newborn pups and territorial males at the representative population (Special Study Beach) for the years of aerial surveys of South Georgia.Top and bottom panels show observed counts in red/orange, and the best fit of a Poisson generalised additive model (blue lines) and an overdispersed Poisson generalised additive mixed model (green lines and grey polygon).Central panels show daily newborn pups (blue) with a fitted logistic model (black) of the cumulative daily number.Vertical solid red lines delimit the survey windows and dashed red lines the estimated peak pupping day.a) Increases in the Southern Annular Mode (SAM) index, with the green line showing the coefficient of variation of a 6year sliding window; South Georgia sea surface temperature (SST; dashed lines for baseline means); mean krill size in the seal's diet; and commercial krill extraction, with the collapse of the USSR fishery marked by the vertical dashed black line.Lines in the krill data indicate significant trends.(b) Median abundance mature females and breeding females (i.e.pup production) at Bird Island (BI) in northwest South Georgia, and at its representative study population (Special Study Beach).The grey polygon is a moving average filter of order three; and the cyan line is a segmented regression model with significant (solid lines) and non-significant trends (dashed lines).Vertical bars are 95% credible intervals.(c-e) The relationships between fur seal population growth rate (ln(λ)) and krill size, SST and SAM; grey polygons show estimates from fitted generalised additive models.(f-h) Time-varying relationships between SAM, SST, krill size and ln(λ) obtained with a wavelet transform multivariate analysis; contours show levels of synchrony (ρ ϵ [0, 1]) in cycles for combinations of variables.Lighter shades are cones of influence containing less reliable values.

F
I G U R E 5 Distribution and range expansion of the female Antarctic fur seal breeding population of South Georgia based on available comprehensive surveys in 1975-1976 [(a) Payne (1977)], 1990-1991 [(b) Boyd (1993)] and 2007-2009 [(c) this analysis].In (a), the location arrow points to Bird Island, identified in the early 1930s as the probable epicentre of the recovery of the species after sealing.Values for panels (a) (*) and (b) (**) are from the literature cited.
and has also been observed to a lesser extent in Bouvet Island(Hofmeyr, 2016), although present data for Bouvet are lacking.Nonetheless, this increases conservation concerns for the Southern Ocean's South Western Atlantic region, which includes three of the four main island groups that have distinct population genetic structures(Humble et al., 2018;Paijmans et al., 2020), and where the species is globally most abundant, but declining.
scientists on Bird Island and the South Georgia mainland for collecting tissue samples and demographic, biometric and count data over four decades.We thank research assistants Emma Foster, Alice Ramsey, Rebecca Hewitt, Eilidh Siegal and other personnel involved in the counting of aerial survey images.Andreas Cziferszky developed specific image counting software.The crew of Royal Navy HMS Endurance (A171) and pilots of the ice-modified Lynx helicopters were instrumental in completing the surveys of South Georgia.Anthony R. Martin contributed to initial planning and development of the surveys.This work contributes to the Ecosystems project of the BAS, Natural Environment Research Council, UKRI, and is part of the Polar Science for Planet Earth programme.The genetic work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the framework of a Sonderforschungsbereich (NC3 project numbers 316099922 and 396774617-TRR 212) and the priority programme 'Antarctic