Warming sea surface temperatures are linked to lower shorebird migratory fuel loads

Warming sea surface temperatures (SSTs) are altering the biological structure of intertidal wetlands at a global scale, with potentially serious physiological and demographic consequences for migratory shorebird populations that depend on intertidal sites. The effects of mediating factors, such as age-related foraging skill, in shaping the consequences of warming SSTs on shorebird populations, however, remain largely unknown. Using morphological measurements of Dunlin fuelling for a >3000 km transoceanic migration, we assessed the influence of climatic conditions and age on individuals’ migratory fuel loads and performance. We found that juveniles were often at risk of exhausting their fuel loads en route to primary wintering grounds, especially following high June SSTs in the previous year; the lagged nature of which suggests SSTs acted on juvenile loads by altering the availability of critical prey. Up to 45% fewer juveniles may have reached wintering grounds via a non-stop flight under recent high SSTs compared to the long-term trend. Adults, by contrast, were highly capable of reaching wintering grounds in non-stop flight across years. Our findings suggest that juveniles were disproportionately impacted by apparent SST-related declines in critical prey, and illustrate a general mechanism by which climate change may shape migratory shorebird populations worldwide.


Introduction
Warming sea surface temperatures (SSTs) and increasingly variable precipitation patterns are altering species interactions at global scales.These new interactions are impacting physiological conditions and processes of individuals that scale up to shift population demographic trends [1,2].In seasonal environments,

Study population
The pacifica subspecies of Dunlin is a population of small (50−60 g; [32]) migratory sandpiper endemic to the Pacific Americas Flyway [30,33].Individuals arrive on YKD breeding grounds in early May, coincident with the availability of snow-free habitat [34], and nest on the ground in low-lying wet graminoid meadows from late May to late June [35].Eggs typically hatch after 21 days of incubation, and chicks fledge at 21−26 days old [36].In July, fledged juveniles and post-breeding adults relocate to the YKD coast where they feed extensively on Limecola balthica clams available in the region's vast intertidal mudflats [37], and moult their body feathers (juveniles and adults) and the previous year's flight feathers (adults only; [38]).By late August, more than 100 000 juvenile and adult Dunlin have relocated to the YKD coast [22], completed or nearly completed their moult [38], and begun fuelling for a >3000 km migration across the Gulf of Alaska to wintering grounds on the Pacific coast of North America (figure 1; [29][30][31]).During the fuelling period, a portion of the arcticola subspecies of Dunlin that breeds in northern Alaska and winters in East Asia also stages in the YKD [30,39].Juveniles and adults stage into late September and early October and depart synchronously [22,31], with heavier individuals tending to depart 2.5 days earlier per additional gram of body mass at time of capture [31].

Statistical analyses 2.2.1. Pre-migratory body condition
Dunlin body mass varies among and within individuals at different points in the year.Part of the variation in mass is owing to: (i) differences in individuals' structural size [32], and part to (ii) individuals seasonally gaining and losing fat and fat-free mass (i.e.lean body mass) in response to life history demands (e.g.migration; [37]).To assess annual variation in individuals' migratory fuel loads-i.e. the amount of residual lean mass and fat individuals were carrying at the time they were captured and weighed-we therefore had to first control for differences in individuals' structural size, or structural lean body mass-i.e. the mass of individuals' structural body components (e.g.skeleton and respiratory system)-and control for differences in individuals' capture dates.
To estimate and control for differences in individuals' structural lean body masses, we first fitted linear mixed models with mass at capture as the response variable; culmen length, culmen length + age or an interaction between culmen length and age as explanatory variables; and year as a random intercept using the R package lme4 (v.1.1−34; [41]).We restricted this analysis to captures that occurred during 1−21 August because this period precedes completion of the pre-basic moult and is before fuelling for southbound migration begins, thereby allowing us to assume that captured individuals had minimal amounts of residual lean mass and fat [38,42].We chose to regress mass by culmen length because other morphological measurements were not always collected, and because previous studies show that culmen length alone is highly correlated with structural size across Dunlin subspecies and sexes [42,43].We included age as an explanatory variable to account for potential differences between recently born juveniles and adults in culmen length-structural lean body mass relationship.Finally, we restricted this analysis and all other analyses to captures of individuals with a culmen length >31.5 mm, thus removing males of the arcticola subspecies (<2% of captures) that have a structural lean body mass-culmen length relationship different than that of individuals of the pacifica subspecies [32,33].Individuals of the arcticola subspecies with a culmen length >31.5 mm could not be separated from pacifica individuals and were, therefore, included in our analyses.Prior analyses indicated that there was not a difference in fuel loads of individuals of unknown subspecies (i.e.individuals with a culmen length >31.5 mm and <40.9 mm) and those of female pacifica with a culmen length >40.8 mm, indicating that our results are not biased by our inability to remove certain individuals of the arcticola subspecies (electronic supplementary material, figure S1).
royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 240324 Using Akaike information criterion (AIC), we found that the linear model containing an interaction between culmen length and age was the most parsimonious (electronic supplementary material, table S1, figure S2).We then used the parameter estimates from this linear model to estimate structural lean body masses and, therefrom, fuel loads of individuals captured during the migratory fuelling period (22 August−22 September).Specifically, we estimated an individual's migratory fuel load by subtracting their structural lean body mass from their body mass upon being captured.This approach assumed that the culmen length-structural lean body mass relationship identified in pre-fuelling individuals (captured 1−21 August; electronic supplementary material, figure S2) was identical to that of fuelling individuals (captured 22 August−22 September).We standardized estimated fuel loads by dividing an individual's fuel load by their structural lean body mass.Fuel load estimates, therefore, represent an approximation of the relative fuel load individuals were carrying at the time they were captured and weighed (as a % of their structural lean body mass [%lbm]).
To control for differences in individuals' capture dates, we next calculated mean daily fuel loads during the migratory fuelling period by fitting gamma-distributed generalized linear mixed models [41] with a log link for juvenile and adult Dunlin using migratory fuel load at capture as the response variable, date as an explanatory variable and year as a random intercept.To ensure juveniles and adults had equal representation across years and sampling periods in these baseline models, we modified the datasets used to fit our models as follows: (i) we divided the migratory fuelling period into four 8 day sampling periods; (ii) we identified the minimum number of age-specific captures per sampling period-year combination; and (iii) we subsampled (without replacement) the age class with  S2).We then estimated individuals' pre-migratory body conditions-i.e. the extent to which individuals' fuel loads were above or below mean fuel loads, given individuals' capture dates-by subtracting the mean fuel load on an individual's capture date (given our baseline models) from their estimated fuel load.Positive values indicated an above-average fuel load after accounting for an individual's age, structural size and capture date, and negative values indicated a below-average fuel load [44].

Climatic effects
To explore climatic effects on annual variation in individuals' migratory fuel loads, we conducted sliding window analyses using the R package climwin (v.1.2.3; [45,46]).In short, a sliding window analysis systematically explores the influence of a weather variable on a biological response across a range of defined time periods.We used a linear mixed model with pre-migratory body condition as the response variable and year as a random intercept as our baseline model.We then identified the time period candidate weather variables (average daily air temperature [°C; avg air temp], average daily SST [°C; avg SST], avg SST in the previous year [year i-1 ], and average daily precipitation [mm; avg precip]) explained the most variation in Dunlin pre-migratory body condition by iteratively adding to the baseline model the mean of a given weather variable within a given candidate climate window, and ranking models using AIC corrected (AICc) for small sample sizes.Among our candidate weather variables, we assessed SSTs from the year before birds were captured, in addition to SSTs in the year of capture, to assess the hypothesis that SSTs in the previous year could affect the availability of critical prey (i.e. 1 year old L. balthica clams) and, thereby, affect Dunlin pre-migratory body conditions.This hypothesis is based on the following observations from the YKD and the Wadden Sea region of the North Sea: (i) in Angyoyaravak Bay, YKD, Dunlin fuelling for southbound migration feed extensively on L. balthica clams (e.g.97% of prey items in 38 stomachs examined; [37]); (ii) annual recruitment of 0 year old L. balthica clams varies widely between years and is associated with temperature-driven variation in predator abundances and timing of phytoplankton blooms [47,48]; (iii) recruitment of 0 year old L. balthica clams drives the subsequent availability of 1 year old clams [49,50]; and (iv) small (1 year old) L. balthica clams have greater energetic value and are preferentially eaten at higher rates by a closely related shorebird species than larger (older) L. balthica clams [51,52].We allowed the start and end dates of candidate climate windows to vary (or 'slide') at weekly intervals across pre-breeding, breeding and molting periods (1 May−21 August; candidate window range: 1−17 weeks in duration), thereby capturing the environmental and biological processes that likely affect Dunlin fuel loads, while also assessing identical climatic conditions across our sample of individuals and reducing the effects of anomalous weather events and spurious correlations that may occur over short time periods [46,53].We fitted candidate models with linear and quadratic response functions, and considered models with ΔAICc < 2 to have equal support.Comparing a large set of candidate models may lead to spurious results and uncertainty in selecting a best model because many models may have similar levels of support [45,46].To assess the likelihood that top models represented spurious false-positive results, we used the climwin randomization function to compare observed levels of model support with levels of model support from 10 randomized versions of our dataset.We considered models with a P C value <0.10 to represent a real climate signal [46].We then built a global model, and all possible subset models (i.e.'all subsets' models), with pre-migratory body condition as the response variable, influential weather variables and secondorder interactions as explanatory variables, and year as a random intercept.We assessed collinearity between influential weather variables using Pearson's product-moment correlations.We excluded subset models with correlated linear terms (r > 0.7) but retained models with correlated interaction terms.Ultimately, our 'all subsets' approach identified a single model, with just one influential weather variable (SST year i-1 ), that performed substantially better than all other models (i.e.ΔAICc > 2).We then returned to the 'sliding window' model set (see above) for the weather variable in the top all subsets model and averaged candidate models that had an Akaike model weight > 0.05 (i.e. the top all subsets model and corresponding sliding window models with an Akaike model weight > 0.05).Daily air and SSTs and precipitation values for the central YKD (60.5°−62°N, 164.75°−166.75°W; hourly data on single levels; 0.25° × 0.25° gridded spatial resolution) were acquired from The European Center for Medium-Range Weather Forecasts ERA5 reanalysis data set [54].We restricted analyses of climatic effects to juveniles because adult sample sizes were not adequate for conducting sliding window royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 240324 analyses (i.e.there were <9 years with >25 captures per year; electronic supplementary material, table S3; [55]).

Flight range predictions
To understand how annual variation in individuals' migratory fuel loads might affect individuals' migratory performance, we conducted flight range simulations using program Flight [56,57] implemented in the R package FlyingR (v 0.2.2; [58]).In short, program Flight combines a biomechanical flight model with a bird's physical and physiological characteristics to predict its maximum flight range (electronic supplementary material, table S4).One key variable in flight range prediction is departure fuel load-i.e.residual lean mass and fat upon departure.We predicted juvenile and adult departure fuel loads using capture records from early October; the period of the year when individuals are at their heaviest and most likely to depart on southbound migration [22,37].Specifically, we drew 100 departure fuel loads for juveniles and adults from normal distributions with means and standard deviations informed by October capture records (n = 5 and 5, respectively; electronic supplementary material, table S3).Although limited, the number of October capture records we used represent, to the best of our knowledge, the only available information on near-departure body masses of Dunlin from the YKD.Because all October captures occurred in the same year (1980), we then adjusted departure fuel loads according to variation in mean annual pre-migratory body conditions and mean climatic effects identified in sliding window analyses (electronic supplementary material, table S5).This approach assumed: (i) that body conditions during the fuelling period reflected body conditions at departure, and (ii) departure timing was fixed (independent of body condition), both of which are simplifications of complex processes.These simplifications were necessary for assessing the influence of climatic conditions and age on individuals' ability to undergo the population's preferred migratory route.The potential impact of these assumptions on interpretation is addressed in §4.
Another key factor in flight range prediction is wind assistance, as Dunlin likely time their southbound migration to coincide with supportive wind conditions [37].We accounted for wind assistance in our flight range predictions following a 3-step process.First, we calculated the average bearing of the shortest route to the population's primary wintering grounds (i.e.112.5°; electronic supplementary material, figure S3).Second, we calculated average daily wind assistance values (given a bearing of 112.5°) for the entire migratory corridor (spanning ~250 km on either side of the shortest migratory route; electronic supplementary material, figure S3) from 29 September to 11 October 1978−1980, 2004−2006 and 2008−2010 [37].Average daily wind assistance values were estimated at Pre-migratory body conditions (i.e. the extent to which fuel loads were above or below average after correcting for individuals' age, structural size and capture date [% structural lean body mass; %lbm]) were estimated using age-specific linear and generalized linear mixed models and, therefore, are not equivalent between juveniles and adults (see text and figure 1).

Climatic effects
Average SST in the previous year (SST year i-1 11-24 June) and average precipitation in the year of capture (4 June−29 July) were influential predictors of juvenile pre-migratory body condition and were negatively correlated with each other during our study (r = −0.87).For both weather variables, quadratic response functions fit better than linear responses (electronic supplementary material, table S6).Comparison of the top model for each weather variable and a model including their interaction suggested SST in the previous year (SST year i-1 11-24 June) alone was the best predictor of juvenile pre-migratory body condition (electronic supplementary material, table S7).Four SST year i-1 sliding window models had Akaike model weights >0.05 and were averaged (electronic supplementary material, table S8).Intermediate SST year i-1 values were associated with above-average pre-migratory body conditions, and high SST year i-1 values were associated with the lowest average body conditions (figure 3). .The horizontal dotted line represents the minimum distance between Angyoyaravak Bay and the northern extent of the species' primary wintering grounds on the Pacific coast of North America (see figure 1).

Flight range predictions
predictions under long-term (1978−2022) SST year i-1 conditions ranged from 2493 to 5470 km and were adequate for reaching primary wintering grounds on the Pacific coast of North America under typical long-term climatic conditions.Median flight range predictions fell below the minimum migration distance, however, when SSTs were much warmer (>10.0°C)than climatological means.Up to 45% fewer juveniles may have reached primary wintering grounds via a non-stop flight under recent high SSTs compared to the long-term trend (i.e.10.3°C versus 4.8°C; figure 3).

Discussion
Warming SSTs are altering the biological structure of intertidal wetlands [65] with potentially serious physiological and demographic consequences for migratory shorebird populations that depend on intertidal sites [16,18].Using body mass measurements of Dunlin fuelling for a >3000 km transoceanic migration from the YKD, Alaska, to wintering grounds on the Pacific coast of North America, we show that juvenile pre-migratory body conditions were correlated with June SSTs in the previous year, with high June SSTs driving low migratory fuel loads the following year.The clear lagged nature of the observed effect suggests SSTs likely acted on fuel loads by altering the availability of critical prey.We propose that the lagged effect reflects a link between SSTs and the recruitment of L. balthica clams.
Limecola balthica clams are eaten extensively by Dunlin fuelling in the YKD (e.g.97% of prey items in 38 stomachs examined; [37]). 1 year old L. balthica clams have high energetic value and are likely selected by Dunlin at higher rates than larger (older) L. balthica clams [51,52].Studies suggest the availability of 1 year old L. balthica clams is driven by recruitment of 0 year old clams in the previous year [49,50], which varies widely according to temperature-driven variation in predator abundances [47] and timing of L. balthica spawning activity in relation to spring phytoplankton blooms [48].In the western Wadden Sea region of the North Sea, for example high SSTs are associated with earlier spawning activity of L. balthica clams, creating a mismatch between spawning and spring phytoplankton blooms.Reduction in phytoplankton available to clam larvae caused by this phenological mismatch is believed to limit food and delay the age and decrease the size at which larvae metamorphose into bottom-dwelling recruits, increasing the overall mortality of 0 year old clams [48,66].
In the Bering Sea region of the YKD, June SST is correlated with the timing of spring sea ice retreat and phytoplankton blooms [25,67].In cold years, late ice retreat leads to early ice-associated blooms in cold water (i.e.May), whereas warm years lead to early ice retreat and later blooms in warm water (i.e.June [25,68]; see also [69]).High SSTs and late phytoplankton blooms (i.e.June bloom productivity) are likely better aligned with spawning activity of L. balthica clams than the low SSTs early blooms (electronic supplementary material, figure S5).Despite the possibility of a better match between spawning activity and phytoplankton blooms under high SSTs, studies suggest less late-season phytoplankton ultimately enters the benthic food web owing to higher rates of grazing by zooplankton compared to years with colder SSTs and earlier blooms [26,68] (see also [70]).High SSTs in the YKD, therefore, may result in less phytoplankton available to developing larvae and 0 year old L. balthica clams, decreasing larval recruitment, and subsequently limiting Dunlin fuelling for southbound migration the following year by reducing the availability of 1 year old clams [71].Particularly cold (i.e.very cold) SSTs could also reduce the availability of 1 year old clams via later ice retreat, poorer blooming conditions and comparably weaker phytoplankton blooms [26], suggesting the observed quadratic effect of intermediate (i.e.cold) SSTs being associated with above-average fuel loads.
The YKD, and the high-quality prey therein, represents a potential bottleneck for many southbound-migrating shorebird populations in the North Pacific [22,72].Changes in availability of critical prey within the YKD, such as 1 year old L. balthica clams, could alter competitive interactions between shorebird species, sexes and age classes and have far-reaching population-level effects.For example, when critical prey resources initially become scarce, better competitors (i.e.adults with better foraging skills) are usually able to acquire adequate resources, whereas poorer competitors (i.e.juveniles) are often disproportionately impacted [8,9].Through an increase in juveniles failing to recruit to the breeding population, changes in prey availability can have a significant effect on population structure and size [20,73].Among Dunlin in the YKD, we found that juveniles and adults had similar fuel loads at the onset of the migratory fuelling period, but adults attained higher fuelling rates and larger migratory fuel loads.Interannual differences in juvenile fuel loads were in turn explained by June SSTs in the previous year.Taken together, these findings suggest that juveniles were outperformed and/or royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 240324 excluded from preferred feeding areas by adults, and were more likely to suffer negative effects from apparent SST-related declines in critical prey [73].
Consequences of prey limitation can be direct, such as immediate starvation, or indirect such as individuals lacking the fuel loads required to migrate in an optimal way.We generated maximum flight range predictions using body mass measurements that accounted for annual variation in migratory fuel loads and found that adults were highly capable of departing the YKD, crossing the Gulf of Alaska, and reaching primary wintering grounds in non-stop flight, whereas juveniles were often at risk of exhausting their fuel loads en route, especially following high June SSTs in the previous year.Poorer foraging skills and/or exclusion from preferred feeding areas, therefore, may have caused a portion of southbound-migrating juveniles to migrate at the edge of their physical ability.Moreover, relatively small changes in prey availability in the YKD, such as possible reductions in L. balthica clams with ongoing climate change [74,75], may push juveniles to migrate at the edge of their ability at higher rates (e.g.[76]).
The warmest years of our study (2003−2005) coincided with record late sea ice formation and early spring sea ice retreat in the Bering Sea region of the YKD [67].Recent significant marine heat waves in the region (2017−2020; [24]) have since resulted in record low sea ice extents, likely reductions in ice-associated phytoplankton blooms, northward shifts in the distribution of crab and fish populations, unprecedented seabird die-offs and widespread reproductive failures and marine mammal mortality events [77]; all indications that ongoing climate change in the eastern Bering Sea could drive long-term changes in energy flow and ecosystem structure [26,77,78].
We found that recent high SSTs in the Bering Sea region of the YKD (2017−2020) may have resulted in up to 45% fewer juveniles reaching primary wintering grounds via a non-stop flight.Given that Dunlin show an endogenous drive to take a direct migratory route over the Gulf of Alaska [29], and that there is no evidence of similarly favourable alternative routes [37], this finding suggests that warming SSTs in the YKD with ongoing climate change [27] could trigger higher juvenile mortality during southbound migration and a lower stable Dunlin population size [20,79] (but see [10,80]).There is ample evidence, however, that migratory behaviours are flexible [81] and shorebirds may compensate behaviourally for smaller fuel loads.Later arrival of harsh winter conditions in the YKD with ongoing climate change [24], for example, could allow populations to stage later and compensate for lower fuelling rates and early October fuels loads [82].Populations may also stop opportunistically at ancillary stopover sites along the Alaska Peninsula and acquire additional fuel [31,83].
The extent to which juvenile shorebirds may flexibly adjust their endogenous migration programs, however, is largely unknown [84].Studies suggest that juvenile birds, in general, modify their endogenous migration programs as they gain experience [84,85], and have especially high levels of mortality during their first migratory flights [85,86]; levels of mortality that may be moderated by annual variation in climatic conditions and fuel loads [76,87] (see also [13]).At a minimum, smaller fuel loads with warming SSTs may push juvenile Dunlin to take a longer, less efficient and possibly more dangerous coastal migration route rather than a direct over-ocean migratory route [21,88,89], which would represent a significant, climate-induced change in their migratory behaviour.
The extent to which staging shorebird populations are negatively impacted by ongoing climate change may largely depend on the ability of juveniles to flexibly adjust their endogenous migration programs to conditions encountered en route [81,85].Given potentially significant consequences of warming SSTs on staging shorebird populations, additional studies are needed to develop a strong empirical understanding of the links between SSTs, intertidal invertebrate population dynamics and shorebird fuel loads and demographic trends.

Ethics.
Research activities were approved by the U.S. Geological Survey Bird Banding Laboratory (permit no.20022), U.S. Fish and Wildlife Service, Alaska Department of Fish and Game and the U.S. Geological Survey Alaska Science Center's Animal Care and Use Committee.

Figure 1 .
Figure 1.(a) Location where Dunlin were captured (Angyoyaravak Bay) and the minimum migratory distance between the YKD, Alaska, and the northern extent of the species' primary wintering grounds on the Pacific coast of North America (Grays Harbor, Washington).Temporal change in mean fuel loads of (b) juvenile (JUV; n = 590) and (c) adult (ADU; n = 590) Dunlin captured during the southbound migratory fuelling period (22 August−22 September 1978−1980, 2008−2010) in Angyoyaravak Bay.Fuel loads (i.e. the amount of residual lean mass and fat individuals were carrying at the time they were captured and weighed, after correcting for differences in individuals' structural size [% structural lean body mass; %lbm]) were estimated using a linear mixed model (see text).Shaded areas represent 95% confidence intervals (CIs).(d) Mean log-fuel accumulation rate (log[%lbm]/day; ±95% CI) by age.

Figure 2 .
Figure 2. Mean annual pre-migratory body condition of juvenile (JUV; n = 1459) and adult (ADU; n = 1133) Dunlin (±95% CI) captured during the southbound migratory fuelling period (22 August−22 September) in Angyoyaravak Bay, YKD, Alaska.Pre-migratory body conditions (i.e. the extent to which fuel loads were above or below average after correcting for individuals' age, structural size and capture date [% structural lean body mass; %lbm]) were estimated using age-specific linear and generalized linear mixed models and, therefore, are not equivalent between juveniles and adults (see text and figure1).

Figure 3 .
Figure 3. (a) Effect of average SST (11−24 June) in the previous year (year i-1 ) on juvenile Dunlin (n = 1459) pre-migratory body condition (mean ± 95% CI) during the southbound migratory fuelling period (22 August−22 September 1978−1980, 2004−2006, 2008−2010) in Angyoyaravak Bay, YKD, Alaska.Pre-migratory body conditions (i.e. the extent to which fuel loads were above or below average after correcting for individuals' age, structural size and capture date [% structural lean body mass; %lbm]) were estimated using linear and generalized linear mixed models (see text).(b) Average SST in the central YKD by year.The horizontal line represents the median annual SST value from 1977 to 2022.(c) Mean effect of average SST year i-1 on median predicted flight ranges of juvenile Dunlin departing Angyoyaravak Bay on southbound migration (1978−2022).The shaded area represents the interquartile range of predicted flight ranges.Flight range predictions were generated using program FLIGHT and account for wind assistance (see text).

Figure 4 .
Figure 4. Median predicted flight ranges of juvenile (JUV) and adult (ADU) Dunlin departing Angyoyaravak Bay, YKD, Alaska, on southbound migration.Boxes and whiskers represent the inner 50% and 90% of predicted flight ranges.Flight range predictions were generated using program FLIGHT and account for wind assistance (see text).The horizontal dotted line represents the minimum distance between Angyoyaravak Bay and the northern extent of the species' primary wintering grounds on the Pacific coast of North America (see figure1).