Data on environmental conditions, breeding phenology and breeding success were collected in 2019 and 2020 during the breeding season, from May to August, in Parco Naturale Val Troncea, 44°57′28″N, 6°56′28″E, Western Alps, Piedmont, Italy (Figure 1A). The elevation of the study area varied from 1,560 m a.s.l. to 2,700 m a.s.l. and was located on a south-eastern slope. The landscape structure is characterised by forests (predominantly larch Larix decidua) with shrubby patches at lower elevations, and larches interspersed with open grasslands at around 2,200 m a.s.l., followed by rocky areas, above the treeline. The valley is used seasonally as pasture, and grazing cows influence the distribution and availability of potential open grassland habitat across all elevations (Jähnig et al., 2018).

To investigate the longer-term trends in the annual date of snowmelt (“the snow free date”), we used the remotely sensed MODIS (MOD10A1.006) Terra Snow Cover Daily Global 500 m dataset with recordings starting in 2000 (Hall et al., 2016). Daily values for all pixels (N = 68) within the study area were extracted from Google Earth Engine (Gorelick et al., 2017) via the R package rgee (Aybar, 2020). For each pixel and each year, we identified the snow free date as the first date without snow that was followed by at least 60 days (approximately 2 months) of snow free conditions, ignoring potential snow free periods within the winter.

To derive higher-resolution snow free dates for the actual nest sites in the study area during 2019 and 2020, we used the remotely sensed Sentinel-2 MSI MultiSpectral Instrument Level-2A dataset with a spatial resolution of 10 m (Drusch et al., 2012). On the Google Earth Engine Server, we first selected all images taken in 2019 and 2020 intersecting at least one of the nest sites. Next, we masked cloud pixels using the Sentinel-2 Cloud Probability layer with the standard probability threshold of 60. We then calculated the normalised differenced snow index (NDSI, Band 8 and 4) and the normalised vegetation index (NDVI, Band 3 and 11). For each nest site and each image, we extracted the mean NDSI and NDVI values of all pixels within a 10 m radius (N = 2–4) around the nest. For each nest site and for both years, we fitted a smooth line (loess fit with span = 0.2 from the R package stats) to the NDVI and NDSI values over the day of the year, and extracted the day when the model fit fell below 0.1 NDSI (less than 10% of snow cover was defined as the snow free date) and exceeded 80% of the NDVI annual amplitude (the defined threshold when the peak in spring green-up was reached).

In the field, snow cover and grass height, as proxies for the spring green-up were monitored throughout the two breeding seasons (see Supplementary Table 1 for details and sample sizes). Data were collected using a standardised protocol for territory mapping. In summary, we defined four sectors above the treeline (total area: 161 ha) that were characterised by open grassland habitat, variation in elevation (∼ 700 m) and a high density of wheatear territories. We defined a route (∼ 6 km in each sector) along which we defined habitat survey points separated by a distance of 200 m (6–8 points each sector, 30 in total). Territory mapping and habitat point surveys were conducted over five visits between 27th May and 21st July 2019. Positions of singing males were recorded with a GPS device during the period between sunrise and 12 am and during suitable weather conditions. Following the protocol introduced by Südbeck and Weick (2005), territories were confirmed when a singing male or other territorial behaviour was recorded during at least two visits at the same point. Habitat parameters were recorded at each of these points (covering an area of 200 m × 200 m at habitat survey points and a 100 m radius around territory points) at each of the five visits and were all included in analyses to track changes in environmental conditions across the season. The territory mapping as well as the habitat point survey was repeated in 2020 over five visits between 14th May and 6th July.

Males and females were trapped either at their nest during chick feeding or with spring traps (using mealworm bait) when the nest location was unknown, not yet established, or chicks had fledged. In 2019, 35 adults were tagged, using 15 (eight males, seven females) light-level geolocators manufactured by Migrate Technology Ltd. [Intigeo P65, light stalk, 0.8 g, fitted with a leg-loop harness (Rappole and Tripton, 1991) and 20 (10 males, 10 females) geolocators manufactured by the Swiss Ornithological Institute (SOI, SOI-PAM, light stalk, recording pressure, 1.1 g)]. Device weights were on average 4.1% of the individuals’ body mass (2.8–5.2%, see Supplementary Table 2). All captured individuals were marked with three plastic colour rings and a metal ring of the Italian bird ringing scheme. As a control for the effect of geolocator devices on return rates, an additional 16 individuals were ringed, but were not fitted with a geolocator. All animal handling and protocols were carried out in accordance with relevant guidelines and regulations under licences issued by Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA, protocol nos. 27303 and 2463) and Città di Torino (licence nos. 56.1433/2019 and 197.4765/2019).

Light-intensity data were recorded at 5 min intervals and analysed using a threshold method (Lisovski et al., 2020a). Sunrise and sunset events (twilight events) were identified on log-transformed light data and a threshold of 1 log lux, using the R package TwGeos (Lisovski et al., 2015).

Twilight events recorded at the known deployment site during periods after deployment and before recapture were used for calibration, i.e., to estimate the error distribution of twilight events and the individual reference sun elevation angle (position of the sun when twilight events were detected). For recordings from Migrate Technology Ltd., tags, migratory movements were identified as sudden and directed changes of consecutive sunrise or sunset events (for details, see the “invMovement” method description in Lisovski et al., 2020b). Stationary periods with a minimum of two consecutive twilight events were then identified as periods surrounded by migratory movements. For data recordings from Swiss Ornithological Institute tags, we used the additional air pressure measurements to identify movement periods with a change in a pressure threshold of 4 (change in the altitudinal position of the bird), using the method implemented and described in the R package PAMLr (Dhanjal-Adams, 2019).

For the final track estimation, we used a Bayesian approach from the R Package SGAT (Wotherspoon et al., 2016) allowing the incorporation of twilight events, their error distribution (gamma density distribution), the information on periods of movement and residency, a spatial probability mask and expected flight speed distribution (gamma distribution with shape = 2.2 and rate = 0.08).

The method provides the most likely paths with credible intervals derived via Markov Chain Monte Carlo (MCMC) simulations. The applied groupedThresholdModel provides a single location estimate for stationary periods as well as one location estimate for each twilight during movement. A spatial mask was applied for the estimate of stationary locations, restricting positions to be estimated on land. During the MCMC simulation, the first and last site of residency (in case the logger was still recording light on return from the wintering grounds) were fixed to the deployment site.

We first ran a modifiedGamma model (relaxed assumptions) for 1,000 iterations to initiate the model, before tuning the model with final assumptions/priors with three runs, each containing 300 iterations. Finally, the model was run for 2,000 iterations to ensure convergence. Median location estimates and 95% CrI were calculated using the entire final MCMC chain (i.e., each location estimate was based on 2,000 estimates from within the MCMC chains). Locations from the most likely track were used for plotting.

We used three different approaches to estimate the timing of departure and arrival at the breeding grounds in 2020 from the geolocator data. First, we used the dates from the above-described analysis of the stationary periods (first day at breeding ground after last northward migration stopover site, Supplementary Table 3). Second, we performed a separate analysis on the defined twilight times to estimate the probability that each sunrise and sunset time were recorded within the larger area around the breeding site (radius = 250 km). To this end, we used the calibrated individual gamma distribution of the twilight error and simulated 500 locations for each twilight time. The resulting locations reflect the expected precision of the location estimates (see Lisovski et al., 2018). The probability was then defined as the percentage of location estimates per twilight falling into the defined radius around the breeding site. A symmetrical Gaussian curve was then fitted to the probabilities over time, and departure and arrival were calculated as the dates when the curve fell, respectively, below or exceeded a probability of 0.5 (as well as the 0.025 and 0.975 probability to report the 95% confidence intervals). The idea behind the third approach emerged while inspecting the light recordings. During the sunrise periods at the beginning and the end of the recordings, thus at the breeding site, we discovered a time delay in the daily maximum light values. This was caused by the south-eastern slope of the breeding area and the fact that the direct sunlight was blocked by the mountain ridge in the east (Figure 1A). To make use of this pattern in the identification of the exact arrival and departure time at the breeding site, we applied an additional calibration with a higher light intensity threshold (9 log lux), slightly below the maximum light values that the tags can resolve (11 log lux), resulting in a gamma distribution of the delayed maximum light intensities. Due to the changes in the position of the sun, the delay was expected to change over the course of the year, thus we estimated the gamma distribution parameters for autumn and spring separately. Next, we calculated the probability that the daily timing of the maximum light intensity during sunrise would belong to the derived gamma distribution (the first half of the light recordings were compared to the autumn gamma distribution and the second half to the spring gamma distribution). A symmetrical Gaussian curve was then fitted to the probabilities over time, and departure and arrival was calculated as the date when the curve fell, respectively, below or exceeded a probability of 0.5 (as well as the 0.025 and 0.975 probability to report the 95% confidence intervals).

To estimate the general migration transition of wheatears migrating northward from African wintering ranges into the regions of Piedmont and Aosta Valley, and thus elaborate a proxy for the arrival in the wider area of the breeding site in the Western Alps in 2019 and 2020, we additionally used citizen science data of wheatear sightings during spring submitted via the ornitho.it platform. The data includes the number of all observations of wheatears per day in the period March–May. We furthermore used data from the years 2010 to 2020 to describe changes over a longer period, being aware of the bias in the data that might include not only local breeding wheatears, but also individuals migrating further north. It is not possible to formally account for observer effort in the ornito.it data, as observers often submit occasional observations of specific species. However, based on both the total number of observations of any species submitted in each year (2019: n = 151.799, 2020: n = 168.180) and the number of individual wheatear observations (2019: n = 101, 2020: n = 93), there appears to be negligible differences between years, hence we consider the data appropriate to indicate the general temporal migration patterns in the larger area of the study site.

From the fourth week of May until the first week of August (2019, 2020), we searched for wheatear nests in all four sectors as well as adjacent accessible areas, including open grassland sites at lower elevations (pastures). Each identified nest was marked with bamboo sticks at a distance of at least 5 m, and consecutively monitored until the chicks fledged or the nest failed.

Depending on the stage (e.g., longer breaks during incubation when date of hatching could be estimated from date of laying the last egg, following Conder, 1989), the status of the nests was recorded on average every 5 days. In 2020, at each visit, habitat parameters at the nest were recorded. We were able to record the exact dates for building, laying, hatching, and fledging for some nests in 2019 and 2020, and subsequently estimate the duration of each stage for our population (see for details Supplementary Table 4). We used the average of these individual durations of each stage (Supplementary Table 4) to estimate breeding phenology dates for nests where direct observation was not possible in all their stages (back or forward in time from a known breeding phenology date). If the number of eggs laid was unknown, we assumed a clutch size of five eggs, which is the average number within our population (nests reaching the stage of incubation and no. of eggs known, 2019: 36, 2020: 22), to estimate the lay date. Breeding phenology dates were calculated only for stages reached by the individual nest, e.g., we did not estimate a fledge date if the nest failed before. The arrival-breeding interval (A-B interval) was estimated by calculating the difference between the median arrival date and lay date, providing two estimates of A-B interval for 2020 (using geolocation and citizen science observations).

To estimate and compare the breeding success of the population between the two consecutive years 2019 and 2020, we considered three different measures: (a) the hatching success rate (eggs that hatched), (b) the fledging success rate (hatchlings that fledged), and (c) the nest success rate (successful nests/all nest attempts). Measures a and b allowed us to separate the breeding success into two states of a nest, i.e., nests during incubation and nests during chick rearing, by considering the success rates of individual eggs and hatchlings estimated as probabilities of reaching the next state (“period” hatching and fledging). We excluded nests where number of eggs or chicks was unknown or incomplete. For the vast majority of the nests, we could not determine whether this was a first attempt or a replacement of a failed attempt due to predation, abandonment or hatching failure. Thus, we included all nests monitored in our analyses, except for the nests known to be replacements.

We applied a general linear mixed model with binomial errors to quantify the success rate of hatching and fledging (a, b) of each egg laid within the population, with the fixed effects “year” (2019, 2020) and “period” (hatching, fledging), adding “nest id” as a random effect [using the function glmer from package lme4 (Bates et al., 2007)]. To quantify the nest success rate (c), we applied a generalised linear model with the fixed effect “year” (using the function glm from package stats). Model assumptions were tested by deriving residual plots and checking for overdispersion [function dispersion_glmer from package blmeco (Korner-Nievergelt et al., 2019)]. Conclusions from test statistics were based on the credibility intervals from the model predictions calculated using the function sim from package arm (Gelman and Hill, 2007), and by drawing 2,000 random samples from the posterior distribution of the model parameters to obtain their 95% credible intervals from these simulations. An effect was considered significant if the credible interval of the model estimate did not include zero (Korner-Nievergelt et al., 2015).

The time series of the snow free date within the study area (Figure 1A) from 2000 to 2020 showed high interannual variation, but no significant trend over time (t₁₁₉₁ = 1.88, p = 0.059, 0.20 days/year; −0.002 to 0.419 95% CI) (Figure 1B). In 2019, the median snow free date was close to the long-term average, whereas the snow free date in 2020 was the earliest for over a decade (Figure 1B). On average and at the finer scale around the detected wheatear nests, the area was snow free more than 3 weeks earlier in 2020 (median = 130 doy; 40th to 60th quantile 125–143) compared to 2019 (155; 148–164) (Figure 2A). Similarly, vegetation growth reached the 80% threshold 17 days earlier in 2020 (162; 156–167) compared to 2019 (181, 178–184), and was strongly correlated with snow free dates (R² = 0.43, p < 0.001, t₂₅₄ = 7.68) (Table 1). Snow cover and vegetation height data from direct observations in the field matched remote sensed images (Figure 2A, date of snowmelt and 80%-NDVI amplitude).

In 2020, we retrieved 11 geolocators, of which nine contained sufficient data for further analyses (6 Migrate Technology, 3 SOI, of five males and four females). Fourteen tagged individuals were resighted (n = 35, return rate: 40%) of which 13 were recaptured. Two individuals lost the device. Nine adults of the control group were resighted (n = 16, return rate: 56%). The frequency of return rate between tagged and control groups was not significant (Fisher exact test, P = 0.126). Six tags recorded light levels for the full annual cycle. One tag stopped recording after arrival at the main wintering site (24EA). Two tags stopped recording during spring migration (24DK, 24EP). In total, we derived arrival dates at the breeding site from six tracked individuals in 2020 (Figure 2B and Table 1; see for details Supplementary Figures 1A–C and Supplementary Table 5).

Geolocation data reveals that in 2020 the arrival period at the breeding site for the six tagged individuals ranged from 13th April to 6th May 2020 (Figure 2B, Supplementary Table 5, and Supplementary Figure 1B). The three males (BT795, BT800, and BT801) were on average 13 days earlier than the three females (Figure 2B and Supplementary Table 6). Interannual variation in the timing of spring arrival/transition in the larger region of the study area (citizen science observations of wheatears) was high over the last 10 years (R² < 0.001), with an average arrival in 2019 (median doy = 119) and a rather late arrival in 2020 (median doy = 130) (Figure 2B, Table 1, and Supplementary Figure 2). For 2020, observation dates differed from the geolocator arrival dates by 12 days, with earlier arrivals of tagged birds (median doy = 118, Table 1).

Individual tracks revealed, compared to further stopovers in the same year, relatively long stopovers (approx. 20–30 days) before crossing the Sahara desert in autumn and shorter stopovers (approx. 5–10 days) just before reaching the wintering sites (Figure 3; see for details Supplementary Table 3 and Supplementary Figures 1A,C). Spring migration appeared to be faster, with fewer and shorter stopovers (approx. 5–10 days) (Figure 3, Supplementary Table 3, and Supplementary Figures 1A–C) until the coast of the Mediterranean Sea and Sardinia was reached. Here, most (4 out of 6) wheatears stopped again for a longer period (approx. 15–20 days, Figure 3, Supplementary Table 3, and Supplementary Figure 1A). With larger sample sizes, it will be possible to link data on migration routes, stopovers and arrival times more strongly to the breeding season, so we will have a much better idea of how non-breeding and breeding seasons are linked.

In 2019, we monitored 74 nests, of which 27 failed. In 2020, 53 nests were monitored, of which 23 failed. Nest failure was caused mainly by predation (23 nests in 2019, 15 in 2020), but also by nest abandonment or hatching failure (4 and 7, respectively). The cause of failure could not be identified in all cases (2 and 1, respectively). The first and last lay dates in 2019 were 28th May and July 18th, and in 2020 May 15th and July 14th. The number of monitored nests is assumed to be representative of this population, as they were located over an elevational gradient of 680 m which encompassed most of the elevation range of the species in the local area (attempts at higher or lower elevations are rare). According to results of the territory mapping (spanning an area of 1.15 km²), the study area holds 60 breeding pairs per km² and includes a range of different habitats (mountain grassland, grazed shrubland, and scree).

In 2020, for five nests out of all 53 attempts, we had enough information to distinguish replacements from first attempts (identified by colour ringed adults). For 2019, only one out of 75 nests could be identified as a replacement nest. For the lay date and the interval between arrival and breeding (arrival-breeding interval) estimates, we therefore excluded the replacement nests.

Lay dates in 2020 were 8 days earlier on average, although there was much variation, compared to 2019 (Figure 2B; median doy 2019: 158, range: 147–198, n = 73; median doy 2020: 150, range: 135–195; N = 48). However, due to later arrival in 2020 (citizen science observations), the arrival-breeding interval was shorter in the early year of 2020 (20 days) than in 2019 (39 days) (Figure 2B and Table 1). The median hatch date was 8 days earlier in 2020 (median doy: 167) than in 2019 (median doy: 175), but in 2020, we observed greater variability in hatch dates (range: 152–212) than in 2019 (range: 164–214). The difference in days between the peak of the green-up (NDVI 80% of amplitude) and hatching was similar in both years (difference in median dates: −6 days in 2019, +5 days in 2020), but in 2020, hatching was 5 days behind the peak (Table 1 and Figure 2B). Breeding phenology dates of tagged individuals did not differ from those of the studied population, but sample sizes were low due to nest failures in 2020 (Table 1).

The success rates (a, b) were lower in 2020 than in 2019. Even though the difference in hatching success between years was negligible, the probability that a hatchling would fledge differed significantly and was lower in 2020 (Figure 4; see for sample size and model output Supplementary Table 7). The nest success model (c) showed a non-significant lower success rate in 2020 (Figure 4 and Supplementary Table 7).