In the current, warming climate many organisms in seasonal environments advance their timing of reproduction to benefit from resource peaks earlier in spring. For migrants, the potential to advance reproduction may be constrained by their migration strategies, notably their ability to advance arrival on the breeding grounds. Recent studies show various changes in migration strategies, including wintering closer to the breeding grounds, earlier departure from the wintering grounds or faster travels by spending less time on stopover sites. However, whether such changes lead to earlier arrival or earlier breeding remain open questions. We studied changes in migration and reproduction timing in 12 populations of nine migratory birds, including seabirds, shorebirds, birds of prey and waterfowl breeding on Arctic sites bordering the Greenland and Barents Sea, a region undergoing rapid climate warming. Timing of migration and reproduction was derived from tracking and field data, and analysed to study (1) how timing has changed in response to changing moment of snowmelt at the breeding grounds, and (2) what adjustments in migration strategies this involved. We found that in years with early snowmelt egg-laying in multiple populations advanced, but only two waterfowl populations also advanced arrival in the Arctic. In contrast, arrival in the Arctic generally advanced with time, even when snowmelt or egg-laying dates did not advance. Earlier arrival with time was mostly explained by populations travelling to the Arctic faster, likely spending less time at stopover sites. Inability to forecast conditions in the Arctic may limit birds to adjust migration timing to annually varying snowmelt, but we show that several species, particularly waterfowl, are able to travel faster and advance timing of migration over the years. The question remains whether this reflects adaptations to Arctic climate change, or other factors, for example environmental changes along the migratory route.

We gathered data on the timing of spring migration and egg-laying to local snowmelt timing for a suite of arctic-breeding bird species, using a large set of tracking data supplemented by field data on egg-laying dates. Rather than using tracking data of individuals tracked from the non-breeding grounds and breeding across the entire breeding range of a species, the focus was placed on tracking data of individuals breeding at specific Arctic breeding sites. In this way, we aimed to reduce variation in the timing of arrival in the Arctic caused by variation in breeding location (e.g., Conklin et al., 2010).

Study species and sites

Our dataset contained tracking data (519 spring migrations of 274 individual birds) from 12 populations and nine migratory bird species, breeding at seven Arctic study sites in the Barents and Greenland Sea region from the period 2007 - 2023. For this paper, we treat each species at a separate study site as a separate population. The species and study sites included greater white-fronted goose Anser albifrons (Kolguyev Island, north-western Russia), pink-footed goose Anser brachyrhynchus (Adventdalen, Svalbard), barnacle goose Branta leucopsis (Kolguyev Island, north-western Russia and Kongsfjorden, Svalbard), tundra swan Cygnus colombianus (Malozemelskaya tundra, north-western Russia), red-necked phalarope Phalaropus lobatus (Ammarnäs, northern Sweden and Slettnes, northern Norway), sanderling Calidris alba (Zackenberg, north-eastern Greenland), Arctic skua Stercorarius parasiticus (Slettnes, northern Norway and Kongsfjorden, Svalbard), long-tailed skua Stercorarius longicaudus (Ammarnäs, northern Sweden) and rough-legged buzzard (Kolguyev Island, north-western Russia),. Birds were tracked using GPS-loggers, GPS-GSM-transmitters, satellite transmitters (GPS-PTT) or geolocators (Global Location Sensing, GLS). We tracked only females of the three goose species and both females and males of the other species. Most birds were captured and fitted with tracking devices on their breeding grounds, and the study sites of these populations, defined as the area in which birds were captured, were relatively small (8 – 56 km²). Only greater white-fronted geese, barnacle geese (north-western Russian population) and tundra swans were captured and tracked from their non-breeding grounds in the Netherlands and Germany. For these populations, we selected larger areas as study sites (6,061 – 15,132 km²) known to host major concentrations of these species (Glazov et al., 2021; Kondratyev et al., 2013; Rees, 2010), namely Kolguyev Island (greater white-fronted geese, barnacle geese) and the Malozemelskaya tundra region, west of the Pechora Sea (tundra swans). Only individuals that initiated a nest in these study sites were included (see methods section ‘Date of egg laying’). Whether or not tundra swans showed nesting behaviour was more difficult to detect from tracking data, and we developed a more comprehensive method to evaluate this, making additional use of acceleration data (supplemental materials). To increase the sample size for this species, we also included male individual swans that stayed in the Malozemelskaya tundra study site for the entire month of June. Rough-legged buzzards were captured and tracked from various locations on Kolguyev Island (Pokrovsky et al. 2024), and we therefore used the study area of Kolguyev Island also for this population. Detailed information on tracking methods can be found in the species-specific studies.

Tracking data preprocessing

Tracking data were acquired from Movebank (Kays et al. 2022) and datasets from the authors. For Arctic skuas, long-tailed skuas, red-necked phalaropes, and sanderlings, we used previously obtained tracks from GLS data as described in (van Bemmelen et al., 2017; Reneerkens et al., 2019; van Bemmelen et al., 2019; van Bemmelen et al., 2024). All studies used a light intensity threshold of 2.0 to determine twilight events. To avoid large errors around the equinox, latitude data were deleted from 14 days before 18 days after the spring equinox. For greater white-fronted geese (Kölzsch et al., 2016), barnacle geese (Boom, Lameris, et al., 2023; Kölzsch et al., 2015), pink-footed geese (Schreven et al., 2021), tundra swans (Linssen et al., 2023; R. Nuijten & Nolet, 2020) and rough-legged buzzards (Pokrovsky et al., 2021), previously published GPS-tracks were downloaded from Movebank. GPS-tracks were filtered to exclude outliers, identified as positions where the speed required to travel from that to the next positions was larger than 120 km/h (which is 10 km/h faster than recorded flight speeds for these species, (Alerstam et al., 2007; Miller et al., 2005)). We only included GPS- and geolocator-tracks that contained at least one position per day from January until arrival at the Arctic Circle.

Migration timing, migration distance, and travel speed

From tracking data, we extracted individual data on spring migration timing, including departure from, non-breeding area and arrival in the breeding area, as well as total migration distance and travel speed. First, individual non-breeding sites were calculated as mid-winter geographical centroids: the median position during January and February, a period during which all populations are considered to be present in their non-breeding regions. The calculation of departure date from the non-breeding site differed between geolocation- and GPS-based tracks. For GPS tracks, we defined the departure date as the date with the last position in spring within 200 km of the mid-winter centroid. For GLS-tracks, which suffer from a larger positional error, especially in latitude during the period around the equinox (Lisovski et al., 2012), we closely inspected raw position estimates to find the first date with consistent directional movement away from the non-breeding site.

Arrival date in the Arctic was determined as the date at which birds crossed the Arctic Circle (latitude 66.33°N). For GPS data, this was defined as the first date above 66.33°N. GLS-tracks are calculated based on sunrise and sunset events, and therefore cannot be calculated for positions above the Arctic Circle during the polar day with 24 hours of continuous light (Lisovski, 2018). For these tracks, we defined crossing of the Arctic Circle as the first day after the date of the last night (which usually occurred at c. 60˚ N latitude). An earlier study found that the crossing from these last positions to the Arctic is generally rapid without stops (van Bemmelen et al. 2024).

Total migration distance was calculated as the cumulative sum of great circle distances between daily averaged positions, from departure from the non-breeding site to arrival in the Arctic (at the Arctic Circle). Travel speed was calculated as the total migration distance divided by the time between departure and arrival in the Arctic in km/day. Due to high flight speeds, the actual time spent flying is relatively short, therefore, variation in migration speed mostly represents variation in the amount of time spent on stopover sites along the migratory routes (Alerstam & Bäckman, 2018). We emphasise that our measure of travel speed is not the same as migration speed, which includes the time spent fuelling at non-breeding sites before departure (Alerstam, 2003).

Date of egg-laying

Date of egg-laying, defined as the date at which the first egg was laid, was determined from (i) field observations, (ii) patterns in light measurements, and (iii) from GPS tracking data. (i) In field observations, date of egg-laying for individuals with and without tracking devices were either observed directly, back-calculated from incomplete clutches in geese (van der Jeugd et al., 2009) and shorebirds (Liebezeit et al., 2014; Reneerkens et al., 2016), back-calculated from observed hatching dates (by subtracting the period required for incubation, based on (Cramp & Perrins, 1988) or assessed by egg floatation (Liebezeit et al., 2007). (ii) For GLS-data, laying dates were determined based on regular periods of darkness of at least 1h while the bird was in an area with continuous daylight (Verhoeven et al., 2020). (iii) For GPS-data of geese (of which only females were tagged), timing and position of nesting was determined from GPS-locations after arrival on the breeding site, where the nest location was determined as a position where the daily standard deviation in latitude was less than 25.4 m (Schreven et al. 2021) for at least three consecutive days. The first of these three days was then determined as the date of egg-laying. For greater white-fronted geese and barnacle geese from Kolguyev Island, for which both population-average laying dates from field data as well as extracted from tracking data were available, yearly averaged laying date correlated strongly between both methods (Pearson’s correlation = 0.83, t  = 4.38, p = 0.002, n = 9 years). In contrast to geese, where only the female incubates the eggs, in tundra swans, females and males take turns to sit on the nest. This makes it more difficult to distinguish nest positions, and we therefore used a modified method that also used accelerometer data were also used (for details, see Supplementary methods). For GPS data of rough-legged buzzards, we determined the start of nesting as the first day in a period for which a bird stayed within one position, i.e., a difference between GPS positions of less than 3 m, for more than 1 day. This method was verified using accelerometer data (for details, see Curk et al., 2022). For long-tailed skuas, red-necked phalaropes, and sanderlings, laying dates could not (with a few exceptions) be determined for tracked individuals, as geolocators stopped functioning before breeding, or birds did not breed. To analyse trends in laying dates for these populations, we supplemented our dataset with laying dates determined from field data for the same study sites and years.

Date of snowmelt

Snow cover for all study sites was estimated using satellite images of the 500m resolution MODIS Terra Surface Reflectance Daily Global product (MOD09GA, v6.1, Vermote & Wolfe, 2021). The analysis was conducted in Google Earth Engine (Gorelick et al., 2017), using the R package RGEE (Aybar et al., 2020) and the automated workflows for the quantification of snowmelt by Versluijs (2024).

We manually mapped spatial overlays of study sites. Although sites differed strongly in size, all sites represented areas of mostly homogeneous elevation with little expected variation in snowmelt, thus giving an average across each of the study sites. We extracted satellite images between March 15^th and September 15^th for each study site, covering the years 2000 - 2023. Within each extracted image, we masked pixels that were classified as clouds by the 1000m MOD35 cloud mask product (MODIS Atmosphere Science Team, 2012), and as waterbodies (i.e. oceans, rivers, lakes and large ponds) by the 250m Terra Land Water Mask dataset (MOD44W, v6.0, (Carroll et al., 2017). We subsequently calculated the normalized difference snow index (i.e., NDSI; Dietz et al., 2012; Dozier, 1989) value per pixel, as the normalized difference between the green and the short-wave infrared band:

For each satellite image and within each study site, we then calculated the percentage of snow cover as the fraction of pixels with NDSI values larger than 0.4 (Dozier, 1989; Hall et al., 1995). We fit a general additive model (GAM; Wood, S., 2023) to these annual time series of snow cover and extracted the moment the GAM first dropped below 50% snow cover, which we used in subsequent analyses as the date of snowmelt.

Statistical analyses

For analyses, we used generalized linear mixed models (GLMMs, using the R-package ‘lme4’(Bates et al., 2012)) were used including a random intercept for individual birds in models with year as an independent variable, and random intercepts for individual birds and year in all other models. The year and date of snowmelt were centered by subtracting population-specific means. As strong interannual variation in environmental conditions (incl. date of snowmelt) could distort temporal trends, we excluded short time series (less than 5 years) from analyses of trends over time (i.e., over the years). Below, the steps of the full analysis are described in detail:

First, we analysed trends in date of snowmelt over time (using the specific time period for which tracking data for the population was available) using linear models (LMs) including date of snowmelt as dependent variable and year, population, and their interaction as independent variables (i).
(i) date of snowmelt (D_sm) = intercept (intrcpt) + year (Y) + population (P) + P x Y

Thereafter, we analysed the relationships of arrival date and laying date with date of snowmelt (ii, iv) and trends over time (iii, v). We used GLMMs for all species combined, including arrival date or laying date as a dependent variable and either date of snowmelt, population, and their interaction, or year, population, and their interaction as independent variables. From these LMMs, population-specific trends were inspected post-hoc using confidence intervals in the ‘emtrends’ function in the package ‘emmeans’ (Lenth 2017).
(ii) arrival date (Da) = D_sm  + P + D_sm x P + (1|ID) + (1|Y)
(iii) Da = Y + P + Y x P + (1|ID)
(iv) laying date (Dl) = D_sm + P + D_sm  x P + (1|ID) + (1|Y)
(v) Dl = Y + P + Y x P + (1|ID)

Subsequently, we analysed whether population-level relationships between arrival date and snowmelt, as well as population-level trends over time, could be explained by similar changes in departure date, migration distance, and travel speed about snowmelt or over time. To do so, we extracted slopes from population-specific GLMMs (vi) with arrival date, departure date, migration distance, or travel speed as dependent variables and date of snowmelt as an independent variable, including random intercepts as described above. In these GLMMs, arrival date, departure date, migration distance, and travel speed were standardized by subtracting the overall mean and dividing by the standard deviation to allow for later comparison of slopes. We then fitted LMs (vii) for all populations combined with arrival date slope as dependent variable and departure date, migration distance, and travel speed slopes as independent variables. We followed the same procedure for population-level trends with time, extracting slopes from population-specific GLMMs with time as an independent variable (viii) and fitting LMs for all populations combined (ix).
(vi) Da/departure date (Dd) / migration distance (MD) / travel speed (V) = D_sm + (1|ID) + (1|Y)
(vii) slope arrival date ~ date of snowmelt (Da_sm) = intrcpt + slope departure date (Dd_sm) + slope migration distance (MD_sm) + slope travel speed (V_sm)
(viii) Da / Dd / MD / V = Y + (1|ID)
(ix) slope arrival date ~ year (Da_y) = intrcpt + Dd_y + MD_y + V_y

Finally, we analysed whether variation in laying dates between individuals and among years could be explained either by local date of snowmelt or individual arrival date, using GLMMs (x) including date of egg laying as dependent variable and date of snowmelt, arrival date and population, as well as interactions date of snowmelt x population and arrival date x population, as independent variables, along with random intercepts described above. In these analyses, laying and arrival dates were centred by subtracting population-specific means.
(x) Dl = Dsm + Da + P + Dsm x P + Da x P + (1|Y) + (1|P)

We used AICc values to compare the performance of models, including all possible combinations of independent variables, including an intercept-only model. For most analyses, this meant comparing the performance of models including either date of snowmelt or year as independent variables with intercept-only models. We selected the model with the lowest AICc, but models within 2 ΔAICc of the top supported model were also considered informative, except when containing additional parameters (Arnold, 2010). All analyses were conducted using R (version 4.3.1).

Data from literature

We collated papers that reported on arrival and laying date relative to the date of snowmelt in the Web of Science database, conducted on 7 February 2024. We used the search term: bird AND Arctic AND (snowmelt OR snow cover OR snow) AND (migration timing OR arrival OR laying date OR nest initiation OR egg) using the option ‘all fields’. This resulted in 213 papers, which were scanned for statistics and data on the timing of migratory arrival (on northern stopover sites or breeding sites), measured for at least three years, concerning any metric of snow cover. We omitted studies where it was not possible to extract slopes, arrival, laying, or snowmelt dates from either the text, tables, or graphs. Data were found in six papers (Table S2). Three papers (which were not identified in the search) known to contain similar data were added (Ely et al., 2018; Hupp et al., 2018; Lameris et al., 2018). Unpublished data for three additional species were supplemented by S. Volkov, with the data provided in the dataset linked to this paper (Lameris et al. 2025) and details on data collection in the supplemental material. Among all considered studies, the snow cover fraction used or extracted as a measure of ‘date of snowmelt’ varied between 60% and 0%. In general, different fractions of snow cover will be highly correlated (e.g., 0.99 Pearson’s correlation between 0.25 and 0.75 fraction of snow cover in our dataset), although the moment of 0% snow cover may show less of a correlation with mid-point snow cover. In total, the dataset from literature contained 18 species (five species of waterfowl, 10 shorebirds, two passerines, and one gull) from six different locations in the Holarctic. We extracted slopes of the change in arrival dates, laying dates, and date of snowmelt over time, and slopes of the change in arrival dates and laying dates about date of snowmelt. We then used t-tests to analyse whether the slope in change in arrival date and laying date, with time as well as concerning date of snowmelt, differed from zero. In addition, we used LMs to analyse whether changes in arrival and laying date over time correlated with changes in date of snowmelt over time, comparing models with either slopes of arrival or laying date as dependent variables and the slope of date of snowmelt as independent variables, and performing model selection as described above.