) Apparent breeding success drives long-term population dynamics of a migratory swan

and found a significant association between juvenile survival both with the water level in lakes during autumn migration, which affects food accessibility for the swans, and with summer temperatures. Such associations are important for understanding the dynamics of species with fluctuating population sizes, and thus for informing management and conservation decisions.


Introduction
In our rapidly changing world, many species need to be able to adapt promptly to environmental change if they are to continue to be viable (Williams et al. 2008, Chen et al. 2011, Hoffmann and Sgró 2011).This is particularly challenging for migratory species, because of the potential mismatch between the (rates of ) change in conditions encountered by individuals at sites used along their migration routes (Visser et al. 2009, Bairlein 2016).The ability of a species to adapt is ultimately reflected in its vital rates (i.e.survival and reproductive success), which together determine a population's increase or decline.Understanding how species respond and how population trends differ under varying environmental conditions therefore is of great interest, from both a scientific and a conservation perspective (Oppel et al. 2014).
Population counts often form the basis for conservation and management strategies in an attempt to conserve or control, either positively or negatively, the abundance of a species.However, for many species, their conservation status remains uncertain, despite large-scale and labour-intensive monitoring programmes (Tempel et al. 2014).This may reflect the elusiveness of the species, or its slow life-history traits which would require long-term monitoring to determine the reasons for population change (Cardillo et al. 2005), making it difficult to assess its status effectively.Fully understanding the dynamics of these populations can be challenging (Robinson et al. 2014), although analysing the underlying vital rates for such species can provide valuable insights (Lebreton and Clobert 1991, Besbeas et al. 2002, Katzner et al. 2006).Determining the main demographic processes steering population dynamics is also valuable for helping to identify the location, time of year (e.g.breeding or non-breeding season) and spatial scale at which conservation strategies need to be implemented, especially for migratory species (Schaub et al. 2012).
Integrated population models (IPMs) were developed initially to exploit the fact that census data contain information on the survival and reproduction of individuals and vice versa (Besbeas et al. 2002, Schaub andAbadi 2011).Rather than analysing census data and demographic data separately (such as survival analyses based on capture-mark-recapture data), integrating these measures permits information flow between demographic processes, accounting for imprecision of data and correlation between estimators (Besbeas et al. 2002, Schaub andAbadi 2011).The method improves the accuracy of parameter estimation and allows estimation of parameters for which no explicit data are available (Tavecchia et al. 2009, Schaub and Abadi 2011, Altwegg et al. 2014).IPMs in the Bayesian framework are based on a state-space model (Brooks et al. 2004): time-series models in which the latent state process (i.e. the true population size) is linked to population counts via an observation process, thus accounting for observation error (Brooks et al. 2004).With these advances in the modelling process, IPMs are increasingly used in conservation science (Zipkin and Saunders 2018).
In long-lived species, adult survival is often the main contributor to population growth (Saether and Bakke 2000) and this relative contribution, in comparison with other vital rates, increases with generation time (i.e. the average difference in age between parents and their newly-hatched offspring; Koons et al. 2014).At the same time, there is often little variation found in the particular traits that make the largest contribution to the rate of population growth (Gaillard et al. 1998, Pfister 1998, Saether and Bakke 2000, Koons et al. 2014, but see Morris andDoak 2004, Jongejans et al. 2010).Greater longevity leads to greater complexity in the age structure of the population, resulting in population growth rates for species with slow life histories (e.g.long life span and generation time, and high adult survival rates) being more sensitive to changes in their stable age distribution than shorter-lived species (Koons et al. 2014).
The Bewick's swan Cygnus columbianus bewickii is a large, long-lived migratory bird species that breeds on tundra habitat across the Russian Arctic (Rees 2006).The northwestern European population, which winters in northwest Europe, has shown a variable trend in numbers since the mid-20th century.Between 1970 and 1995 numbers increased rapidly from ca 15 000 to about 30 000 birds (Beekman et al. 2019).Similarly increasing trends observed for other European swan and goose populations have been attributed to the birds' increased use of agricultural land improving their winter survival and breeding success, combined with a reduction in mortality resulting from hunting bans (Van Eerden et al. 2005, Fox et al. 2010, Rees and Beekman 2010).In contrast to many of these populations, Bewick's swan numbers decreased again, to 21 500 birds in 2005 and 18 000 in 2010 (Rees andBeekman 2010, Beekman et al. 2019).The most recent international count indicated an increase in population size with around 20 500 birds in 2015, but results of the January 2020 census are required to evaluate the population trajectory (Beekman et al. 2019).The population is thought to be sensitive to climate and habitat changes, and to illegal or accidental shooting, though the causes of these population trends and the mechanisms underlying the population dynamics remain unclear (Newth et al. 2011, Nagy et al. 2012).
Given its small population size and recent declining population trend, the northwestern European population of Bewick's swans is now considered endangered and is included in the European Red List for birds of conservation concern (Birdlife International 2015, Rees et al. 2019).The species is also listed as a target within the Natura 2000 framework (Annex I Birds Directive 2009/147/EC of the European Parliament).In response to the decline, experts on the species from across the flyway met to discuss and develop an international single species action plan (ISSAP) for the population; the resulting ISSAP was adopted by the AEWA in 2012, with the overall aim of halting the ongoing decline in the short-term and maintaining the population minimally at its 2000 level (23 000 birds) in the long-term (Nagy et al. 2012).Several hypotheses for the decline are listed in the ISSAP, and the plan also calls for continued international collaborative monitoring and research into demographic parameters and population change, in order to understand and thus address the causes of the decline.
Previous efforts to study the population dynamics of the northwestern European population of Bewick's swans found no trends in breeding success (neither in the percentage of juveniles in the population nor in brood size) over the period 1964-2015, although fewer young were observed in winters following a colder summer or a higher abundance of Arctic fox Vulpes lagopus in the breeding area (Wood et al. 2016).Similarly, no consistent temporal trend was found in survival over the years 1970-2015 (Wood et al. 2018).Survival rates over the study period varied between decades, with mean survival rates being highest in the 1980s and lowest in the 2010s (Wood et al. 2018) in accordance with the population trends described by the count data.However, neither study provided a comprehensive explanation for the population trend.
In this study, we therefore developed an integrated population model (Abadi et al. 2010, Kéry andSchaub 2011) to study the Bewick's swans' population dynamics and underlying vital rates.By combining three different long-term datasets we aimed to obtain more precise estimates of the swans' survival and breeding success during periods of positive and negative population growth.In addition, we performed a transient life table response experiment to assess which changes in age structure and vital rates contributed most to variation in population growth.Finally, we also explored the environmental variables that might explain variation in the Bewick's swans' vital rates, to determine where conservation action might be required.

Data collection
Three different types of data sources were used in this analysis: census data, capture-mark-resighting (CMR) data and juvenile proportions.All data were recorded in the Bewick's swans' wintering range (northwest Europe), as described below.

Census data
Internationally coordinated Bewick's swan counts have been carried out since the early 1970s, initially through the International Waterbird Censuses coordinated by Wetlands International in 1973International in , 1976International in , 1979, then as special surveys focussing on Bewick's swans in 1984 and 1987 (Dirksen and Beekman 1991, Rees 2006, Nagy et al. 2012).From 1990 onwards, the Wetlands International/IUCN SSC Swan Specialist Group has organised an International Swan Census every five years in mid-January (Beekman et al. 2019, Laubek et al. 2019).To estimate vital rates annually in the IPM, numbers in the population between census years were imputed, as IPMs without imputed population counts did not converge nor produced reliable year-specific estimates.Imputed values were drawn from a normal distribution with mean and standard deviation based on counts in the nearest two census years to mimic natural fluctuations in population size (Zhao et al. 2019).

Capture-mark-resighting data
Data from two different ringing schemes were used in this study.In the first, undertaken by the Wildfowl and Wetlands Trust (WWT) in the United Kingdom since the mid-1960s, Bewick's swans are ringed with a plastic leg ring, each engraved with a unique alphanumeric code (Ogilvie 1972).The second scheme was initiated by author TH (in 1989) and continued by the Netherlands Institute of Ecology (NIOO-KNAW) (since 2005).In this scheme, Bewick's swans are ringed with a neck band with alphanumeric codes (and especially since 2005 also with leg rings with identical codes).The WWT dataset used in this study  consists of 830 individuals marked as juveniles (i.e.first winter birds, determined by their grey plumage), 529 marked as yearlings (birds in their second winter, identified by traces of grey juvenile feathering in otherwise white adult plumage) and 1833 marked as adults (in their third winter or older).The NIOO dataset (1989-2015) contained 79 individuals marked as juveniles, 33 as yearlings and 432 as adults.Ringed birds have been reported by volunteer bird watchers from across northwest Europe.For all years following capture and marking, we determined for each individual whether it had been resighted at least once per winter, between the months of December and February inclusive.This resulted in individual leg-ringed birds being recorded for a total of 13 358 winters and neckbanded birds recorded for 2201 winters over the study period.From these encounter histories, resighting probabilities were modelled for leg rings (p l ) and neck bands (p n ) separately, because the former were known to have lower resightings rates (Wood et al. 2018;Supplementary information).

Winter reproduction data
Long-term reproduction data from the breeding grounds is not available for Bewick's swans.Therefore, we used winter brood size (i.e. the average number of juveniles per successful breeding pair) and proportion of juveniles in the population, recorded in winter, as a proxy for productivity.Both variables were monitored annually in the UK (Wood et al. 2016) and in the Netherlands (Poorter 1991, Hornman et al. 2020) in mid-January and late November/early December, respectively.These two countries receive 60-90% of the total population in mid-winter, and no consistent geographical trend in juvenile proportions has been detected (Beekman et al. 2019).We therefore took the mean of their national productivity estimates, for each of the two variables, as representative of annual breeding success for the population as a whole.The annual proportion of juveniles was used to create age-structured counts.Winter brood size was used to create an informative prior for the estimate of apparent breeding success κ (below).

Integrated population modelling
We developed a female-based, age-structured population model with three age classes (i.e.juvenile, yearling and adult).We modelled (female) apparent breeding success κ as a latent parameter.Because counts contain information about survival and reproduction, integrating counts, annual proportion of juveniles and age class-specific survival in a single model allowed us to estimate apparent breeding success, despite the lack of reproduction data from the breeding grounds.As such, κ reflects breeding propensity, breeding success and post-fledging survival until winter, all in one parameter.
The CMR data were analysed as individual encounter histories and modelled with a Cormack-Jolly-Seber (CJS) model (Lebreton et al. 1992).Apparent survival from year t to t + 1 was modelled on the logit scale as a function of age class and time.Resighting probability in year t was modelled on the logit scale as a function of marker type and time.In both rates (survival and resighting probability), environmental stochasticity was included by considering time as a random year effect.
The census data were analysed using a state-space model (Brooks et al. 2004).The state-space model consists of a state process and an observation process.The state process describes the true but unknown population trajectory and accounts for demographic stochasticity by describing the number of juveniles via a Poisson distribution and the number of yearlings and adults via a binomial distribution.The observation process describes the link between the true and the observed population size (De Valpine and Hastings 2002).Using the data on the proportion of juveniles, we divided the counts into juveniles and older individuals.Assuming an even sex ratio, we halved counts before implementation in the femalebased model.
The population transition matrix of the female-based, agestructured population model is the following: where W j , W y and W a are the number of juvenile, yearling and adult females in winter, respectively, φ x is the survival rate between year t and t + 1 for females in age class x, and κ is the apparent breeding success in year t: the average number of female fledglings per successful breeding female, that survived until winter.
To estimate age class-specific numbers and demographic parameters, the joint likelihood is composed of the likelihood of the state-space model for the population census data and the CJS model for the individual encounter histories.The IPM is analysed in the Bayesian framework, which combines the joint likelihood with prior probability distributions to obtain posterior distributions of the parameters.We specified non-informative priors for all parameters but apparent breeding success.The prior for apparent breeding success was a uniform distribution between 0 and 2 (i.e. the average winter brood size (Wood et al. 2016) that is set to be the upper limit of κ).Markov chain Monte Carlo (MCMC) sampling was used to simulate observations from the posterior distributions with JAGS ver.4.3.0(Plummer 2003) run from R with jagsUI ver.1.5.1 (Kellner 2019).We ran three chains of 200 000 iterations with a burn-in of 100 000 and a thinning rate of 50, resulting in a total of 6000 posterior samples.Convergence of the MCMC chains was evaluated by ensuring that the Brooks-Rubin-Gelman diagnostic R for each parameter was below 1.1 (Brooks and Gelman 1998).
A more detailed description of the IPM, including likelihoods and priors, can be found in Supplementary information.R and JAGS code of the IPM can be found in Supplementary information.

Retrospective perturbation analysis
Retrospective perturbation analyses aim to decompose the effects of observed variation in demographic parameters on past population dynamics (Caswell 2000).Traditionally, the focus has been placed on identifying the demographic rate with the strongest influence on the asymptotic growth rate λ (Koons et al. 2016).These analyses, however, ignore the feedback between population structure and demographic rates (Tuljapurkar 1990) and short-term transient dynamics (Stott et al. 2011) that are common for many contemporary populations (Koons et al. 2016).Here, we performed a transient life table response experiment (LTRE) (Koons et al. 2016(Koons et al. , 2017) ) to determine the contributions of changing vital rates and population structure to the variation in the realized population growth rate, λ t = W t+1 /W t .A more detailed description of the transient LTRE can be found in Supplementary information.

Environmental variables
We collated data on several potential explanatory variables for the variation in vital rates based on previous research on Arctic-nesting waterfowl (Skinner et al. 1998, Trinder et al. 2009, Morrissette et al. 2010).In order to prevent collinearity, we condensed variables wherever possible (Morrissette et al. 2010).Weather-related variables included were the cumulative negative degree days (ndd, °C) and cumulative degree days (cdd, °C) in Naryan-Mar (Russia, i.e. the breeding range) and average tailwind in the Baltic region (58°N, 23°E).Ndd was comprised of the summed mean daily negative air temperatures in late spring (1 May-30 June), corresponding to the nesting and incubation period (Rees 2006).Cdd was comprised of the summed mean daily positive air temperatures during summer and autumn (1 July-31 October), corresponding to the period between hatching and fledging and set to include the frost days in autumn (Rees 2006).Average tailwind was collated for autumn (18 September-31 October), corresponding to the autumn migration period (Rees 2006), and calculated from average north-south and east-west windspeeds, multiplied with the cosines of the heading in degrees (225°).Temperature data were acquired from the NOAA Global Summary of the Day dataset (NOAA National Climatic Data Centre 2019) and wind data from the NOAA NCEP/NCAR reanalysis project (Kalnay et al. 1996).
To account for a potential effect of density dependence in the population, we included winter totals as estimated by the IPM as a predictor.In addition, we included whooper swan population estimates for the northwest mainland European population (Netherlands, Germany and Denmark combined), derived from the International Waterbird Census over the study period (Wetlands International 2019), as it was hypothesised that competition with this species might be an influence on the declining population trend recorded for Bewick's swans since 1995 (Nagy et al. 2012).The Lake Peipsi (Estonian-Russian border, 58.5°N, 27.5°E) water levels for each year of the study (average September-October), measured by the Estonian Meteorological and Hydrological Institute, were also included, as water levels determine the availability of macrophytes in this key stop-over site and other lakes in the region for Bewick's swans arriving from the Russian Arctic during autumn migration (Luigujõe et al. 1996).
We ran a series of generalized linear models (GLMs) with the posterior means of the IPM output parameters (adult, yearling and juvenile survival, apparent breeding success) as response variables, and the environmental variables as predictors.We included year as a predictor variable in the analysis to cover all environmental variables that changed monotonically over time but that we were unable to capture in data.Most predictor variables were not correlated (Pearson correlation coefficient r < 0.3), apart from year, Bewick's swan estimates and the whooper swan estimates, which were highly correlated (year and whooper swan r = 0.96; year and Bewick's swan r = 0.60; Bewick's swan and whooper swan r = 0.60).We therefore ran all models three times: once with year, once with the whooper swan estimates, and once with Bewick's swan estimates included.All variables were centered and scaled to a mean of 0 and sd of 1 before running the model.We compared the AIC of the models to choose the most parsimonious (i.e. the one with the lowest AIC).All analyses were run in R ver.3.6.2(<www.r-project.org>).

Contributions to variation in population growth rate
The mean realized population growth rate was 1.035 (±0.052SD) in period 1 and 0.983 (±0.041SD) in period 2. The two demographic rates that contributed most to between-year variation in the population growth rate λ t (var(λ t ): 0.003), were apparent breeding success κ and adult survival φ a (Fig. 2).The contribution of juvenile survival φ j and yearling survival φ y to variation in λ t was positive but small, whilst the effect of the population structure components were negative and also small (Fig. 2).Based on non-overlap of the posterior means with zero, all variables contributed significantly to the variation in realized population growth rate, except for φ j , and φ y .Vital rates contributed 94.6% of the total contribution.Supplementary information for numerical estimates.

Environmental variables
Comparing the model outputs of the three different sets of GLMs (Supplementary information) based on their AIC values (Burnham and Anderson 2002), showed that, although differences were minor and often within ΔAIC of two, the model including whooper swan estimate was the most parsimonious for adult survival and apparent breeding success.The model with Bewick's swan estimate and year were the most parsimonious for the GLMs with yearling survival and juvenile survival as the response variable, respectively (Supplementary information).For adult and yearling survival, we found no significant effects of any of the predictor variables.Year, the water level in Lake Peipsi and ndd showed a significant positive association with juvenile survival (higher juvenile survival in more recent years, with higher water levels in autumn and with colder summers).Apparent breeding success was negatively associated with whooper swan estimate (Table 1).

Discussion
By integrating three long-term datasets available for the northwestern European Bewick's swan population (census data, CMR data and the proportion of juveniles present in winter) we were able to obtain precise estimates for the survival rates for three different age classes (Abadi et al. 2010).In addition, we used the possibility of the IPM framework to synthesize parameters for which no specific data were available (Abadi et al. 2010, Altwegg et al. 2014), namely breeding propensity, hatching success and post-fledging survival (until winter), to generate an estimate of annual apparent breeding success κ.We show that adult survival and κ were the vital rates that contributed most to variation in the population trends (growth and subsequent decline) observed over the last five decades (Fig. 2).
Although population growth rates of long-lived species are in theory particularly sensitive to changes in adult survival, this vital rate is often stable over time.Other demographic rates are more sensitive to environmental fluctuations and thus, in effect, contribute more to the population growth rate (Koons et al. 2014).In our species we indeed find that adult survival made a high contribution to variation in the rates of population growth, but varied little over time.This is in line with the Demographic Buffering Hypothesis (Boyce et al. 2006, Hilde et al. 2020), but note that no proof for this hypothesis was found in a comparative study among waterfowl species (Koons et al. 2014).More between-year variation was visible in the juvenile survival and κ estimates, for the latter especially between the two periods with contrasting growth rates, i.e. when the population was increasing during the 1970s to mid-1990s and the decrease in numbers thereafter (Fig. 1a,  d).The same pattern (i.e.high contribution of adult survival to the population growth rate but low variation in this parameter, and more variation in the fecundity-related parameters such as nest success and juvenile survival) has also been found in other studies of waterfowl (Cooch et al. Hoekman et al. 2002, Bentzen and Powell 2012, Wilson et al. 2012).On considering potential explanatory variables for the variation in demographic rates, we found that juvenile survival increased significantly over time (year).The year variable however can also be seen as a proxy for all environmental variables that have changed over the years without being captured in our set of variables.This includes for example the migration distance of the northwestern European Bewick's swan population, which has decreased gradually since the early 1970s (Nuijten et al. 2020b).Other variables that might have changed over the same period, such as the abundance of predators on the Arctic breeding and moulting grounds (Bêty et al. 2002, Smith et al. 2010, Nolet et al. 2013), could also play a role.The increase in juvenile survival coincides with the swans switching from feeding on aquatic resources to agricultural resources (Merne 1972, Mullié andPoorter 1977).Juveniles have lower intake rates than adults when feeding on below-ground part of aquatic macrophytes, and hence are thought to benefit from this change in diet, more so than adult birds (Nolet et al. 2014).The increase in juvenile survival occurred in the period of population increase (< 1995).One factor known to have changed over time is whooper swan abundance (Supplementary information).The northwestern mainland European population of whooper swans has been increasing steadily in recent decades (Laubek et al. 2019).Being larger, the whooper swan is likely to win competitive encounters with Bewick's swans over food (Black and Rees 1984) or nesting places.Competition for food outside the breeding season may lead to a lower breeding success, by influencing body condition prior to the breeding season (Hoye et al. 2012, Koyama et al. 2013).It has been found that food abundance in the wintering area is not limiting, even when competition is present (Wood et al. 2019), but food limitation is more likely at migratory sites where swans concentrate and food depletion has been demonstrated (Nolet and Drent 1998).Previous research has shown that the two species often overlap on their spring staging areas in Estonia (Luigujõe et al. 2002), and potentially also on the Gulf of Finland near St Petersburg (Zaynagutdinova et al. 2019).Aerial surveys undertaken between 1980 and 2000 of the Pechora Delta and the adjacent Russkiy Zavorot Peninsula, considered to be one of the main breeding areas for Bewick's swans in European Russia (Rees 2006), found an increase in the number of territorial whooper swan pairs breeding in the region, whereas the number of Bewick's swan pairs was stable (Shchadilov 2002).Further studies are needed to investigate in greater detail the potential effects on Bewick's swans of interspecific competition for resources in different parts of the northwestern European flyway.
In addition, both the water level in Lake Peipsi in autumn and the cumulative negative degree days (ndd) in summer (both continuous variables) were significantly and positively associated with annual estimates of juvenile survival.Lake Peipsi comprises an important stop-over site for Bewick's swans in autumn (Luigujõe et al. 1996).To be able to reach the belowground parts of macrophytes in the substrate of the lake, the water level must not be too high (Luigujõe and Kuresoo 2007).Observations in 2005 and 2006, years with a very high and very low water table respectively, convincingly showed the strong relationship of Bewick's swan numbers with water tables: in 2005 smaller groups visited far fewer sites in autumn, while in 2006 the total autumn count reached up to 3000 on both sides of the lake (Luigujõe and Kuresoo 2007).Autumn water levels are heavily dependent on rainfall, and high water levels in Lake Peipsi are likely to be representative of high water levels in surrounding lakes; hence, aquatic foraging conditions would be poor in the whole region when water levels are high, and it was suggested that this would result in slower autumn migration by Bewick's swans through the region (Beekman et al. 2002).We therefore expected that there would be a negative relationship between the water levels and (juvenile) survival, but this hypothesis was not upheld.Our satellite tracking in an autumn with extremely high water levels in Peipsi (autumn 1998) revealed that the Bewick's swans visited many alternative sites, while slowly continuing their migration (Beekman et al. 2002).This might be beneficial for juveniles, perhaps through less competition at specific sites.A similar reasoning could explain the effect of ndd: in years with low summer temperatures, the reduction in apparent breeding success might lead to lower competition between the hatchlings, and to higher juvenile survival rates the following year.Another plausible explanation for the positive association of juvenile survival with ndd, is that in harsh autumns many fledglings die (Wood et al. 2016), leaving only the strongest individuals in the juvenile phase.However, we did not find a negative association between apparent breeding success and ndd to substantiate these possible explanations (Table 1).
Previous studies have found that Bewick's swans seem to profit from favourable wind conditions during spring migration (Klaassen et al. 2004).In line with this, we expected that headwinds encountered during autumn migration would have negative effects on (juvenile) survival.Our analyses however found no correlation between headwinds and time, and hence no evidence of any long-term trend in headwinds (Supplementary information).It remains to be seen, therefore, whether this factor affects the longer-term population dynamics of Bewick's swans in northwest Europe.
Negative density dependence is regarded as an important driver of population growth rates (Koons et al. 2014).If influential, it can have both negative and positive effects on the demographic rates when population numbers are high or low, respectively (Turchin 1995) and effects can vary for different age classes (Layton-Matthews et al. 2019).Results from the LTRE however, show only a small negative effect of components of the population size on the variation in population growth rate over the study period, of which only the effect of adults W a was significant (Fig. 2).In addition, Bewick's swan population estimates did not proof to be an important predictor of the vital rates in the GLM analysis (Table 1).These results are in line with previous results where it was found that densitydependence performed poorly as a predictor of Bewick's swan survival estimates (Wood et al. 2018).
An assumption of the IPM framework is the independence of the underlying datasets (Kéry and Schaub 2011).Although it is possible that this study violates this assumption, because every dataset originates from the same population in the wintering grounds and individuals can appear in multiple datasets, we argue that this assumption is valid for our input data.The international censuses are conducted on a five-year basis, in one weekend in January, the juvenile proportions are determined separately and on an annual basis at a fixed time in winter, and the resightings of marked individuals are collected by individual observers throughout the winter season on any location.Besides, it has been shown that violation of this assumption has little effect on estimation of parameters (Abadi et al. 2010, Schaub andFletcher 2015), in particular when the number of individually marked birds is a small proportion of the studied population (Weegman et al. 2016).
As part of the IPM, we ran a CJS model to obtain survival estimates from the encounter histories.In traditional survival analysis, usually a goodness-of-fit test is conducted, but we are currently unable to find a means of implementing this within the IPM framework (Kéry and Schaub 2011; but see (Besbeas and Morgan 2014).In a previous, traditional survival analysis based on the same data, a goodness-of-fit test revealed evidence for both trap-dependence (trap-happiness) and transience (Wood et al. 2018).However, the ĉ value (2.5; calculated in U-Care, Choquet et al. 2009) was below the maximum ĉ value of 3.0 recommended for reliable survival modelling (Lebreton et al. 1992).In addition, Abadi et al. (2013) showed that capture heterogeneity did not have a large effect on the resulting estimates.
In the model, all individuals in the adult age class considered to be capable of breeding, with a certain breeding propensity (included in κ).We know from observations (individuals that returned to the wintering grounds with juveniles) that, although birds of three years old do reproduce, most individuals start reproduction later in life (4-5 years of age ;Rees 2006).This simplification in the model means that κ underestimates the apparent breeding success of potential breeders (i.e.≥ 4-year-olds), while still providing a useful estimate for the whole population of ≥ 3-year-olds.We do not know whether the distribution of age of first breeding has changed over time, which might occur if an increase in population size and breeding densities is resulting in density-dependent regulation of the population in the breeding range (Koons et al. 2014).
In our study, apparent breeding success κ was modelled as a latent state variable (i.e. it was inferred in our IPM rather than informed by data).This was necessary as there are currently no long-term data on the productivity of the northwest European Bewick's swan population that are measured on the Arctic breeding grounds.Given our finding that apparent breeding success was critical to explaining Bewick's swan population dynamics, our work reinforces the need for improved long-term surveillance of measures of swan productivity in the breeding range, as suggested in the Bewick's swan AEWA action plan (Nagy et al. 2012).Such measures would allow future IPMs to be parameterized using counts, CMR and productivity data sets together, which would improve accuracy of the model estimates.The availability of these measures would also enable researchers to test for temporal trends directly, and would complement the age ratio data that is already collected in the wintering areas (Wood et al. 2016).We think that our measure of apparent breeding success, κ, was more representative of actual breeding success than the productivity parameters recorded in winter because it was based on the summer population of potentially breeding females (Supplementary information).This may explain why we found a significant decline in reproductive success over the years, when previous analyses using the winter measurements have not done so (Wood et al. 2016).The importance of apparent breeding success as indicated by the LTRE (Fig. 2) must be treated with some care, as it is known that latent variables often seem to explain large proportions of the growth rates (Ahrestani et al. 2017, Saracco andRubenstein 2020), potentially because they may absorb systematic bias associated with the different datasets (Tavecchia et al. 2009, Taylor et al. 2018).
In conclusion, we found that declining reproductive success in this long-lived migratory species was the major driver of the observed population decline in north-western Europe.As such, our study advances one of the key aims of the AEWA action plan for Bewick's swans, namely to improve our understanding of the demographic causes of the documented decline in population size (Nagy et al. 2012).Future studies will be required to assess why apparent breeding success may have declined over time, and to identify important environmental drivers.As per the recommendations of the AEWA action plan, such studies ought to address potential habitat changes, predation and interspecific competition with other swan species within the Bewick's swan breeding and staging areas (Nagy et al. 2012).In addition, we found an increase in juvenile survival over time, associated with the water levels at important foraging sites in autumn.Our study highlights the importance of considering all information available on the population level to estimate precise demographic rates which can in turn be used to accurately assess the influence of potential explanatory variables on population growth rates.Such analyses are important in understanding the population dynamics of species with fluctuating population sizes, and to inform management and conservation decisions.

Figure 1 .
Figure 1.Estimates from the integrated population model of northwestern European population of Bewick's swans.(a-c) Annual apparent survival for (a) juveniles (φ j,t ), (b) yearlings (φ y,t ) and (c) adults (φ a,t ).(d) Annual apparent breeding success (κ t ; black points) and the proportion of juveniles in the population during winter (j, grey points).(e) Annual winter proportion of juveniles as estimated by the model (black points) and as observed in winter (grey points).(f ) Annual winter abundance as estimated by the model (black points), with actual (solid grey points) and imputed (open grey points) population counts.Estimates are represented by posterior means and 95% credible intervals.

Figure 2 .
Figure2.Contributions of the vital rates (i.e.apparent breeding success κ, apparent survival for juveniles φ j , yearlings φ y and adults φ a ) and normalised components of the population structure (i.e.winter abundance for juveniles w j , yearlings w y and adults w a ) to variation in realized population growth rate var(λ t ), determined through a transient life table response experiment (LTRE).

Table 1 .
Model output of GLM with demographic parameters estimated by the IPM as response variables (adult, yearling and juvenile survival, and apparent breeding success) and the explanatory variables as predictors (year, whooper swan estimate, Bewick's swan estimate, tailwind in the Baltic sea in autumn, water level in Lake Peipsi in autumn, negative degree days in Naryan-Mar in spring, cumulative degree days in Naryan-Mar in summer and autumn combined).The model including whooper swan estimate was the most parsimonious for adult survival and apparent breeding success.The model with Bewick's swan estimate was the most parsimonious for yearling survival and the model with year was the most parsimonious for juvenile survival.Differences in AIC were however minor (Supplementary information).The result of the other GLMs are presented in Supplementary information.Bold text indicates effects that are significantly different from zero.