The demographic causes of population change vary across four decades in a long‐lived shorebird

Abstract Understanding which factors cause populations to decline begins with identifying which parts of the life cycle, and which vital rates, have changed over time. However, in a world where humans are altering the environment both rapidly and in different ways, the demographic causes of decline likely vary over time. Identifying temporal variation in demographic causes of decline is crucial to assure that conservation actions target current and not past threats. However, this has rarely been studied as it requires long time series. Here we investigate how the demography of a long‐lived shorebird (the Eurasian Oystercatcher Haematopus ostralegus) has changed in the past four decades, resulting in a shift from stable dynamics to strong declines (−9% per year), and recently back to a modest decline. Since individuals of this species are likely to respond differently to environmental change, we captured individual heterogeneity through three state variables: age, breeding status, and lay date (using integral projection models). Timing of egg‐laying explained significant levels of variation in reproduction, with a parabolic relationship of maximal productivity near the average lay date. Reproduction explained most variation in population growth rates, largely due to poor nest success and hatchling survival. However, the demographic causes of decline have also been in flux over the last three decades: hatchling survival was low in the 2000s but improved in the 2010s, while adult survival declined in the 2000s and remains low today. Overall, the joint action of several key demographic variables explain the decline of the oystercatcher, and improvements in a single vital rate cannot halt the decline. Conservations actions will thus need to address threats occurring at different stages of the oystercatcher's life cycle. The dynamic nature of the threat landscape is further supported by the finding that the average individual no longer has the highest performance in the population, and emphasizes how individual heterogeneity in vital rates can play an important role in modulating population growth rates. Our results indicate that understanding population decline in the current era requires disentangling demographic mechanisms, individual variability, and their changes over time.

SECTION S1 -SURVIVAL S(Z) We performed a mark-recapture analysis to investigate survival and whether survival was related to the individual co-variates of age, breeding status and lay. Survival was estimated using a Multi-state model that included both live and dead recoveries (Appendix S1). Our analysis focused on individuals that were ringed on Schiermonnikoog and did not include individuals that had been ringed elsewhere. We only considered age classes and hence whether survival varied among fledglings (first year survival), sub-adult (second and third survival), pre-breeding ages (see Appendix S1) and adults (four year or above; Figure 1). Survival of adults was modelled against the lay date of their breeding attempts, which was estimated as the average lay date during the period of interest. In contrast, sub-adult survival was modelled against the lay date of the nest from which the individual was recruited. In order to know the lay date of the nest, sub-adults had to first survive from fledgling to age one, hence the same individuals were used to estimate fledgling survival (see Fecundity) and sub-adult survival (Appendix S1). As described in Appendix S1, the parameters for survival were obtained by model-averaging all considered models. Hence all age classes included a quadratic relationship with lay date and the logistic equation used in the IPM for survival took the form ^ where ZA are age classes of fledgling, sub-adult (second and third year) and adults, with breeding status B of Pre-Breeder, Breeder or Non-Breeder, which each have their own coefficients whereby 0 is the intercept, 1ZL is the beta coefficient for lay date and 2ZL^2 is the beta coefficients for lay date squared. The coefficients are shown in Appendix S1: Table S3 and S6.

SECTION S2 -GROWTH G(Z',Z)
Growth considers how all three states (ZA, ZL and ZB) develop from t to t + 1. Conditional on an individual surviving, an individual first ages (ZA). The lay date (ZL) would then develop for birds that were breeders, whilst non-breeders kept the same lay date. Finally, the breeding probability (ZB) was determined for the next breeding season.
Age has a simple directional model such that ′ , 1 where ZA age increases by one annually until a maximum age of 40 after which all remaining old-age birds die. Note that age in our model is in years from the census moment (May 1), as opposed to general bird ringing terminology that refers to calendar years (CY). So, for example a fledgling is age 0 (as opposed to 1CY) and enters the population at Age 1 ( Figure  1) when in bird ringing terminology it is 2CY.
We considered alternative model structures for the growth of lay date including whether lay date was a fixed trait (i.e. constant over time), if lay date at t + 1 (Z'L) depended on current lay date (ZL) and if it changed as an individual aged (ZA; Appendix S1). The parameters were determined through model-averaging and consequently included parabolic relationships with both age and lay date (Appendix S1). The growth models were fitted using generalised least squares which allows errors to be correlated or to have unequal variances. We considered different weighting structures to model heteroscedasticity in growth of lay date, with an exponential variance structure providing the best fit. The equation for lay date at time t + 1 (Z'L), which only applies to individuals that were a breeder at time t, (i.e. ZB = 1), was where 0 is the intercept, 1ZL is the slope with lay date at time t, 1ZL^2 is the slope with lay date squared, 3ZA is the slope with age and 4ZA^2 is the slope with lay date squared. Individuals that were non-breeders (i.e. ZB = 0) were multiplied with a diagonal matrix so that the lay date remained the same from t to t + 1. The coefficients for each parameter are shown in Appendix S1: Table S10. To be able to discretise the continuous state variable of lay date in the IPM, we identified the minimum and maximum lay dates to identify the boundary points, increased these by 10% and subsequently created 100 bins using the midpoint rule.
The final step in growth was to estimate breeding probability for t + 1 (ZB'). Prior research indicates that breeding status is not only age dependent, but it also depends on the previous breeding status, whereby breeders are more likely to remain breeders than nonbreeders becoming breeders (Ens et al. 1995). The model selection thus considered alternative model structures containing age and previous breeding status (Appendix S1). The model parameters included prior breeding status (ZB), a quadratic relationship with age (ZA) and an interaction between breeding status and age. The logistic equation for probability of breeding at time t + 1 (ZB'), with the constraint that ZA is three or higher, was , 1 1 where 0 is the intercept, 1ZB is the current breeding status with 0 = non-breeder and 1 = breeder, 2ZA is the slope with a continuous variable of age, 3ZA2 is the slope with age squared, 4ZBZA is the slope for the interaction between current breeding status and age and 5ZBZA2 is the slope for the interaction between current breeding status and age squared. Due to sample size restrictions, it was not possible to estimate separate, age-specific breeding probabilities for non-breeders and pre-breeders. However, since pre-breeders are by default a non-breeding class, the non-breeding breeding probabilities were used for the pre-breeding state as well. The structure of the IPM was such that it was only possible to move from Pre-Breeder to a Breeder, and once a Breeder, it was possible to alternate between Breeder and Non-Breeder (but impossible to return to Pre-Breeder). The parameter coefficients are shown in Appendix S1: Table S8. SECTION S3 -FECUNDITY R(Z) Fecundity, i.e. the number of offspring that were recruited at age one, took the general function of * * * * * * where NS, NH, HS, CL and FS are vital rates for nest success, number of hatchlings, hatchling survival, probability of initiating a replacement clutch and fledgling survival respectively, for the first (1) and second (2) clutches ( Figure 1; Appendix S1). The probability distribution of NH was estimated using a truncated Poisson distribution and took the general form of where P(X) is the probability of success for vital rates NS, HS, CL and FS (Figure 1), and x are the beta estimates for the variables x included in the vital rate function (Appendix S1).
We only considered whether reproduction vital rates were dependent upon the state variable of lay date (ZL), and not on age (ZA) nor breeding probability (ZB), since only breeders could reproduce (ZB = 1) and breeding probability was dependent upon age (ZA). Multiple regression analyses were performed to determine the relationship of reproductive vital rates with lay date (Appendix S1). We considered three model structures, namely an intercept model, lay date and a quadratic relationship with lay date, and model-averaged the coefficients from all considered models (Appendix S1). Almost all regression analyses for fecundity were performed using a single dataset containing 5,852 nest records (Appendix S1). The parameter coefficients for each of the reproductive stages are shown in Appendix S1: Table S12, S14, S16, S18, S20, S22, S24. The only vital rate for fecundity that was not estimated from this data was fledgling survival to age one. Similar to survival for adults, we performed a multi-state live and dead recoveries mark-recapture analysis to estimate survival to age one for individuals ringed during the fledgling phase (Appendix S1).

SECTION S4. DEVELOPMENT D(Z')
Oystercatchers that were recruited from R(Z,t) were automatically assigned an age (ZA = 1) and a pre-breeding state (ZB = 0). Although oystercatchers only begin breeding at a minimum