Fuelling conditions at staging sites can mitigate Arctic warming effects in a migratory bird

Under climate warming, migratory birds should align reproduction dates with advancing plant and arthropod phenology. To arrive on the breeding grounds earlier, migrants may speed up spring migration by curtailing the time spent en route, possibly at the cost of decreased survival rates. Based on a decades-long series of observations along an entire flyway, we show that when refuelling time is limited, variation in food abundance in the spring staging area affects fitness. Bar-tailed godwits migrating from West Africa to the Siberian Arctic reduce refuelling time at their European staging site and thus maintain a close match between breeding and tundra phenology. Annual survival probability decreases with shorter refuelling times, but correlates positively with refuelling rate, which in turn is correlated with food abundance in the staging area. This chain of effects implies that conditions in the temperate zone determine the ability of godwits to cope with climate-related changes in the Arctic.

G lobal warming is not globally uniform. In the Arctic, climatic changes are the strongest 1,2 , with the highest rates of advance in spring phenology, which contrasts with the slower changes in equatorial regions. Migratory bird populations that fail to maintain a match between the timing of breeding and the local phenology of resources have low reproductive output and show declines 3,4 . Some populations appear capable of tracking spring at the breeding grounds, but little is known about the fitness costs that such adjustments imply. We may expect these costs to be especially high in long-distance migrants whose annual cycles include use of widely separated places with different rates of phenological drift.
To examine the fitness trade-off involved in the adjustment of reproductive timing by Arctic breeding birds to advancements in the general phenology of the tundra, we studied bar-tailed godwits (Limosa lapponica taymyrensis; hereafter godwits). This population is among the several long-distance migratory shorebirds that travel from wintering grounds in West Africa to breed in the Russian Arctic with a single refuelling stop (of ca. 25 days) in the Wadden Sea of north-western Europe 5,6 . Godwits arrive in the intertidal areas of the Wadden Sea in late April-early May after a non-stop flight of about 5000 km, and spend most of their available time feeding 7,8 . The birds double their body mass in less than a month, which fuels their next 5000 km long migratory flight to the Arctic breeding grounds 9 .
Our study combines data on godwits and their food resources at the main wintering, spring refuelling, and breeding sites along the flyway. To track in detail how individual birds connect these sites, we instrumented eight godwits with satellite transmitters in 2016. To estimate population trends, we counted birds each winter from 2002 to 2016 at the main wintering area, the Banc d'Arguin, Mauritania, West Africa. For godwits staging in the Wadden Sea during northward migration, we assessed the relationship between the annual refuelling rates and the density of their main prey, adult lugworms (Arenicola marina), using a 21-year dataset and a hierarchical Bayesian model that accounted for year-specific arrival dates and arrival mass of godwits.
Godwits are sexually dimorphic, with females being 17% larger than males 6 . For this reason, we modelled arrival mass and refuelling rates at the Wadden Sea as being sex dependent. The departure dates from the Wadden Sea were derived by subtracting duration of migration (obtained in 2016 with satellite transmitters) from the observed arrival dates on the breeding grounds at the Taimyr Peninsula, Russian Arctic (Supplementary Table 1). For each sex we then estimated the effects of refuelling time and refuelling rate on subsequent survival probability of individually colour-marked birds. To evaluate the effects of climate change on the breeding ecology of godwits, we assessed the influence of dates of snowmelt 10 and emergence of the main food of shorebird chicks (adult crane flies; Tipula sp. 11,12 ) on arrival and breeding dates of godwits over a 25-year period (1992-2016).
The emergence of crane flies in the Arctic is advancing in concert with advancements of snowmelt. The godwits follow these changes and arrive on the breeding grounds earlier and initiate nests earlier too. They manage this by shortening the time for refuelling in the Wadden Sea, but this comes at the cost of lower subsequent survival probability as the shorter time is not fully compensated by higher refuelling rates. To some extent higher food densities in the Wadden Sea would increase the level at which Arctic warming can be mitigated on temperate shores.

Results
Spring migration and breeding phenology of godwits. The date of snowmelt at the godwit breeding grounds on Taimyr has shifted forward by 0.73 ± 0.16 s.e.m. d year −1 (P < 0.001, Table 1, Fig. 1a). The main prey for godwit chicks, the adult Tipula crane fly, has responded to changes in snowmelt by advancing its first emergence time by 0.38 ± 0.14 (P = 0.018) days per day of advancing snowmelt (Table 1, row 2, Fig. 1b). The godwits advanced their arrival on the breeding grounds by 0.22 ± 0.071 days (P = 0.006, Table 1, row 3) and advanced clutch initiation by 0.56 ± 0.17 days (P = 0.006, Table 1, row 4) per day of snowmelt advancement. The path analysis of the statistical causality in the phenological data (Fig. 1b and Supplementary  Table 2), revealed that breeding date of godwits was driven by arrival date, and not by the timing of crane fly appearance, while godwit arrival on the breeding grounds, as well as crane fly appearance, was driven by timing of snowmelt.
Satellite tracking of godwits confirmed the previously suggested migration strategy 6 of swift and direct migrations between West Africa and the Wadden Sea and between the Wadden Sea and Taimyr (Fig. 1c). Individual birds stopped only on six occasions during ten migration legs (two from West Africa to Wadden Sea and eight from Wadden Sea to Taimyr) for a mean period of 1.8 d (maximum 2.8 d) and completed the 5000 km migration leg between the Wadden Sea and Taimyr in ca. 5 Annual survival probability of godwits. The fuel load accumulated in the Wadden Sea is the product of refuelling time and refuelling rate 13 . This relationship means that birds can reach the same level of body stores in a shorter time by increasing their refuelling rates. However, in association with the reduction in refuelling time, during the study period annual survival probability of godwits decreased by 2% (−0.08 ± 0.01 on logit scale, ΔQAICc = 35.11, Table 1 Fig. 2e and g). Even though lugworm densities did not increase statistically significantly during the study period (0.15 ± 0.14, P = 0.29, Table 1, row 17), the refuelling rates of godwits did (0.07 ± 0.02 g d −1 year −1 in both sexes, P = 0.006, Table 1, row 18, Fig. 2d and f). The increased refuelling rate partly offset the effect of reduced refuelling time on survival probability. In males the annual survival probability in 2015 was 3% higher than expected if refuelling rates would not have changed since 1995. In females, the compensation in survival was 2%. The decrease in annual survival probability caused by changes during refuelling in the Wadden Sea would thus have contributed to the 4.0 ± 1.6% per year population decline (P(λ < 1) = 0.994, Table 1, row 20) revealed by midwinter counts on the Mauritanian wintering grounds (Fig. 1d).
Thus, to maintain the annual survival probability as initially observed, godwits would need to refuel at higher rates than measured in our study, and therefore would require even higher densities of lugworms. For full compensation (i.e. a refuelling rate of 6.6 g d −1 ), males would need an average density of 20 lugworms m −2 . These densities occur in the Wadden Sea occasionally. Females, however, for full compensation would need to refuel at 9.9 g d −1 , which would require 32.7 lugworms m −2 , an average density never encountered by our monitoring effort in the Wadden Sea during the study period (Fig. 2h).

Discussion
Our long-term, hemispheric scale observations suggest an important and previously unrecognized mechanism by which migratory birds cope with global change. Rather than the use of multiple sites simply being a liability 3,14 , it may provide opportunities for among-season compensation 15 . In contrast to many other species [16][17][18] , bar-tailed godwits adjusted their arrival on the breeding grounds and the onset of breeding, thereby tracking the seasonal advancement of their main arthropod prey on the breeding grounds. They achieved this by shortening their refuelling period in the Wadden Sea, albeit at the cost of lower survival especially in years of low lugworm densities.
Even though godwits were able to compensate for the reduced refuelling time by increasing refuelling rates, these rates were insufficient in years when lugworm abundance was low. In such years, godwits left the Wadden Sea with lower body stores that compromised their subsequent survival probability. Females were more sensitive to the shorter refuelling times than males, perhaps because they are larger and have higher energetic requirements such as the need to produce eggs after arrival on the tundra breeding grounds 19 . According to our calculations, refuelling females need lugworm densities 2.2 times higher than the average observed in the Wadden Sea over the last two decades to fully compensate for the observed shorter staging duration (Fig. 2h).  Refuelling rate required to maintain annual survival rates at 1995 in 2015 5.6 g d - 1 4.9 g d -1 9.9 g d -1 6.6 g d -1 9.9 g d -1 9.9 g d -1 6.6 g d -1 6.6 g d -1 1995 =0.93 1995 =0.89 The refuelling conditions in the Wadden Sea are critical for godwits to cope with earlier snowmelt on the Arctic breeding grounds. Even though refuelling rates are determined not only by food availability, but may be limited physiologically 20,21 and by other environmental factors, such as disturbance rates 22 , improvement of food stocks at staging sites can increase survival rates in godwits. Thus, to mitigate negative climate-change effects on godwit population travelling through the Wadden Sea, we need to maintain fuelling conditions for them. As a first step we may suggest the suspension of mechanical lugworm harvesting practices in the Dutch Wadden Sea 23 .
The population of bar-tailed godwits we studied is just one of many long-distance migrant birds challenged by the rapid global warming of the Arctic 24 . Food-related limitations on refuelling rate are likely to be a common problem in populations that need to keep up with advancing springs. It is a sobering realization that such refuelling areas are poorly protected in some parts of the world 25 with many being entirely lost or reduced in quality by urban and industrial developments [26][27][28] . Proactive, international collaborations focused on maximizing resources at staging sites for Arctic breeding migratory birds could help maintain the connectivity between the worlds' variably changing biomes.

Methods
Snowmelt dates from remote sensing data. We used the NOAA Climate Data Record (CDR) estimates of extent of snow cover on a 100 × 100 km grid based on remote sensing data 10 . Weekly snow cover for 1992-2017 was overlain with the breeding range of taymyrensis bar-tailed godwit 29 . For each grid cell snowmelt date was estimated as the next day after continuous snow cover period. Annual snowmelt date for the breeding range was estimated as the mean date across all grid cells and annual snowmelt date at the breeding site (South-Eastern Taimyr) was estimated as the mean of the two closest NOAA grid cells. Because rangewide trend in snowmelt dates did not significantly differ from the trend at the field site (the slope of the linear regression of overall snowmelt date on that at South-Eastern Taimyr did not differ significantly from unity 1.00 ± 0.76, t = 1.32, d.f. = 35, P = 0.20) we used the snowmelt dates at the site in the analyses below.
NOAA CDR data have estimates of uncertainty within 3-5%, but since snow cover was also recorded daily at the field site 30 (72.8°N, 106.0°E, Supplementary Table 1), we checked whether remote sensing data matched our field observations. The slope of the linear regression for the two NOAA grid cells closest to the field site on the snowmelt dates derived from field measurements did not differ significantly from unity (0.97 ± 0.129, t = −0.2, d.f. = 16, P = 0.84), and the intercept (0.18 ± 1.05) did not differ significantly from 0 (t = 0.17, d.f. = 16, P = 0.87). Details on snow data processing and analysis are available as Supporting information 31 .  8°N, 106.0°E). A three square km area was systematically searched for nests in each year. Clutch initiation dates were determined either by egg flotation 32 or back-calculation from recorded hatching dates. Clutch initiation dates for all 13 nests found and all phenological data collected at the Taimyr are presented in Supplementary Table 1. Phenological trends over time were determined as slopes of linear regressions of the observed dates over years. All data and analysis code are available as Supporting information 31 .
Structural equations modelling of the phenology data. To estimate the statistical causality among the observed phenological variables we used path analysis, a special case of structural equations modelling framework 33,34 . In the proposed model we estimated the strengths of the potential causal relationships between phenological variables (dates of arrival to the Wadden Sea, arrival to Taimyr, clutch initiation, snowmelt, and crane fly emergence). The model contained only one independent variable, time (year). The model structure is presented in Fig. 1b. Variables were assumed to have latent state and variable-specific normally distributed errors: Arrival Taimyr Crane fly emergence Clutch initiation Arrival dates in the Wadden Sea were hypothesised to affect arrival dates to Taimyr, and snowmelt dates on Taimyr to affect all phenological events except for arrival dates in the Wadden Sea. Clutch initiation dates were hypothesised to depend on all variables except arrivals in the Wadden Sea. We used mean estimates for dates of arrivals to Taimyr and snowmelt without uncertainties to ease their combination with other phenological observations. All the dates were centred to have zero mean but not scaled. Effects of variables were estimated as maximum probability of parameter to be strictly positive or strictly negative. The model parameters were estimated with MCMC JAGS sampler 35 via the R2jags 36 interface from the R computing environment 37 , the model data and code are available in supporting information.
Spring migration parameters from satellite tracking data. In 2016 we deployed 5-g solar-powered Argos PTT-100 transmitters (Microwave Telemetry Inc., Maryland, USA) on eight female godwits (two on the wintering grounds in Mauritania and six during the refuelling period in the Wadden Sea). We used legloop harnesses weighing ca. 1 g, similar to the ones used successfully for blacktailed godwits (Limosa limosa) 38 . The total attachment mass was ca. 1.2% of the body mass for females departing from the Wadden Sea and 1.5% for females departing from Mauritania 6 . We tracked birds via the Argos system (CLS France, http://www.argos-system.org/) and removed occasional outliers in the Argos data with a hybrid filter 39 . We estimated migration duration as total time that it took a bird to get from wintering to refuelling sites and from refuelling sites to the breeding grounds. Stopovers were defined as time periods during which the birds were located within a 25-km radius for at least 24 h. The northward migratory tracks of the eight godwits are shown in Fig. 1c.
Arrival dates to the Wadden Sea using citizen science data. The arrival of godwits to the Dutch coastline has been monitored by a citizen science project http://www.trektellen.nl from 1992 to 2016. In this project, experienced observers count migrating birds at established sites along the Dutch coast 40 . To estimate mean godwit arrival dates, we used observations between 10 April and 25 May (the period when godwits are known to arrive 8 ) from sites that had over 100 records of at least a single godwit (n = 7). The final dataset contained almost 400,000 godwits recorded during 2318 counting sessions.
To estimate annual mean arrival date to the Wadden Sea T0 k , we extended the model by Lindén and Mäntyniemi 41 . We assumed that at year k true daily number of arriving godwits N ijk is distributed normally over days with mean arrival date T0 k and standard deviation σT0 k . Each observation site i has its own multiplicative effect on the number of godwits that does not shift T0 k . Counts at the sites are proportional to the daily arrival but are outcomes from a random observation process with negative binomial error.
The expected number of birds N ijk on a day j at a site i and a year k can be calculated then with the following equation: where T ijk is the observation day; offset ijk is the natural log of observation duration (in hours); a 0 is the overall baselinenatural log of average number of birds passing through an average site in average year; total k is the annual random effect; site i is the random effect for observation site.
The observed count Y ijk was assumed to be a result of a random observation process with the error following the negative binomial distribution.
Godwits arrive in small flocks 42 and, therefore, the negative binomial distribution of counts was preferred over Poisson or quasi-Poisson distributions 41 . The model was estimated with JAGS sampler 35 via the R2jags 36 interface from the R computing environment 37 and reasonably fitted the data 43 (Bayesian P-value 0.59).
Data collection on food abundance in the Wadden Sea. We obtained lugworm (Arenicola marina) densities from 1996 to 2016 on the basis of the biannual sampling effort at 15 permanent sampling stations located at Balgzand in the western Dutch Wadden Sea 44,45 (52.9°N, 4.8°E). We used the densities of adult lugworms in late winter (Feb-Mar, mean of all stations) as a proxy of their abundance in April and May 46 .
During refuelling in the Wadden Sea, most godwits were captured, colourmarked, and resighted near Terschelling island (53.40°N, 5.34°E), 50 km from the lugworm sampling area at Balgzand 8 . To justify the use of the Balgzand data for statistically explaining refuelling rates of the Terschelling-captured godwits, we compared lugworm densities between these two areas using data from 2008 to 2014 collected by the Synoptic Intertidal Benthic Survey (SIBES) program, a large-scale grid sampling of benthic fauna in the Wadden Sea 47,48 . We selected stations within 500 m of the shoreline of Terschelling, as most godwits fuel within this range 8 , and compared mean lugworms densities for each year between areas. Lugworm densities were highly correlated (Pearson's r = 0.86) between the two areas.
Data collection on refuelling godwits in the Wadden Sea. To estimate population-level refuelling rates of godwits in the Wadden Sea, we used sex and body mass records from the birds captured at the two main sites: Castricum and the Dutch Wadden Sea. Godwits arriving from wintering areas in West Africa migrate along the Dutch coast before landing in the Wadden Sea. With song playbacks and decoys, overflying birds were lured into landing at a site near Castricum and captured with modified double clapnets for finches powered by elastic cords. There are no foraging areas for migrating godwits near the Castricum catching site thus birds caught there provided samples representative of migrants arriving after a long non-stop flight 49 . Between 1992 and 2016, members of the Castricum Ringing Group captured, ringed and measured 2722 adult godwits 50 . These data were used to estimate annual sex-specific body mass at arrival. Between 1992 and 2016, 6251 refuelling adult godwits were captured across the Dutch Wadden Sea (4.74-6.21°E), mostly around high tide, with wind driven and pulled wilsternets 9 .
Arrival mass and refuelling rates and time. The details on how we estimated annual refuelling time and refuelling rate are outlined in Fig. 3. With the measurements obtained from arriving godwits captured at the Castricum ringing station we estimated mean arrival mass of godwits W0 ks for each year k and sex s and sex-specific standard deviation of residuals σW0 s using lme4 package 51 .
Refuelling time RT k ð Þ was obtained as the difference between year-specific times of departure TD k ð Þfrom and arrival T0 k ð Þto the Wadden Sea, Fig. 3a. TD k was calculated as observed arrival dates at the breeding grounds on the Taimyr Peninsula minus 5.5 days (estimated duration of migration from satellite telemetry data; see Results) and T0 k is the arrival date to the Wadden Sea, estimated with the arrival model (eq. 6).
Godwits are sexually dimorphic, with males being smaller than females, and thus population-level refuelling rates α ks were modelled for each year k and sex s using the yearly arrival dates (T0 k ), arrival mass at Castricum W0 ks (Fig. 3b and d) and the mass of refuelling godwits captured in the Wadden Sea ( Fig. 3c and e). We also assessed linear effects of abundance of lugworms on refuelling rates: α ks 2 Norm μα ks ; σα 2 where μα ks is the expected and α ks are realized refuelling rates. Individual mass at capture in the Wadden Sea W i at Julian day of capture (T i ) was assumed to be influenced by arrival mass W0 ks , the product of refuelling rate α ks , and the time since arrival T i À T0 k : where W0 iks 2 Norm W0 ks ; σW0 2 Original data and model predictions are presented in Supplementary Fig. 1  and 2.
Effects of refuelling time and rate on survival probability. To estimate effects of departure fuel load on subsequent survival, we individually colour-marked and resighted godwits in the Wadden Sea from 2001 9,52 . We used data from adult birds captured in April-May in The Netherlands in 2002-2016 (3995 individuals). Of these birds, 2252 individuals were resighted at least once, yielding a total of 4813 resightings. It has been suggested that departure fuel load will affect subsequent survival (and even reproduction) in migratory birds 53 . During the spring refuelling period godwits are time limited 7 , hence the amount of fuel stored for migration is a product of the refuelling time and refuelling rate 13 . We estimated the effect of fuel load on subsequent survival in the migrants by regressing annual survival Φ ik over year-specific mean refuelling time (RT k ) and mean refuelling rate (μα ks , eq. 3).
To model multiplicative effect of these two variables in the additive framework we used their natural logarithms.
We used the Cormack-Jolly-Seber (CJS) model to estimate apparent survival probability independently from the resighting probability 54 . Survival probability Φ k was modelled as a function of time since marking (TSM, with two classes-'first year after marking' and 'later years'), sex, refuelling rate and refuelling time and their two way interactions. We used mean estimated refuelling rates and durations without accounting for their uncertainty. Because refuelling rate and time are multiplicative we used their natural logarithms to model them in an additive framework: Resighting probability P was modelled as a function of year The set of nested models was estimated in program MARK 55 , using the RMark 56 interface. The model results were corrected for overdispersion of 1.118 ± 0.005, estimated by the median c-hat test in program MARK. Apparent survival probabilities from the CJS model differ from true survival probabilities as they include permanent emigration, but in the case of the godwits, absence of refuelling sites other than the Wadden Sea and thus absence of emigration opportunities makes these two probabilities the same.
Godwit population trends from winter counts in West Africa. Every November or December, 2002-2016, we counted godwits at six high-tide roosts near the village of Iwik in Parc National du Banc d'Arguin, Mauritania (19.8°N, 16.3°W). These counts represent only a portion of the wintering population, and exchanges of individuals may occur with other wintering sites along the African coastline 52,57 . The wintering population counts had high variance not permitting integration with annual survival probabilities, so we used them only to model time-independent population growth λ Thus, we modelled the change in local numbers in a statespace model similar to the Kéry & Schaub approach 43 . The approach assumes that there were only random, but no systematic, shifts in the distribution over the observation period (for which there is no evidence 57 ).
As in equation 2, we assumed a random observation process error following a negative binomial distribution but with a site-dependent error: A negative binomial error distribution was chosen to allow for overdispersion, as birds were counted at high tide roosts in large flocks 41 .
Ethical statement. This research complied with the ethical guidelines of the Dutch law on animal experiments, and was supervised by the Dutch Central Committee for Animal Experimentation (CCD) under guidance of protocol AVD8020020171505 to NIOZ.
Code availability. All data used, R code to perform statistical analyses and model results 31  Arrival date to Wadden Sea ± s.d.
Refuelling rate females (g d -1 ) Refuelling rate males (g d -1 ) Arrival mass ± s.d., females (g) Arrival mass ± s.d., males (g) 1998 Fig. 3 Estimation of refuelling time and refuelling rate of godwits in the Wadden Sea for a sample year. a Annual refuelling time for year k (k = 1998 in the figure) RT k in the Wadden Sea was estimated as a difference between average arrival and departure dates. Mean arrival date to the Wadden Sea T0 k , and its standard deviation σT0 k were estimated from citizen science data on arrival date accounting for observation duration and for variation in observation efficiency between observation sites. Dates of departure from the Wadden Sea were obtained by subtracting the estimated time taken by the migration between Wadden Sea and Taimyr (5.5 days) from dates of first arrival at Taimyr. b Annual arrival mass for females W0 k female was estimated from godwits captured immediately upon arrival from West Africa birds in Castricum. c Population-level female annual refuelling rate α k female estimation combined arrival date T0 k and arrival mass W0 k female estimates with body mass values obtained from godwits refuelling in the Wadden Sea. d For males, arrival mass W0 k male and e refuelling rate α k male were estimated separately as males fuel up slower but are lighter and need to accumulate less fuel for migration