Unravelling processes between phenotypic plasticity and population dynamics in migratory birds

Populations can rapidly respond to environmental change via adaptive phenotypic plasticity, which can also modify interactions between individuals and their environment, affecting population dynamics. Bird migration is a highly plastic resource-tracking strategy in seasonal environments. However, the link between the population dynamics of migratory birds and migration strategy plasticity is not well understood. The quality of staging habitats affects individuals’ migration timing and energy budgets in the course of migration, and can consequently affect individuals’ breeding and overwintering performance, and impact population dynamics. Given staging habitats being lost in many parts of the world, our goal is to investigate responses of individual migration strategies and population dynamics in the face of loss of staging habitat, and to identify the key processes connecting them. We started by constructed and analysed a general full-annual-cycle individual-based model with a stylized migratory population to generate hypotheses on how changes in the size of staging habitat might drive changes in individual stopover duration and population dynamics, and to identify the key processes connecting them. Next, through the interrogation of census data, we tested these hypotheses by analysing population trends and stopover duration of migratory waterbirds experiencing loss of staging habitat. We found empirical support for our modelling-identified hypotheses: the loss of staging habitat generates plasticity in migration strategies, with individuals remaining on the staging habitat for longer to obtain food due to a reduction in per capita food availability. The subsequent increasing population density on the staging habitat has knock on effects on population dynamics in the breeding and overwintering stage. Our results demonstrate how environmental change that impacts one energetically costly life history stage in migratory birds can have population dynamics impacts across the entire annual cycle via phenotypic plasticity.


Introduction
their migration route, they must carefully manage the timing of departure and arrival 73 (Alerstam and Lindström 1990, Alerstam et al. 2003, Winkler et al. 2014). In general, 74 individuals that arrive at breeding grounds earlier have higher reproductive success 75 than those that arrive later (Marra et al. 1998, Norris et al. 2004, and selection 76 favours individuals that minimize the time spent travelling during the northward 77 migration (Lindstrom and Alerstam 1992). Migratory birds usually spend much longer 78 accumulating energy reserves in staging areas than in flying (Hedenström and 79 Alerstam 1997). Therefore, the total time spent on migration is consequently strongly 80 influenced by the quality of, and an individual's behaviour at, staging areas 81 (Hedenström andAlerstam 1997, Erni et al. 2002). where density-dependence is strongest (Sutherland 1996b). Habitat loss in staging 99 habitats can reduce food resources, decrease foraging and fat-accumulation rates of  In this study, we use theoretical modelling to generate hypotheses that we next 109 test with empirical data to examine the effects of habitat loss within staging areas on 110 individual migration strategies and population dynamics, and the process that connect dynamics in the staging, breeding, and overwintering grounds. We did this by 118 examined where in the life cycle density-dependence operated most strongly. Next, 119 using 13 years of census data, we examined whether observed empirical trends were 120 consistent with hypotheses we generated from our individual-based model. Our can be ignored such that we construct a female-only model.

139
Model description 140 Our IBMs includes three types of habitats in the model landscape, which are 141 wintering habitat, breeding habitat and staging habitat. The staging habitat split into 142 two habitat types -S1 and S2 area. In order to simulate the process of habitat loss in 143 the staging habitat, the size of S1 area remained constant across all simulations, while 144 the habitat size of S2 area was adjusted (Fig.S1a-c). The rest of the grid cells in the 145 landscape are non-habitat, where individuals pass by during migration but do not stop.

146
Each grid cell within those three types of habitats contained renewing food resources, 147 food resources renewed each time step at the habitat-specified food recovery rate after    Table S1 and Table S2 in Appendix1,   Table S1 in Appendix S1). More specifically: 197  Model 1 is the null model, with the same size and the food recovery rate 198 at all three habitats, and a "capital breeding" strategy for reproduction.  Appendix S1).

223
The average daily population density was the average of daily bird number in the 224 habitat during the period of the current life cycle stage each year, it was recorded at 225 the S1 area. The total number of individuals was recorded for the population and three 226 types of habitats each year respectively. Individuals were recorded in age classes and 227 reproduction status (juveniles, breeding adults, nonbreeding adult).

228
The per capita reproduction rate for adults, the reproduction rate among breeding

238
Mean individual stopover duration in the S1 area was recorded each year.   Statistical analyses for bird number trends 257 The abundance of all waterbirds, the abundance of the most common species, 258 and waterbird species richness were analyzed in this study. Survey sites that were 259 surveyed on less than 30% of survey dates were excluded from the analysis.  268 We calculated survey effort as the number of observers each day at each site as 269 Ej,k. We transformed date to Julian date t, and calculated t 2 as visual examination of 270 the data revealed a quadratic relationship. We treated 'year' (y) as a continuous 271 variable in our models, and we treated 'survey site' (k) as a categorical variable.
Statistical analyses for stopover duration trends 284 We estimated the stopover duration for each common species by using our 285 census data. Since the quadratic relationship between Julian date and bird abundance 286 revealed from our census data (Fig. S4 a), we first estimated normal distributions of 287 bird abundance within the period of stopover for each species each year (Fig. S4b). To 288 do this, we extracted the mean date and variance in date for each species each year 289 from census data, and scaled the curves by the bird number from the census data.

290
Then we estimated the date by which each quantile of the distribution of bird 291 abundance is reached, date at 2.5% quantile and 97.5% quantile was the date of 292 arrival at staging habitat and the date of departure from staging habitat for each 293 species each year within 95% confidential interval respectively (Fig. S4c). We 294 calculated the stopover duration based on the arrival date and departure date for each 295 species each year, denoting it as Ti,y (details are provided in S5.4 Appendix S5).

296
The distribution of Ti,y was normal distribution, so we fitted Linear model in 297 program R to test the stopover duration trends across years for each common species, The correlation between stopover duration and bird abundance 301 We fitted a regression between bird abundance and stopover duration to test the 302 correlation between them. We calculated the mean abundance during the stopover  309 We used the ANOVA command in R to assess the significance of each variable, 310 used adjusted-R 2 to assess the goodness of fit of the linear model, and used 1-311 (Residual Deviance/Null Deviance) to assess the goodness of fit of the GLMs. As the size of S2 area decreased through eight scenarios in the model 2 316 simulations, there was a decrease in the total number of individuals, and an increase in 317 average daily population density at the S1 area (Fig. 1a). The curve of daily 318 population density at the S1 area became steeper in the early and late phases of the 319 stopover period, and higher and wider when the population reached peak numbers 320 (Fig 1b), showing that the population took less time to reach the peak number and 321 remained at the peak number for a longer period of time. Individual's annual stopover 322 duration at the S1 area increased (Fig 1c), suggesting individuals stayed longer in the 323 S1 area as the size of S2 decreased. accumulation rate when birds remained in the staging habitat decreased (Fig 2a), the 335 distribution shifted towards lower energy reserves with fewer individuals having 336 reached the energy threshold for departure from the staging habitat (Fig 2b). The 337 energy reserves when individuals left the staging habitat was decreased, however, the 338 departure energy reserves from the breeding habitat and wintering habitat was 339 increased, both in adults and juveniles (Fig 2c).   results are listed, "density" is the average daily population density, "duration" is the 380 individual stopover duration, "energy" is the energy reserves when leaving the habitat,

381
"individual" is the number of individuals, "proportion*100" is the one hundred fold proportion of individuals of each age class in the population, "rate*100" is the one 383 hundred fold vital rates. Four types of habitat are listed, including the S1 area, the 384 staging habitat, the breeding habitat, and the wintering habitat. Age classes includes 385 juveniles and adults, adults were split into two types: breeding adults and non-   (Fig. S8b), species (F24,18373 = 80.14, p < 0.001) and the interaction between 418 year and species (F24,18328 = 6.96, p < 0.001). There was a quadratic association  Empirical data: stopover duration 436 Analysis of the stopover duration of common species revealed a significant increasing 437 trend at the rate of 1.39 days per year (F1,245 = 14.58, p < 0.001) (Fig.5d). The 438 stopover duration (F24,221 = 7.20, p < 0.001) and its temporal trends (F24,197 = 1.90, p = 439 0.009) were differed among species (Fig S8d). The adjust R-square of the linear 440 model was 0.428.

441
In the analysis of the correlation between bird abundance and the stopover duration, correlation between bird numbers and stopover duration (Fig. 1b, Fig. S9 in Appendix 495 S1). This contrasts the "buffer effect" hypothesis, which proposes that population 496 density changes only in space (Brown 1969, Sutherland 1996a, Gill et al. 2001). Our 497 results show that in addition to birds becoming more concentrated in the remaining 498 area, a decrease in the extent of the staging area can also result in an increase in time 499 spent there, leading to intensified competition during staging period. Therefore, 500 during the time-limited northward migration, the reduced staging area not only leads 501 to a higher population density in the staging habitat as individual birds stay longer, but 502 also that the high population density is maintained for a longer period of time during 503 this part of the life cycle (Fig.6a). And this can alter the strength of population 504 regulation in either the wintering or breeding areas (Fig. 6b). How do we reach this 505 conclusions?

506
In our models, individuals adjust their migration strategy to prolong the stopover 507 duration when food becomes scarce. As a consequence, individuals often arrive at the 508 breeding ground late and with low energy reserves, or they fail to reach the breeding 509 grounds. This results in fewer adults arriving in the breeding ground, fewer adults 510 reproducing, and a decline in total reproductive output. However, those individuals 511 that do arrive experience less competition, and have higher departure energy reserves 512 at the breeding habitat, breeding adults have higher per capita reproductive rate. All

521
All of this is due to increased competition in the staging habitat, which reduces 522 per capita food availability, and generates behavioural plasticity in migratory strategy.

523
As individual behaviour changes, so to do the population dynamics. The altered 524 population dynamics, in turn, affect individual behaviour (Miner et al. 2005). This 525 process continues until an equilibrium is reached, when the behaviour settles down to 526 equilibrium, as do the population dynamics (Fig. 6b). The connection between 527 individual phenotypic plasticity and population dynamics is the result of feedback 528 process across the annual cycle that can result in eco-evolutionary feedbacks are 529 argued by Coulson (2020). 530 If the part of the annual cycle that determines carrying capacity is not switched to 531 staging habitat, patterns and processes of individual strategies and population 532 dynamics can be different (Fig. 4). The staging habitat was no longer the stage