Transmitter attachment
To monitor godwit movements during northward migration, we captured godwits during the non-breeding season on the Chiloé Archipelago, Chile (41°49′S–73°37′W and 42°28′S-73°41′W) prior to migratory departure (see details in [26]) and affixed solar-powered satellite transmitters to 54 adults. We chose both males and females of larger body sizes to control for the effects of transmitter weight; in all cases, the combined weight of transmitters and harnesses comprised < 3% of body mass at capture. In total, we attached 29 Argos Solar Platform Terminal Transmitters (PTTs) weighing 5 g (Microwave Technologies Ltd.) and 25 Pinpoint GPS/Argos Solar S Transmitters weighing 6.6 g (Lotek) using leg-loop harnesses made of nylon or silicon cord.
Both the frequency and precision of location estimates varied between transmitter models. The 5-g devices transmitted Argos locations according to a pre-programmed duty cycle (5-h transmitting/24-h charging); the 6.6-g devices transmitted GPS locations every 2 h and Argos locations opportunistically. To reduce spatial error, we retained standard-quality Argos locations (classes 3, 2, and 1) and passed those of lower-quality (classes 0, A, and B) through a Douglas-Kalman hybrid filter with a maximum realistic rate of movement of 130 km h−1 [27] and a maximum redundant distance of 5 km [28]. For GPS locations with a ± 10-m resolution, we applied a rate-based filtering algorithm that removed locations requiring similarly implausible flight speeds (> 130 km h−1).
Preparing flight tracks
Godwits from Chiloé mainly breed in Alaska. To reach their breeding sites, godwits make a series of flights beginning with a transhemispheric non-stop flight that proceeds over the Pacific Ocean, Central America, and the Gulf of Mexico along a south/north axis. After crossing the Gulf of Mexico, the remainder of their migratory journey occurs over land (Fig. 1; [29]). We focused on godwit behavior during the initial flight stage—which spans from Chiloé to the northern Gulf Coast (hereafter, ‘NGC’) and is generally flown non-stop [25]—for several reasons. First, flights of this nature are seldom featured in behavioral studies (but see [30]), particularly at this level of precision. Second, this flight also improves interpretability of our analyses, as both vertical and lateral movements here are believed to be influenced primarily by wind conditions aloft. By contrast, overland flights after godwits cross into North America are more complex, often involving lengthy searches for spatially dispersed wetland stopover habitat. Behavior over land thus may be influenced by a combination of wind conditions and underlying habitat suitability, which we do not consider here.
To establish godwits’ intended migratory direction, we used their convergence on a main stopover region in midcontinental North America. The region is relatively narrow—encompassing the eastern portions of Kansas, Nebraska, South Dakota, and North Dakota—and contains a high concentration of the ephemeral seasonal wetlands that godwits use to refuel. While godwits do not revisit particular stopover sites, their convergence on and passage through this region is nonetheless highly predictable [25]. We adapted methods from Pearse et al. [31] to identify where godwit tracks exhibited the least latitudinal dispersion in the approach to and while crossing this main stopover region (see Additional file 1 for details). The western and eastern boundaries of this area, as well as their geographic midpoint, served as approximations of the godwits’ overall preferred direction. Dispersion distances and midpoints were calculated on an ellipsoid using functions from the geosphere package [32].
We selected only movements within godwit flight tracks that were associated with their initial flight. From the first offshore location indicative of directional movement away from Chiloé, we used current ground speed and straight-line flight distance from the last onshore location to estimate an hourly time window during which an individual likely departed. All locations prior to this time were removed. Next, tracks were truncated at the first overland location after crossing the Gulf of Mexico, where many individuals began searching for stopover sites. In between these endpoints, we calculated distance, ground speed, turning angle, and heading from each location to the next consecutive location using functions from the R package move [33]. Overland locations where ground speeds fell below 3 ms−1 [34] and consecutive movements traversed < 15 km were indicative of either stopover or transmitter failure and removed from flight tracks. The result was a northward track for each individual, comprising either Argos or Argos and GPS locations, and reflecting directional migratory flight from Chiloé to the NGC [35].
Finally, we performed two additional data thinning steps in preparation for behavioral analyses. To reduce the influence of error and spatial dependence in closely spaced locations, we thinned tracks to a maximum of one location per hour (i.e., the temporal resolution of our wind data), removing locations with lower quality or, where quality was equivalent, locations at random. To reduce the likelihood of dramatic shifts in wind conditions between very widely separated locations, we also removed locations preceding a reporting gap > 12 h. This excluded locations immediately preceding the regular 24-h “off” period in 5-g transmitters and locations preceding consecutive incidental gaps in 6.6-g transmitters, which became especially pronounced in one transmitter deployed continuously for several years. The resulting tracks (n = 28, with n = 689 locations total) contained, on average, one location every 3.25 ± 0.73 h for the 5-g devices and one every 3.15 ± 0.41 h for the 6.6-g devices (Fig. 1).
Wind conditions during flight
We linked each location with the range of possible wind conditions that godwits were likely experiencing at the time. Conditions were retrieved from the European Center for Medium-Range Weather Forecasts ERA5 climate reanalysis dataset, retrieved from the Copernicus Climate Data Store [36]. For each godwit location, we extracted the east/west (u) and north/south (v) components of the wind vector at the nearest hour via bilinear interpolation of the nearest points, and for seven altitudinal pressure levels within the known range of avian flight heights: 1000, 925, 850, 775, 700, 600, and 500 hPa, corresponding to altitudes of roughly 100, 750, 1500, 2250, 3000, 3400 and 5500 m above sea level (m.a.s.l.), respectively. We then calculated the total magnitude (ms−1) and mathematical direction of the wind flow (i.e., the direction the wind is flowing towards) for each location at each altitude using vector trigonometry.
The exact wind conditions that godwits experienced at a location depended on their flight altitude. Though we lacked direct information on altitude, migrating birds are known to alter their flight height dynamically in response to atmospheric conditions [23, 37]. In particular, Black-tailed godwits (Limosa limosa) migrating from Europe to Africa change altitude in order to exploit optimal wind conditions and reduce the thermoregulatory burden of powered flight over the Sahara [27]. Since our godwits likely did not experience similarly extreme temperatures during this stage of migration—they travel during the austral fall/boreal spring along a largely transoceanic corridor—we considered wind support alone to be the best predictor of altitude during this period of continuous flight.
To validate this assumption and define the range of possible flight altitudes, we analyzed observed ground speeds in relation to wind conditions aloft. Ground speed remains among the best indicators of the winds birds are experiencing during flight and is commonly used to estimate flight altitude [38]. In total, we built seven linear mixed-effects models (LMMs) to evaluate how ground speed is affected by wind conditions at various altitudes. Three models allowed individuals to fly at a constant altitude (either 100, 750, or 1500 m.a.s.l.). The remaining four models allowed individuals to fly at the altitude offering optimal support, calculated either in the preferred direction of flight (θd) or the direction of the next location along the realized flight track (θt), selected from either the lower five altitudes or all seven altitudes. Wind support (ws, ms−1) was calculated as:
$$ws = \cos \left( {atan2(u,v) - \frac{\theta \pi }{{180}}} \right) \times \sqrt {u^{2} \times v^{2} }$$
In this equation, \(\theta\) refers to the direction of travel, while u and v refer to the east/west and north/south components of the wind vector, respectively. These models, as well as all LMMs hereafter, initially included year and individual as random intercepts to account for spatial dependence and individual differences [39]. Models were constructed with the R package lme4 [40] and compared using corrected Aikake’s Information Criterion (AICc) scores [41]. We selected the model with the lowest score or, where competing models demonstrated no substantial difference (ΔAICc < 2), the simpler of the models.
With this information, we assigned the most likely wind conditions to each godwit location. To evaluate changes in these conditions, we organized the migratory corridor into seven latitudinal bands (~ 10° each) and used circular statistics to evaluate within-band scalar means of wind speed, circular mean directions (\(\overline{\theta }\)), and the mean resultant length (\(\overline{R}\)) of flow vectors, where \(\overline{R}\) is a statistic between 0 and 1 that indicates the spread of the mean direction (i.e., 0 = large spread; 1 = concentrated at a single value). We confirmed the non-uniformity of wind flow within each band using Rayleigh tests. To assess the robustness of our assigned wind conditions to imprecision in altitude estimation, we used a Pearson’s product-moment correlation for wind direction and wind support at two different altitudes: the altitude offering the maximum wind support and the altitude offering the next-highest wind support. Analyses of circular data were performed with functions from the circular package [42].
Classifying in-flight behavior
To test our behavioral predictions, we assigned one behavior to each godwit location based on the direction of wind flow and the godwit’s subsequent movement. While several methodologies have already been developed for behavioral classifications (e.g., [13, 43]), they generally require a single location that can serve as the migratory destination, such as a fixed breeding or stopover site toward which an individual consistently navigates. In our case, and indeed for many other species that exhibit more flexible movements during the early stages of migration [44], these methodologies may be problematic. Godwits navigate toward a range of suitable stopover destinations in midcontinental North American [25]. To preserve this range in our calculations, we developed a new workflow for assigning in-flight behaviors, which utilizes an intuitive circular framework to compare the directions of movement vectors.
We first confirmed unimodality in the travel headings of our tracked godwits, both by visual inspection and by comparing results from maximum likelihood models of circular orientation developed for animal movement via the CircMLE package [45]. Next, we defined two normalized (i.e., direction-only) movement vectors for each location: the total wind flow (\(\hat{w}\)), and the realized direction of godwit travel to the next location (i.e., the track direction, or where a godwit actually flew; \(\hat{r}\)). To allow for orientation to a geographic region rather than a particular site, and to account for the fact that the proportion of the flight horizon that encompasses this region will increase with greater proximity, we defined a population-wide range of preferred directions. This range was bounded by the western- and easternmost points in the narrowest region of the migratory corridor (\(\left[ {\widehat{{d_{w} }},\widehat{{d_{e} }}} \right]\)) and had as its center the angular midpoint of this range (\(\hat{d}\)).
We then calculated the angles separating these normalized vectors (i.e., their “closeness”), focusing on the angles between the realized travel direction and the wind (θrw), the realized travel direction and the midpoint of the preferred range (θrd), the midpoint of the preferred range and the wind (θwd), and the midpoint and edge of the preferred range (θdd), where θ ∈ [θ, π]. To account for imprecision in these measurements, we incorporated two tolerance values: 1) a larger tolerance value (τ1 = 0.2 rad or ~ 11.46°) to account for the combined imprecision of the realized and wind directions, and 2) a smaller tolerance value (τ2 = 0.1 rad or ~ 5.73°) to account for the imprecision of the realized direction alone, which is generally small for infrequently sampled, fast-flying species ([34]; see Additional file 1 for details on tolerance value estimation). Tolerance values in this way allowed us to account for known and variable imprecision in our data without adding additional unknown sources of imprecision. Together, this information enabled us to differentiate between fully drifting flight (realized direction ≈ direction of the wind flow, but not within the preferred range), supported flight (realized direction ≈ direction of the wind flow and within the preferred range), and fully compensating flight (realized direction ≠ direction of the wind flow and within the preferred range).
If none of these categories applied, we inspected the arrangement of vectors (i.e., their “betweenness”). The cross products \(\left\| {\widehat{r} \times \widehat{d}} \right\|\) and \(\left\| {\widehat{w} \times \widehat{d}} \right\|\) returned orthogonal vectors with the sign determined by the degree of rotation from one angle to another. Comparing the signs of cross products thus revealed which vector lay between the other vectors on a semi-circle. Different signs indicated overcompensating flight (midpoint of the preferred range lies between the realized direction and wind flow). Similar signs indicated either partially compensating flight (realized direction lies between the wind flow and the midpoint of the preferred range) or overdrifting flight (defined here as when the wind flow lies between the realized direction and the midpoint of the preferred range). Together, these steps produced an easily implemented, hierarchical decision list that assigned to each in-flight location one prevailing behavior (Fig. 2).
To test our predictions that drift would prevail in the early stages of the flight from Chiloé to the NGC, and particularly over open water, we first examined conditional associations between assigned behaviors and latitudinal bands. We applied an omnibus χ2 test for independence and used a post-hoc inspection of the adjusted standardized residuals with a Bonferroni adjustment (⍺ = 0.0014, critical value = ± 2.98; [46] to evaluate cell contributions to significance, using functions from the stats [47] and vcd [48] packages. Because overdrift was exceptionally rare, we verified that it was distributed across several latitudinal bands and excluded it from the χ2 test.
Since categorizing behavior and binning by latitudinal bands may obscure more subtle patterns, we also performed a companion analysis investigating ongoing changes in godwit divergence from the preferred range (θrd) with changes in latitude, wind support, and crosswind using multivariate adaptive regression spline (MARS) models implemented via the earth package [49]. MARS models use partitioning and iterative pruning to capture possible interactions and nonlinearities in predictor variables. Because MARS models do not have the ability to deal directly with outliers, we inspected for and removed any extreme outliers (θrd > 1.5 rad or 90°, n = 1) prior to analysis. We centered and standardized our predictor variables prior to inclusion and set minspan = 20 to limit over-fitting in sections of the migratory route with few observations. We tested for the possibility of interactive or additive effects and allowed variables to enter linearly or with nonlinear change-points (i.e., hinge functions). Results obtained with MARS and similar machine-learning methods are generally robust to collinearity in the set of predictors [50]; nonetheless, we verified that the variance inflation factors of variables were suitably low (VIF < 2). We report the unscaled, uncentered MARS model parameter estimates.
Finally, to compare our results with patterns observed in another individual tracking study, we adapted methods from Klaassen et al. [13] to evaluate the effect of the forward component of the wind vector in the preferred direction (wind support) on the forward rate of godwit movement in the preferred direction (rv), as well as the effect of the lateral component of the wind vector (crosswind) with the lateral rate of godwit movement (ru). Individual and year were again included as random effects. We also adapted their calculation for total drift, which is the ratio between the slopes of these models (βcw/βws). Residuals from these models informed where godwit lateral movements were most strongly influenced by crosswind.
Crossing over the Northern Gulf Coast
We predicted that the cumulative effects of wind flow during flight would influence where godwits crossed into North America. To test this, we identified the location of NGC crossover (i.e., crossing over land, but not necessarily stopping) for each complete flight track. When a track contained an onshore location estimate < 10 km from the coast, we used this location. For other tracks, we interpolated a plausible crossover location along the straight-line intersection of the flight track with a coastline shapefile (Natural Earth, 10-m global coastline). For each individual, we calculated its total longitudinal displacement by subtracting its actual longitude at the point of NGC crossover from the crossover longitude of a shortest-distance (i.e., Great Circle) line extending from Caulín on Chiloé to the geographic midpoint of the preferred range. We used a linear regression to assess the extent to which this displacement was shaped by mean crosswinds aggregated over all preceding locations along the flight track. To control for displacement associated with habitat selection during unexpected stopover events in Central America, individuals with complete tracks that stopped prior to the NGC (n = 3) were removed from this analysis.
We present results as means ± standard deviations or medians and interquartile ranges, unless otherwise stated. The global significance of the fixed effects in all models was approximated with parametric bootstrapped p values (1000 simulations) obtained with the package afex [51]. We evaluated the uncertainty of parameter estimates using 95% parametric bootstrapped confidence intervals (1000 simulations) obtained with the ‘confint.merMod’ function in lme4 [40]. Estimations of the partial variance explained by fixed effects (R2LMM(f)) were calculated using ‘rsq.lmm’ function via the rsq package [52]. Random intercepts for year and individual were retained if their variance contributions were greater than (near) 0 and fits did not result in singularity or non-convergence [53]. All analyses were conducted in the R programming environment, version 4.0.3 [47].