Study site and population
The Manawatu River estuary (40.47°S, 175.22°E; Fig. 1) is a small (ca. 1 × 2 km) intertidal mudflat area on the west coast of the North Island, New Zealand. Since 2006, we have captured bar-tailed godwits by cannon-net or mist-net, and marked individuals with a numbered metal band, plus either a unique combination of the white flag and four colorbands, or an engraved white flag with a field-readable three-digit alphabetical code. During 2008–2020, 35% (range: 23–48% per year) of adult (i.e., migratory) individuals were marked in the estuary’s highly site-faithful population of 170–280 godwits.
At capture, godwits were aged (≥/<3 years) by plumage or state of primary feather molt; because godwits ≥3 years cannot be precisely aged, the great majority of birds in our studies are of unknown age. Due to sexual dimorphism (females > males), most (89%) individuals were sexed by bill length (length of exposed culmen: >99 mm = female; <88 mm = male), but intermediate birds (88–99 mm) cannot be sexed by this measure25; however, conspicuous dimorphism in plumage before departure from New Zealand allowed unambiguous sexing of the remaining individuals. In this population, size and migration timing in both sexes vary along a latitudinal cline in the Alaska breeding range (60–71°N; higher latitude = smaller and later migration)24,25.
Directly observed migratory departures
The local godwit population is highly approachable and can usually be observed entirely from one of several vantage points along the perimeter of the mudflat. Each year, a single observer (occasionally two) monitored the flock daily from ca. 1 March until all migratory individuals had departed (29 Mar–5 Apr). During surveys, the observer watched and listened for migratory behavior, which included distinct vocalizations and low circling flights expressing intent to depart. Flocks typically engaged in conspicuous calling, preening, bathing, and short exploratory flights for 0.5–4 h before actual departure15; during this time, the observer recorded all marked individuals involved, using a 20–60× spotting scope and digital camera with a 400 mm lens. All flying flocks were watched and/or photographed until they resettled or disappeared from sight. After each departure, the observer quickly surveyed the estuary to account for all remaining marked godwits. In addition, we conducted daily high-tide surveys to confirm the size and composition of the remaining flock; the average daily resighting probability of marked godwits was 88% (range 51–100% per individual) before departure.
We recorded departures of 44–81 marked individuals per year (174 total individuals). In 86% of these departures, the bird was directly observed preparing and/or departing with an observed flock, and therefore exact time and flock size were known. The remaining individuals were assigned a departure date based on the last day they were recorded at the estuary. For 10% of cases, this coincided with an observed departing flock of partially or completely unknown individual composition; therefore, we considered a time and flock size to be known. For the remaining 4% of cases, the individual’s disappearance coincided with a decrease in local flock size unexplained by observed departures; we considered these departures unobserved, and calculated flock size based on successive high-tide counts. 98% of observed departures occurred during 13:00–21:00, and unobserved departures likely occurred at night or early morning. By including unobserved departures, annual totals represent a virtually complete accounting of flocks and individuals migrating from the site.
We are confident that observed movements out of the estuary represented a migratory departure from New Zealand. Site fidelity of marked godwits is extremely high for the entire non-breeding season (September–March), and nonmigratory movements in and out of the estuary were rare. Departures were easily distinguished from local movements by both altitude and direction; departing flocks always flew NW/NNW and slowly ascended toward the ocean before disappearing from view. Furthermore, for all departures captured by both direct observation and geolocator tracking (n = 51; see below), conductivity (wetness) data recorded by geolocators confirmed a ca. 7-day dry period (indicating a non-stop flight to Asia) starting on the same day (0–2.5 h difference) as the observed departure64.
All statistics were calculated in the software R 3.5.165. Using departures of all godwits (n = 128–251 total individuals in 12–17 flocks on 10–14 days per year), we analyzed phenology change across years using a linear model in the R package lme4 v.1.1.2166.
Including only individuals observed departing New Zealand in ≥3 years (n = 124 individuals; 44–74 per year), we distinguished the contributions of within- and between-individual variation using within-subject centering31 in a linear mixed model in lme4. To explore within-population variation in these effects, we performed the same analysis separately for the two regions in the Alaska breeding range (North, South; Fig. 3), including 34 individuals with known breeding locations from geolocator tracking (see below).
To quantify the proportion of individuals (≥3 observations) that showed an advancement in their departure timing, we obtained individual posterior distributions via direct simulations of 1000 values from the joint posterior distribution of ordinary linear model parameters using the function sim from the R package arm v.1.11-267. The 2.5% and 97.5% quantiles of the simulated values across individuals were used to define the range of slopes. The proportion of negative slopes within-individual simulations was used to estimate the percentage of individuals showing advanced departure dates.
We calculated repeatability (intra-class correlation coefficient, r ± SE)68 of departure date from New Zealand for 124 individuals observed 3–13 years each.
Geolocator tracking and analysis
In 2008–2009 and 2013–2014, a subset of individuals was tracked for the full migration to Alaska using light-level geolocators. The units (1.0–1.5 g) were mounted to a colorband on the tibia; attachment plus all individual markings represented 1–2% of lean body mass. We deployed geolocators in March 2008 (n = 17 MK14, British Antarctic Survey, UK), October 2008 (n = 19 MK14), February 2013 (n = 15 MK4093, Biotrack, UK, and n = 24 Intigeo-C65K, Migrate Technology, UK), and November 2013 (n = 22 Intigeo-C65K) and recaptured birds during the following non-breeding season(s) to retrieve the units. This resulted in the following number of complete northward tracks (including breeding location; see below): 12 in 2008, 13 in 2009, 13 in 2013, and 12 in 2014; 25 individuals were tracked in 1 year, 11 in 2 years, and 1 in 3 years.
Godwits depart New Zealand near the austral autumnal equinox, when light-level geolocation is least effective for describing movements, particularly for flights in a north-south direction69. However, geolocators also recorded conductivity (i.e., contact with salt water); wet-dry transitions can thus identify extended periods of flight between intertidal habitats with high precision64. From conductivity data, we derived four parameters for each northward track: dates of departure from New Zealand (NZdep), arrival and departure from the Yellow Sea region of Asia (YSarr, YSdep), and first arrival in Alaska (AKarr). Most godwits spend up to two weeks in coastal SW Alaska before moving to breeding sites22,24, but some individuals show no discernible shift in light or conductivity data during this period. Therefore, we did not analyze arrival at the ultimate breeding site. We derived a fifth timing parameter, duration of the stopover in the Yellow Sea (YSdur), as the difference (in days) between YSarr and YSdep.
Because the geolocators were leg-mounted, light-level data in Alaska also indicated periods of nest incubation, when the unit was shaded by the sitting bird24,27. During the breeding season, geolocators registered nights as regular, clearly demarcated periods of darkness 0–4.5 h in length; these did not appear at all if birds were north of 64°N. Days appeared as continuous light, irregularly broken by a brief (<1 h) shading events, most likely corresponding to behaviors such as wading or sitting. Within 6–25 days of apparent arrival on breeding grounds, most godwits (41 of 50 total tracks) displayed a conspicuous pattern of incubation, in which semi-regular shading events of 4–13 h were overlaid on the day/night pattern for periods up to 25 d. Bar-tailed godwits incubate tundra nests bi-parentally, and so both sexes demonstrate this pattern. We considered the first day of this period to be the start of incubation (IncSt), and analyzed this timing parameter along with the five migration parameters (above).
For godwits tracked in 2008–2009 using MK14 geolocators, breeding locations are published previously24 (Supplementary Fig. 3a). We used BASTrak70 software (British Antarctic Survey, UK) to produce twice-daily location estimates, and calculated mean latitude and longitude (plus 95% range of estimates) during stationary periods in the breeding season, after removing light-dark transitions clearly affected by behavioral shading events such as incubation. Based on periods when the bird’s true location was known (at the deployment site in New Zealand, and in some cases additionally in Asia), we used sun angles of 3.0–3.6° below the horizon to represent twilight for calibration of light-level data for location estimates. For three individuals that traveled north of the Arctic Circle, narrow-range light sensitivity of MK14 geolocators precluded location estimation in 24-h daylight; however, longitude estimates immediately before and after the breeding season were east of 160°W, strongly indicating breeding destinations on the central North Slope of Alaska. Therefore, we assumed a breeding latitude of 70.2°N (the mid-point of the known breeding range in this region; see Supplementary Fig. 3a) for these individuals. For individuals tracked in multiple years (n = 6), the location derived from the year with the best available data (i.e., longer uninterrupted stationary breeding period) is shown in Supplementary Fig. 3a.
For godwits tracked in 2013–2014 using Intigeo-C65K geolocators, which recorded full-range light levels, we estimated breeding locations using the R package PolarGeolocation v.0.1.071. To calibrate the maximum light curve and the error distribution of light recordings from the maximum light curve, we used the longest possible period for each bird, i.e., all days when the bird was at the Manawatu River estuary (range 38–161 days; mean = 127 days). Next, we selected a time period during the breeding season that included a clear pattern of breeding (incubation inferred from shading patterns); although including incubation increases uncertainty in the breeding-site estimate, we chose such periods to ensure that we did not include periods spent away from the breeding site in the location estimation. The duration of the period used to estimate breeding locations was 12–48 days (mean = 28 days). PolarGeolocation estimates a gridded likelihood surface from which we can derive confidence levels of our best location estimate; these confidence intervals depend on the extent of the grid. To make the confidence estimates comparable across individuals, we ran a first estimate using a grid with 50 km resolution and a radius of 1500 km around the approximate center of the breeding range (162°W, 62°N). Next, we re-centered the mask around the best estimate and re-ran the simulation using the same radius of 1500 km, and the same likelihood contour level of 0.05 to describe the confidence intervals for all birds (Supplementary Fig. 3b). For individuals tracked in multiple years (n = 4), the location derived from the first year of tracking is shown in Supplementary Fig. 3b.
To estimate breeding sites from four Biotrack MK4093 geolocators in 2013, we performed the simple threshold approach. First, we defined sunrise and sunset times using a light intensity threshold of 2 arbitrary units, and used the periods the individual bird was at the deployment site in New Zealand for calibration and calculation of a reference zenith angle69 that was then used to estimate locations via the thresholdPath function in the R Package SGAT v.0.1.372. All loggers showed dark periods during the breeding season (i.e., remained south of the Arctic Circle), allowing estimation of locations throughout the year. We extracted the locations for the period including incubation and calculated the median location and the 95% credibility interval.
Bar-tailed godwits breeding in northern and southern Alaska differ in size and migration timing24,25, and potentially experience very different temporal trends in conditions during the breeding season. Therefore, tracked individuals were grouped into two regional groups, comprising those breeding >64°N (North) or <64°N (South). Individuals tracked in multiple years and/or with different geolocator types allow assessment of repeatability of breeding location estimates; however, we note that the uncertainty measures we show in Supplementary Fig. 3 for the three estimation methods are not strictly comparable. For six individuals with breeding locations derived in both 2008 and 2009 from MK14 geolocators, location estimates differed by 24–124 km, respectively. For four individuals tracked with Intigeo units in both 2013 and 2014, location estimates differed by 50–190 km, respectively. For one individual tracked with both MK14 and Intigeo loggers, the best estimates for 2008 and 2014 differed by approximately 78 km. These differences fall within the expected error of geolocation, and in no case affected an individual’s assignment to the northern or southern groups. The final set of tracked godwits with known breeding latitude included 36 individuals (North = 16, South = 20; Fig. 2a and Supplementary Fig. 3).
To describe the rate of change across 2008–2014 in six migration and breeding phenology parameters (NZdep, YSarr, YSdep, AKarr, IncSt, and YSdur; see above), we used linear mixed models including individual as a random effect and breeding region (North or South) as a fixed factor (Fig. 4, Supplementary Table 3). To account for unequal samples of northern and southern breeders across years, we initially included the interaction of year*breeding region in the models; this interaction factor was non-significant for all phenology parameters (p = 0.08–0.85), and so we excluded it from the final models.
Environmental data and analysis
In Arctic-breeding shorebirds, the timing of breeding closely follows the retreat of snow-cover from tundra nest sites in the spring73. We compared two indices of breeding phenology (timing of snowmelt and spring green-up), summarized separately for the two regions (North, South) of the Alaska breeding range of bar-tailed godwits. We did this for two temporal periods: one encompassing the entire study period (2008–2020), and a shorter-term related directly to the period in which we tracked individuals with geolocators (2008–2014).
Snowmelt
Remotely sensed IMS Daily Northern Hemisphere Snow and Ice Analysis data for the period 2008–2020 on a scale of 4 × 4 km were downloaded74. Next, grid cells falling in the subspecies’ breeding range were extracted. Then, data from cells that did not include measures of sea ice were modeled using a maximum likelihood fit (mle2 function from R package bbmle v.1.0.23.1) of the asymmetric Gaussian model function75 with a binomial error distribution. The model fit is a compromise between the first snow-free day and the day when the pixel remained snow-free in late spring and summer (for more details and comparison between different snow-melt indices see Supplementary Methods). Using these year-specific fits, the date of snowmelt was then determined as the date on which the fitted curve predicted 1/3 of the cell area to be snow-free.
Normalized Difference Vegetation Index (NDVI)
From the NOAA STAR Vegetation Health Product data set, weekly noise-removed NDVI data on a scale of 4 × 4 km grid cells were downloaded for the period 2008–202076. Next, cells located in the subspecies’ breeding range were selected for further analysis. To identify the annual day when the vegetation starts to grow (green-up day), a penalized cubic smoothing spline was fitted to the yearly NDVI values, allowing daily interpolation and the calculation of when the curve exceeds 15% of the greatest amplitude during spring77. To ensure a robust fit of the smoothing spline and to remove potential artifacts of single pixels, the NDVI values of the surrounding pixels within a radius of 15 km were taken into account. However, the distance to the focal cell was used to weigh the importance of the model fit (with a Gaussian decline over distance with a standard deviation of 4 km). Because the presence of snow and ice can significantly affect NDVI values, the model was used to interpolate daily NDVI values only for days that were snow-free in at least one of the pixels within the 15 km radius. Next, the day of the maximum NDVI was established, and the day when the increasing NDVI curve (start to maximum) exceeded 15% of the amplitude was extracted as the NDVI start date. More detailed information and a comparison between the estimated NDVI start date and the higher spatial resolution MODIS derived green-up day (also based on a 15% threshold) are in Supplementary Methods.
Snowmelt and NDVI start were highly correlated across the breeding range (t160221 = 664.95, p < 0.001, R2 = 0.856) as well as in the northern (t120852 = 494.19, p < 0.001, R2 = 0.817) and southern (t39367 = 285.65, p < 0.001, R2 = 0.821) regions separately.
Ethics statement
Fieldwork was conducted with Massey University Animal Ethics Committee approval (#07/163, 12/90, and 16/117) and appropriate New Zealand Department of Conservation permits (Banding permit 2007/39, 35503-FAU, 38111-FAU).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
