Warming springs and habitat alteration interact to impact timing of breeding and population dynamics in a migratory bird

In seasonal environments, increasing spring temperatures lead many taxa to advance the timing of reproduction. Species that do not may suffer lower fitness. We investigated why black‐tailed godwits (Limosa limosa limosa), a ground‐breeding agricultural grassland shorebird, have not advanced timing of reproduction during the last three decades in the face of climate change and human‐induced habitat degradation. We used data from an 11‐year field study to parameterize an Integral Projection Model to predict how spring temperature and habitat quality simultaneously influence the timing of reproduction and population dynamics. We found apparent selection for earlier laying, but not a correlation between the laying dates of parents and their offspring. Nevertheless, in warmer springs, laying dates of adults show a stronger positive correlation with laying date in previous springs than in cooler ones, and this leads us to predict a slight advance in the timing of reproduction if spring temperatures continue to increase. We also show that only in landscapes with low agricultural activity, the population can continue to act as a source. This study shows how climate change and declining habitat quality may enhance extinction risk.

, plankton and fish (e.g., Edwards & Richardson, 2004), any phenological mismatch may reduce fitness and population size. This reduction in fitness has been shown to be greater in habitats which have been significantly altered by humans in recent decades (Forister et al., 2010).

Modification of habitats often impacts survival and reproductive
rates (e.g., Allen et al., 2017;Piersma et al., 2016), for example, by altering food availability, introducing predators or making the habitat more suitable for predators. If reproduction becomes sufficiently suppressed that it cannot counter losses due to mortality, populations become "sinks" (Pulliam, 1988). Without a contribution from source areas, where more recruits are produced than adults die, such sink populations will decline in size and eventually go extinct. In addition, rapidly deteriorating habitats may impact environmental cues associated with optimally timed reproduction (Bourgault, Thomas, Perret, & Blondel, 2010;Winkler et al., 2014), making the combination of climate change and altered habitats especially difficult for populations to adapt to. Several integrative studies have disentangled climate or habitat effects on population dynamics (e.g., Gamelon et al., 2017;van der Meer, Jacquemyn, Carey, & Jongejans, 2016;Simmonds & Coulson, 2015) or have studied correlations between timing of reproduction and population growth (e.g., Both et al., 2006;Dunn & Møller, 2014). To our knowledge, the interactive mechanistic effects of habitat quality, climate change and timing of reproduction on population growth have not been examined before.
Agricultural landscapes have seen a higher rate of change in the past decades than any other type of landscape, and this change may well continue given projected human population growth (Foley et al., 2005;Tilman et al., 2001). On temperate agricultural grasslands, landscape-scale changes in agriculture, including increased fertilizer use and drainage of excess rain-and groundwater, combined with increasing spring temperatures have triggered earlier grass growth, insect emergence dates and advanced agricultural schedules (Kleijn et al., 2010).
In this study, we modelled the effects of spring temperature and agricultural land use intensity on laying dates while simultaneously examining how laying date affected fitness, and consequently population growth rates. We parameterised our structured model (Ellner & Rees, 2006;Smallegange & Coulson, 2013) using data from a long-term demographic study (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). We then used the model to (a) predict laying dates at different spring temperatures at different agricultural land use intensities, (b) explain why black-tailed godwits laying dates have not advanced in recent decades, (c) investigate the effects of spring temperature and agricultural land use intensity on population growth, and (d) predict population persistence and laying date in each of the different agricultural land use intensities if spring temperatures would continue to increase.

| Model framework: Integral Projection Model
We constructed an Integral Projection Model (IPM) to investigate how spring temperature and agricultural land use intensity influenced timing of reproduction and population dynamics of black-tailed godwits. An IPM is a discrete-time population projection model structured by a continuous trait (Coulson, 2012). Instead of body size, which is primarily used as continuous trait, we used laying date. We built the IPM from functions describing how laying date (z) in year t influenced: (a) adult survival to t + 1, S(z, t); (b) laying date in year t + 1, G(z′|z, t); (c) recruitment of offspring to the population at year t + 1, R(z, t); and (d) offspring laying date at t + 1, D(z′|z, t), where z′ denotes the laying date at time t + 1. These functions were parameterized for three different habitat types and three spring temperature scenarios using predictions from field data. The IPM can be written as: where n(z, t) is the distribution of character traits z at time t, and n(z ′, t + 1) is the distribution of character traits z′ at time t + 1. From the IPM, we calculated the long run population growth rate and mean laying date, and assessed the relative contribution of each parameter to laying date and growth rate with a sensitivity analysis. As laying date is a labile trait (Childs, Sheldon, & Rees, 2016), we also tested effects of spring temperature in t + 1 for functions G(z′|z, t) and D(z′|z, t).  . To enable early mowing, water tables are kept low, and the vegetation consists of a monoculture of fast-growing rye grasses Lolium sp., onto which fertilizers are applied Howison et al., 2018

| Spring temperatures
To define spring temperatures during field work years, we used data knmi.nl). We used R package "segmented" (Muggeo, 2008) to analyse the start of spring temperature rise in our study area ( Figure 1).

| Timing and demographic field data of godwits
Godwits breed on the ground; their nests were located and positions stored in a GPS. Because nests were almost always found during the incubation stage, we estimated laying date by floating the eggs in water and measured the float angle (Liebezeit et al., 2007). A nest was considered successful if at least one chick was found in the nest, or if we found an indication of successful hatching (Kentie et al., 2015). From 2008 onwards, 1-day-old chicks were marked with a flag-ring with a unique code. Adults were caught on the nest and uniquely marked with a flag and four plastic colour rings and a numbered metal ring. Chicks older than 10 days, including recaptured chicks wearing a flag-ring with a code, were also marked with a flag and four plastic colour rings and a numbered metal ring. Every year we scanned birds for colour rings on a daily basis throughout the breeding period. Colour marked black-tailed godwits were linked to the nest by observing returning parents from a distance, or by placing camera traps at the nest site. We collected a blood sample to genetically determine sex (see Trimbos et al., 2013 for details), but if a blood sample was lacking, we used morphological measurements to sex adults (Schroeder et al., 2008; 8% of 649 individuals).

| Statistical demographic analyses for G(z′|z, t)
Within-individual development of laying dates accounts for the flexibility of laying from 1 year to the next. Godwits can lay repeat clutches if the first failed . Including these repeat attempts would lead to a bias in the relationship between laying date at time t and at time t + 1.
However, usually we were unaware whether we found the first or a repeat nest. From laying dates of confirmed repeat nests (N = 67), we estimated that 95% of the repeat nests occurred after 27 April. in the year of breeding, we included Tsum at t + 1 as direct laying date cue instead of Tsum at t. We tested for effects of laying date at t, and Tsum at t + 1 and its interaction on laying date at t + 1, and included land use intensity as factor. We estimated the variance around the intercept to include in the IPM model.

| Statistical demographic analyses for R(z, t)
To estimate reproductive success, we separately estimated nest survival and first-year survival probability. We assumed 3.7 hatched eggs per successful nest (Kentie, 2015) and equal sex ratios at hatch We used binomial generalized linear models to test models for nest survival using all nests of which we knew when it was laid (N = 4,018, 951 and 966 nests for low, intermediate and high land use intensity, respectively). To account for nests, we have not found before they were lost, nest age when found was included as a covariate. A nest which is found shortly after it was initiated has a higher chance to get predated while being under observation than a nest which is due to hatch soon. We assumed a linear relationship between the covariate nest age when found and nest success, and kept it zero when extracting parameter estimates, so that nest survival was estimated from the moment of laying. We analysed nest survival with covariates laying date, quadratic laying date, land use intensity, Tsum and all two-way interactions.
We used Cormack-Jolly-Seber mark-recapture models with an age structure to estimate survival from hatch day until the following year. We always included land use intensity as predictor based on previous results (Kentie et al., , 2014, and we tested models with laying date of their birth nest, Tsum, and each interaction between laying date, land use intensity and Tsum as predictors for survival rate of the first age class (N = 2239 hatchlings). Survival probability after the first year was kept constant. Based on earlier results (Kentie et al., 2014), resighting probability always contained age class (first year and adult), whether a bird was ringed with a flag-ring or a more visible colour ring combination (3.6% of chicks) and was either constant or could vary between years. A preliminary test showed no evidence of a quadratic effect of laying date (dAIC > 44 than in the best model without quadratic effects), which we therefore not included in the models.

| Statistical demographic analyses for D(z′|z, t)
We tested for a parent-offspring association in laying dates with a linear model using a Gaussian error structure and estimated the variance around the intercept. Because most young godwits start breeding in their second year after hatching, we excluded earlier breeding attempts (N = 8) to exclude potential age effects. We included all chicks to increase the sample size (N females 19, N males 9, N unknown sex 8

| Statistical demographic analyses: model selection and inferences
Statistical analyses were carried out in R 3.4.3 (R Core Team, 2017).
We used RMark (Laake, 2013) for mark-recapture models. Goodness-of-fit was checked with program release from within R (release.gof). We used a variance inflation factor, or ĉ, of 1.5 for the data set with birds ringed as adults as the data were overdispersed ( | 5295 p > 0.5). Goodness-of-fit test showed no signs of overdispersion for the data set of birds marked as chicks (Test 2: χ 2 = 23.0, df = 3, p = 0.9). Model selection was based on Akaike's information criterion adjusted for small sample sizes (AICc; Burnham & Anderson, 2002).
Model selection procedures are described in the Supporting Information. We used Tsum/100 to reach model convergence.

| IPM modelling
As demographic relationships both depend on Tsum at time t and at time t + 1, we built stochastic IPMs with simulated yearly fluctuating Tsums (see Figure 2 for a life-cycle diagram). The IPM contained matrices of 250 discrete bins of laying dates, which ranged between −10 and 80. We built separate IPMs for fields with low, intermediate and high land use intensity, with past, present and future spring temperature scenarios. The three spring temperature scenarios had a mean Tsum of 350, 450 and 550°C, for "past," "present" and "future" respectively, and a standard deviation of 70, which corresponds to the standard deviation of the temperatures between 1900 and 1975 (see Supporting Information Figure S1 for simulated Tsums). Future spring scenario represents 2040 if Tsum will continue to increase at the same rate. The temperature scenarios were generated randomly for a 5,000-year time span. We then started with 100 individuals laying at 24 April, and followed the population trajectories for each. We inspected the results, then discarded the first 100 time steps and stored laying date, population structure and growth rate (natural log lambda) of the remaining 4,900 time steps. Sensitivities of the vital rates were tested for the models with a current temperature scenario, by subsequently changing each parameter by 0.1%.

| Are there population level changes in laying dates?
The annual mode of laying date, which corresponds to peak laying dates (and thus independent of repeat clutches), did not advance

| Demographic and laying date functions
Adult females survived with a probability of 0.86 (0.84-0.87 95% CI) to the next year (  (Figure 4).
Laying dates of individual females were more repeatable between consecutive years in warm springs at t + 1 (i.e., had a slope closer to x = y) and were not correlated with land use intensity (Table 1d, Figure 5a). Laying date of young godwits 2 years after hatching did not correlate with Tsum when hatched, Tsum 2 years after hatching, nor laying dates of their parents (Table 1e; Figure 5b).
We found this lack of relationship after removing a suspicious outlier, which was a nest with a laying date of 25 May. This appeared to be a second breeding attempt, as this individual was seen nest building on 11 April.

| IPM: effect of habitat, spring temperature and laying date on population growth rate and laying date
The results of the IPMs predict that mean laying dates did not differ much between different land use intensities and showed a wide standard deviation of the mean (Figure 6a). The predicted laying dates for past and present spring temperature scenarios are close to the observed peak laying date; a difference of 1-2 days (see also

| Sensitivity analysis
The sensitivity analysis showed that population growth rates were Previous studies have shown that black-tailed godwits have not advanced their laying dates since the early 1980s (Kleijn et al., 2010;Schroeder et al., 2012) or possibly even the 1930s (Meltofte et al., 2018), which contrasts most temperate species that advanced their reproductive timing in response to increasing spring temperatures (Crick & Sparks, 1999;Parmesan, 2007). This is surprising as we showed that godwits would have higher fitness if they bred earlier.
Are they unable to respond to a warming world? Over the past decade, we did not observe advancement of laying dates. However, we did find that the yearly mode of laying dates (the mean or median could be biased by repeated nesting attempts in years with low nest survival) negatively correlated with spring temperatures, although  We predict that, if springs become warmer, mean laying dates will advance by 4-6 days-not enough to match the optimal timing of offspring production. As we did not find a correlation between laying dates of the parents and laying dates of their offspring, the predicted shift in laying dates was entirely caused by phenotypic flexibility of the adults. The low number of recruits prevented us from carrying out more formal tests for heritability of laying dates (Charmantier & Gienapp, 2014).
The subspecies of black-tailed godwits breeding on Iceland (Limosa l. islandica) also faces warming springs and agricultural intensification (Gill et al., 2007;Jóhannesdóttir, Alves, Gill, & Gunnarsson, in press), but spring warming at this more northern latitude and the relative low level of agricultural intensification so far seems to have benefitted the birds. The Icelandic black-tailed godwits seem to be more flexible in general; laying dates within pairs were not repeatable and advanced by 2.5 days per year (and 4 day/°C June temperature) during the last 11-year period, but little is known about how individual arrival affects laying date in this population. Furthermore, Icelandic godwits were able to advance arrival date by 2 weeks over the last 20 years (Gill et al., 2014), whilst the continental West-European godwit population did not advance arrival dates since at least 1970 (Kleijn et al., 2010).
To explain these timing phenomena, Gill et al. (2014) suggested a mechanism that does not require phenotypic plasticity or evolution of phenology: Early hatched offspring would winter at higher quality wintering grounds, which would allow for earlier migration to the breeding grounds. Although wintering location also affects arrival and laying date in Dutch godwits-those who winter at the most southern location arrive and lay earlier than those wintering at northern sites   (Senner, Conklin, & Piersma, 2015).
Mechanisms that influence phenology and habitat choice could thus differ in closely related subspecies. It is also possible that over time, mechanisms influencing timing of reproduction will change, such as shown with an increase in heritability of laying dates in great tits (Parus major; Husby, Visser, & Kruuk, 2011). Our sensitivity analysis showed that changes in the flexibility of adult laying dates have the relative strongest effect on population laying dates and growth rates, which makes this a potential mechanism to select upon. Showing their potential to react to temperature, we observed that during a cold spell in the Netherlands in early spring 2013, black-tailed godwits returning from West-African wintering sites were able to respond by delaying their arrival and laying dates (Senner, Verhoeven, Abad-Gómez, et al., 2015). Yet, perhaps a cold spring causes a direct behavioural response due to an immediate lack of food, while they are less sensitive to slower ecological processes such as a mismatch between reproduction and food peaks of which the outcomes -lower reproductive success-are less directly experienced.

| Population consequences: sources and sinks
The Dutch godwit population started to decline in the same period when spring temperatures started to increase and agricultural intensification accelerated (Harms, Stortelder, & Vos, 1987;Kentie et al., 2016 (Howison et al., 2018), albeit in lower densities. Young godwits disperse larger distances than adults and seem to be less discerning about habitat quality (Kentie et al., 2014), suggesting that intermediate and intensive agricultural areas are sinks, which are maintained by godwits breeding on fields with low land use intensities which act as sources (Pulliam, 1988).
The effects of spring temperature and land use intensity on reproductive output is a combination of mowing dates and food availability for the chicks. On intensive agricultural fields, warmer springs enhance grass growth and mowing starts earlier in the season than in cold springs. Mowing destroys the nests, or, if marked by local volunteers so that farmers can spare them, removes cover and increases depredation probability . Mowing also poses a threat for prefledged chicks, which may partly explain the decrease in first-year survival during the season, as most mortality happens in the first week after hatching (Schekkerman, Teunissen, & Oosterveld, 2009). Additionally, late hatched chicks have a lower body condition (Loonstra, Verhoeven, & Piersma, 2018), possibly caused by a mismatch relative to the insect food peak, or due to a decline in insect biomass caused by mowing (Schekkerman & Beintema, 2007). Godwit families with prefledged chicks are mobile and, if in the vicinity, families from mown fields can move to grassland with delayed mowing regimes. Most grassland with low land use intensity is especially managed to protect breeding birds; here agricultural activity is delayed until at least 15 June regardless of spring temperature, that is, after the nesting phase .
Yet, on fields where mowing is postponed, chicks hatched from nests laid after 26 April (assuming 25 incubation days and 25 prefledging days) also have a high mortality risk due to agricultural activities.
There are only few studies that identify the possible mechanisms that underlie (the lack of) phenological change and integrate this with population dynamics (but see Vedder, Bouwhuis, & Sheldon, 2013;Weegman, Arnold, Dawson, Winkler, & Clark, 2017). Our model suggests that for black-tailed godwits, climate change has considerably stronger negative effects in areas with most human-induced habitat alterations, compared to habitats which are less affected by humans. This stronger effect is caused by demographic differences between these habitats rather than phenological differences. Habitat degradation is a much greater risk, and for blacktailed godwits, as dependent as they are on agricultural landscapes not only during reproduction but also in the rest of the year, a serious threat. Many other species also live in habitats both affected by habitat and climate change, while when modelling extinction risks often only climate change are taken into account (Urban, 2015).
Although species may adapt to climate change, either through phenotypic plasticity or evolution, large-scale human-induced habitat change may accelerate species extinctions even more. It is therefore important to preserve landscapes in which species are able to keep up with global climate change.

ACKNOWLEDGEMENTS
We are grateful to the many farmers, most of them organized in the "Collectief Súdwestkust," and the conservation management organizations "It Fryske Gea" and "Staatsbosbeheer" for cooperation and granting us access to their properties. We thank our field crews and students from 2007 onwards and volunteers of local bird conservation communities for locating nests, and the many volunteers who Welfare Act Articles 9, 10 and 11 of animal experiment documents.

CONFLI CT OF INTEREST
We have no competing interests.