Using stationary vital rates in impact assessments may underestimate potential threat

Population viability analysis (PVA) is commonly used to assess future potential risks to threatened species. These models are typically based on mean vital rates, such as survival and fecundity, with some level of environmental stochasticity. However, the vital rates of wild populations, especially those already exhibiting declining trajectories, may be nonstationary, such that the mean or variance changes over time. In this study, we examined whether including observed temporal trends in vital rates affects the predictive accuracy of PVA, as well as the projected impact associated with a hypothetical threat. To achieve this, we ran a series of simulations using Leslie matrix PVA models that included different combinations of environmental stochasticity, temporal trends in vital rates, and threat. We apply our analysis to a long-lived colonial species of seabird, the black-legged kittiwake Rissa tridactyla, that is classified as globally Vulnerable and is potentially highly sensitive to offshore renewable energy development. We found that including observed temporal trends in vital rates was (i) crucial for the accurate reconstruction of observed population dynamics and (ii) had a dramatic effect on the projected impact from the hypothetical threat. In an era when many animal and plant populations are declining due to long-term trends in their vital rates, we identify that including this demographic structure is essential for robustly evaluating potential threats using PVA models. Omitting observed temporal trends in vital rates from impact assessments is highly likely to yield unreliable results that could misinform conservation and management decision making. This result has immediate application for conducting impact assessments on protected species and populations.

Population viability analysis (PVA) is commonly used to assess future potential risks to threatened 23 species. These models are typically based on mean vital rates, such as survival and fecundity, with 24 some level of environmental stochasticity. However, the vital rates of wild populations, especially those 25 already exhibiting declining trajectories, may be nonstationary, such that the mean or variance changes 26 over time. In this study, we examined whether including observed temporal trends in vital rates affects 27 the predictive accuracy of PVA, as well as the projected impact associated with a hypothetical threat. 28 To achieve this, we ran a series of simulations using Leslie matrix PVA models that included different 29 combinations of environmental stochasticity, temporal trends in vital rates, and threat. We apply our 30 analysis to a long-lived colonial species of seabird, the black-legged kittiwake Rissa tridactyla, that is 31 classified as globally Vulnerable and is potentially highly sensitive to offshore renewable energy 32 development. We found that including observed temporal trends in vital rates was (i) crucial for the 33

Introduction 41
In many countries, it is necessary to conduct a full evaluation of potential impacts to protected wildlife 42 populations to gain consent for certain activities, such as biological resource use, land-use change or 43 energy production (e.g., EU Birds Directive 2009/147/EC, EU Habitat Directive 92/43/EEC, EU EIA 44 Directive 85/337/EEC). A commonly used approach for inferring how populations may respond to these 45 activities is to predict how vital ratessuch as survival or fecundityare likely to be impacted, and 46 then use these altered values to evaluate the population response. Population viability analysis (PVA) 47 4 including this structure in impact assessments is essential for obtaining meaningful estimates of 70 potential threats. 71 72

Methods 73
Prevalence of temporal trends in vital rates 74 To demonstrate the potential prevalence of temporal trends in the vital rates of seabirds, we collated 75 colony-specific fecundity data for kittiwakes breeding in the UK. Here, fecundity rates are defined as 76 the proportion of apparent nesting attempts that produce a successful fledgling each year. Data were 77 extracted from the Seabird Monitoring Programme Database (SMP 2020), and we included colonies 78 with at least eight years of data between 1985 and 2019. To identify colonies with a significant temporal 79 trend in rates of fecundity, we used Poisson generalised linear models fitted to each colony with a log 80 link function. In these models, the response variable was set as the number of chicks alive at fledging 81 for each colony, and the number of nests was included on the log scale as an offset. This analysis was 82 conducted in Program R (v. 4.0.2) (R Core Team 2020). To examine whether kittiwakes demonstrate a 83 temporal trend in rates of fecundity at the national level, we fitted a Poisson generalised linear mixed 84 model to the complete dataset using the R statistical package "lme4" (Bates et al. 2016). 85

86
How do temporal trends in vital rates influence predicted population dynamics? 87 We ran five PVA simulations to examine how observed temporal trends in vital rates may influence the 88 predictive accuracy of PVA, as well as the projected impact associated with a hypothetical threat (Table  89 1, code provided as supplementary information). The first PVA simulation was deterministic, based on 90 mean demographic rates only. The second model incorporated environmental stochasticity by 91 assigning demographic rates from a beta distribution for survival and a gamma distribution for fecundity. 92 To set the expectation and variance of these distributions, we used colony specific values for kittiwakes 93 breeding on Skomer Island, Wales (51.74°N, 5.30°W). The vital rates for this colony, estimated for the 94 period 1989 to 2018, are 0.86 (±0.06 SD) for adult survival ( t  ) and 0.68 (±0.08 SD) for fecundity ( 5% annual mortality to each age class to simulate the hypothetical threat. The fourth model also 97 included environmental stochasticity but added a negative temporal trend on fecundity corresponding Finally, the fifth model included environmental stochasticity, the negative temporal trend on fecundity, 100 and the hypothetical threat, i.e., an additional 5% annual mortality applied to each age class (Table 1). 101 For each PVA simulation, we used a Leslie matrix approach following a pre-breeding census to predict survival rates (i.e., from age 1 to 2 years) to be 87.5% of adults, and survival rates to be comparable to 113 adults from age 2 onwards (Cam et al. 2005). Therefore, we assumed an additive relationship between Here, birds survive the juvenile age class (0-1 years) with probability 2 t t B S if the mother was age 3 121 years and 2 t t f S if the mother was age 4 or older, the immature age class (1-2 years) with probability 122 t I and one pre-breeding year (2-3 years) as well as two breeding age classes (3-4 years and all birds 123 4 years or older) with a survival rate equivalent to that of adults ( t  ).

125
We ran the five PVA simulations for 1000 iterations. Within each iteration and across all five PVA 126 simulations, temporal residuals associated with environmental stochasticity were held constant (Table  127 1) by using the same draws from the beta distribution for survival and gamma distribution for fecundity. To examine predictive accuracy across the five PVA simulations (Table 1), we compared the simulated 141 trajectories to the observed count data for kittiwakes breeding on Skomer Island. To evaluate how 142 observed temporal trends in vital rates may influence the simulated impact associated with a 143 hypothetical threat, we compared the output of the models with and without this demographic structure, 144 as well as the additional 5% mortality, i.e., Model 2 vs. Model 3, and Model 4 vs. Model 5 (Table 1). 145 146

How do temporal trends in vital rates influence predicted population dynamics? 159
We ran five PVA simulations (Table 1) to examine how incorporating observed temporal trends in vital 160 rates may influence predictive accuracy, as well as the projected impact associated with a hypothetical 161 threat. The population trajectory predicted by the deterministic model (Model 1, Table1) declined slightly 162 but not at the rate observed in the colony ( Fig. 2A). Likewise, the mean population trajectory predicted 163 by the stochastic model (Model 2) overestimated colony size, such that the lower limit of the 95% 164 confidence interval was above the observed trajectory ( Fig. 2A). By contrast, the population trajectories 165 simulated by the model including the observed decline in rates of fecundity (Model 4) were much closer 166 to the observed count data (Fig. 2B). As expected, including an additional 5% annual mortality across 167 all age classes increased the rate of population decline (Fig 2). Combining this hypothetical threat with 168 the observed trend in fecundity also meant that local extinction was highly likely to occur within the 30 169 year timeframe of the simulation (Fig. 2B).  Fig. 1). However, this demographic structure is often 175 overlooked in impact assessments that evaluate potential risks to threatened populations (Chaudhary 8 To recreate the observed population trajectory for kittiwakes breeding on Skomer Island, it was 181 necessary to include a negative trend in rates of fecundity. The deterministic model failed to replicate 182 the observed rate of decline, as did the stochastic model with a stable mean rate of fecundity. This 183 result emphasises that an essential prerequisite step to using PVA for impact assessment should be 184 testing the predictive accuracy of models by simulating dynamics for a period with available count data. 185 For populations where vital rate data are available, a mismatch between the simulated and observed 186 trajectories will highlight that the distribution of vital rates (mean or variance) has changed over the 187 study period, and that this demographic structure needs to be incorporated to support reliable 188 assessments of impact. We show that temporal trends in rates of fecundity are common and widespread for kittiwakes breeding 201 in the UK (Fig. 1). A recent meta-analysis demonstrates that this trend is reflected in fish-eating, surface-202 foraging species of seabird, like kittiwakes, across the northern hemisphere (Sydeman et al. 2021). In 203 our case study colony, fecundity is declining at a shallower rate than observed in many of the other UK 204 colonies of kittiwakes. For the colonies that demonstrate a greater rate of change in fecundity, omitting 205 temporal trends from impact assessments is likely to have a greater effect than demonstrated here. 206 However, not all populations will have long-term monitoring data available. When conducting impact 207 assessments on data-limited populations, we recommend an appropriate monitoring study be 208 undertaken to quantify the local population trajectory and stability of demographic rates. If this is not 209 possible, an alternative approach is to reconstruct mean vital rates using predictive methods based on Consequently, there is a pressing need for more research determining the prevalence of density-224 dependence feedbacks across these species, how best to incorporate these structures into impact 225 assessment models and the sensitivity of results to their inclusion. 226

227
In this paper, we demonstrate that temporal trends in rates of fecundity reflecting climate change and 228 habitat degradation are relatively widespread for populations of kittiwakes breeding in the UK. We also 229 show that omitting this structure from PVA models used to assess potential impacts to protected 230 populations can dramatically reduce predictive accuracy and underestimate projected risk. We 231 advocate that a prerequisite step to conducting impact assessments with PVA models involves testing 232 predictive accuracy by simulating population dynamics for a period with available count data.  vital rates without additional mortality, red bold line shows mean trajectory and red shading shows 95% 348 confidence interval) and Model 3 (environmentally stochastic vital rates with additional 5% mortality, 349 blue bold line shows mean trajectory and blue shading shows 95% confidence interval). Panel B shows 350 Model 4 (environmentally stochastic vital rates and trended fecundity without additional mortality, red 351 bold line shows mean trajectory and red shading shows 95% confidence interval) and Model 5 352 (environmentally stochastic vital rates, trended fecundity and additional 5% mortality, blue bold line 353 shows mean trajectory and blue shading shows 95% confidence interval). Observed population 354 trajectory for kittiwakes at Skomer shown as grey dashed line in both panels.