Long-term tracking of an Arctic-breeding seabird indicates high fidelity to pelagic wintering areas

Site fidelity is driven by predictable resource distributions in time and space. However, intrinsic factors related to an individual’s physiology and life-history traits can contribute to consistent foraging behaviour and movement patterns. Using 11 yr of continuous geolocation tracking data (fall 2008 to spring 2019), we investigated spatiotemporal consistency in nonbreeding movements in a pelagic seabird population of black-legged kittiwakes Rissa tridactyla breeding in the High Arctic (Svalbard). Our objective was to assess the relative importance of spatial versus temporal repeatability behind inter-annual movement consistency during winter. Most kitti wakes used pelagic regions of the western North Atlantic. Winter site fidelity was high both within and across individuals and at meso (100−1000 km) and macro scales (>1000 km). Spatial consistency in non-breeding movement was higher within than among individuals, suggesting that site fidelity might emerge from individuals’ memory to return to locations with predictable resource availability. Consistency was also stronger in space than in time, suggesting that it was driven by consistent resource pulses that may vary in time more so than in space. Nonetheless, some individuals displayed more flexibility by adopting a strategy of itinerancy during winter, and the causes of this flexibility are unclear. Specialization for key wintering areas can indicate vulnerability to environmental perturbations, with winter survival and carry-over effects arising from winter conditions as potential drivers of population dynamics.


INTRODUCTION
Site fidelity, the propensity of organisms to return to previously occupied locations, is a common form of behavioural consistency across taxa (Switzer 1993, Börger et al. 2008, Piper 2011. Such spatial consistency is closely related to resource quality and pre-dictability over space and time (Switzer 1993). Site fidelity for foraging patches may be favoured by individuals living in unpredictable environments, notably when resources are patchily distributed but spatially and temporally predictable (Switzer 1993). Returning to the same high-quality patch is reinforced by individuals' acquisition of local knowledge during previ-ous journeys, likely facilitating foraging efficiency (the 'always-stay' strategy; Switzer 1993, Irons 1998, Piper 2011. However, if changes in re source predictability and availability occur, fluctuations in behavioural consistency of individuals are expected to generate a decrease in site fidelity over time (the 'win-stay, lose-shift' strategy; Kamil 1983, Switzer 1993. The spatiotemporal predictability in re sources is also a strong driver of migratory movements during which individuals follow consistent seasonal changes in environmental conditions (Mueller & Fagan 2008). Examples of individuals tracking predictable resources along a migratory pathway, sometimes over remarkable distances, in clude whales following planktonic blooms or geese following the 'green wave' of plant phenological development in spring (Alerstam & Hedenström 1998, Kölzsch et al. 2015, Abrahms et al. 2019. Resource tracking across space and time is likely driven by a trade-off between (1) memory of location and timing of profitable patches and (2) exploration of novel environments when memory is unable to locate profitable patches, with repeatability governed by both extrinsic and intrinsic factors that generate or constrain profitability (Fagan et al. 2013).
In marine systems, top predators, such as seabirds, are dependent on patchily distributed resources (Weimerskirch 2007, Fauchald 2009). Because predictability in pelagic resources is habitat-and scaledependent, seabirds often rely on specific higherquality foraging zones (e.g. shelf edges, frontal zones, upwellings) that are predictable at meso (100− 1000 km) and macro scales (>1000 km ;Weimerskirch 2007). Specialization in foraging behaviour would be optimal when resources are predictable, thus stimulating site fidelity in high-quality patches that show consistency in productivity over time (Barraquand & Benhamou 2008, Carroll et al. 2018. Extrinsic factors related to interactions with other species or the environment (e.g. competition, prey availability, heterogeneity in resources) can therefore have profound influence in determining individuals' movements across landscapes and distribution (Fayet et al. 2017). However, oceans are dynamic environments and currently undergoing major changes that can affect re source predictability (Cury et al. 2008, Polovina et al. 2008, Hoegh-Guldberg & Bruno 2010. Temporal and spatial changes in marine resource predictability and availability can lead to shifts in the distribution of seabirds (Hamer et al. 2001, 2007, Ceia et al. 2014, Orben et al. 2015. The resilience and adaptability of populations to such changes are closely related to phenotypic plasticity and variability in behavioural traits of individuals, with populations showing high specialization in distribution, foraging behaviour, and diets being more sensitive to an alteration of their environment (Canale & Henry 2010, Patrick et al. 2015, de Grissac et al. 2016. Although seabirds often show within-individual wintering site fidelity (Ceia & Ramos 2015 and references therein), some species like the Cory's shearwater Calonectris borealis or the long-tailed skua Stercorarius longicaudus are highly flexible with some individuals shifting their winter distribution at the ocean scale (e.g. from the western to the eastern Atlantic Ocean; Dias et al. 2011, van Bemmelen et al. 2017. Flexibility in movement behaviour that is affected by extrinsic factors may only be measurable over large time scales, stressing the need for long-term movement tracking. Intrinsic factors related to physiological or life-history traits, such as sex, age, nutrition state, or breeding status and investment, can also contribute to individual variation in foraging and movement patterns (Phillips et al. 2017). Failed breeders often depart on migration earlier than successful breeders, and this difference sometimes persists into the overwintering period with the previous breeding status affecting the winter distribution and the timing of the arrival at and departure from the wintering site (Phillips et al. 2005, Bogdanova et al. 2017. Personality can also influence movement patterns, with bold individuals typically displaying higher foraging site fidelity, while shy individuals are shifting sites and travelling farther (Patrick & Weimerskirch 2014, Krüger et al. 2019, Harris et al. 2020. Moulting patterns, flight capabilities, and levels of stress hormones can affect the migration timing and an individual's ability to engage in long migratory flights (Dawson et al. 2000, Guglielmo et al. 2001, Schultner et al. 2014, Cherel et al. 2016. Moreover, important life history stages that compose the annual life cycle such as migration or reproduction are controlled by physiological processes whose timing or duration is re gulated internally by an individual's biological clock (the endogenous processes behind the biological rhythms) in response to temporally fixed cues (e.g. photoperiod; Kumar et al. 2010). Therefore, movement patterns driven by intrinsic factors with strong biological clock control (e.g. changes in physio logy, body condition, moult) will tend to occur consistently be tween years over an individual's life (Kumar et al. 2010, Wascher et al. 2018. Overall, the mechanisms underlying inter-annual movement consistency involve complex interplay between extrinsic and intrinsic factors, and the relative contributions of those factors is likely to vary over time (Phillips et al. 2017).
Using a long-term tracking dataset of 11 yr of continuous non-breeding movement data (fall 2008 to spring 2019) of black-legged kittiwakes Rissa tridactyla breeding in Svalbard (High Arctic Norway), we investigated inter-annual variations in space utilization and movement phenology during the nonbreeding period. The North Atlantic populations of kittiwakes largely congregate in winter at a shared staging area in the Western Atlantic Ocean (Frederik sen et al. 2012). Although the high degree of overlap in winter suggests site fidelity among and within populations, these results rely on 2 yr of data only (Frederiksen et al. 2012). Moreover, other studies showed that kittiwakes may alter their winter distribution according to intrinsic factors (e.g. stress hormone levels, breeding statuses, or reproductive investment; Schultner et al. 2014, Bogdanova et al. 2017, Whelan et al. 2020, which highlights the importance of multiple years of tracking data to assess the consistency of movement over time. We estimated inter-annual consistency in movement among and within individuals (2 to 7 yr repeated per individual) and how this consistency varied over the course of the non-breeding period across a decade. Our main objective was, therefore, to as sess the relative importance of spatial versus tem poral repeatability underlying inter-annual movement consistency. Site fidelity can emerge from resource specialization where resources are predictable so that individuals would show high distribution repeatability for consistent high-quality foraging areas (Switzer 1993). Higher repeatability for space utilization would then suggest that individuals are specialized on certain resources and are more vulnerable to extrinsic factors (e.g. environmental conditions, resource availability). Timing of climatic events (e.g. sea ice breakup, peak in primary productivity) can show high interannual variations associated with fluctuations in strength of large-scale climatic and oceanographic systems like the North Atlantic Oscillation (NAO; Visbeck et al. 2001) or ocean gyres (Polovina et al. 2008). Individual specialization driven by individuals following consistent patches learned on previous winters might thus be expected to show consistency in space over time. In contrast, higher repeatability in time than space would suggest that individual specialization is driven by factors occurring consistently at the same time each year in relation to the individual biological clock. The temporal and spatial consistency in movement patterns can thus inform about the relative importance of extrinsic versus intrinsic parameters behind the site fidelity observed in a given population, but intrinsic and extrinsic fac-tors often remain related to one another (i.e. hormonal clocks may change with foraging success driven by environmental variation), and the dominance of one factor over the other does not preclude the interaction or importance of both.

Logger deployment and geolocation processing
From 2008 to 2018, we deployed 276 geolocators (light-loggers or Global Location Sensors, GLS) on black-legged kittiwakes Rissa tridactyla to track their non-breeding movements (Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m676 p205 _ supp. pdf). Adults were captured using a ny lon loop attached to a fishing rod at the colony site in Kongsfjorden, Svalbard (High Arctic Norway; 78°54' N, 12°12' E) between May (pre-breeding) and July (chick-rearing) and equipped with geolocators. We used mk18 and mk13 (British Antarctic Survey), mk4083 and mk4093 (Biotrack) and Intigeo F100 and C65 (MigrateTechnology) mounted on a Darvic leg band. Devices measured light intensity every minute and recorded the maximum light intensity every 5 or 10 min. They also measured conductivity (as proxy for bird immediate environment, i.e. immersion) every 3 or 30 s and stored the number of wet measurements for every 10 min period. We recaptured 83% of the individuals at their return to the colony and recovered the geolocator (Table S1). Only complete annual tracks were used in the analyses after filtering out partial tracks caused by device failure or battery discharge. Overall, we acquired 200 complete tracks from 130 different individuals (see Fig. 1), covering 11 non-breeding seasons, continuously (fall 2008 to spring 2019). We tracked 33 individuals across multiple years (from 2 to 7 years) providing a total of 104 repeated tracks that were used to investigate within-individual consistency in movement. This repeated tracking was not always continuous, with gaps of 1 to 5 yr for 14 out of 33 individuals.
To infer geographic positions, geolocator data were processed according to the procedure developed for the SEATRACK project (Bråthen et al. 2021) and based on the threshold method calculating positions from twilight events ('coord' function from GeoLight package; Hill & Braun 2001, Lisovski et al. 2020. The procedure automatically identifies twilight events from raw light data ('twi-lightCalc' function from GeoLight package; Lisovski & Hahn 2012) and applies a set of filters to twilight events (removing or moving events from false day/ night detections or noise) and positions (speed, distribution limits, angle filter). Thus, all the geolocator data were processed automatically and consistently for all years of the study. Because light sensors from different geolocator models may differ, each track was calibrated individually. As such, the calibration method avoided systematic bias in latitude related to potential differences in light sensors among geolocator models or years of production. Based on the approach by Hanssen et al. (2016) and van Bemmelen et al. (2019), the calibration method used a set of criteria that allowed calibration of tracks from kittiwakes breeding in the Arctic (79° N), where constant daylight prevents calibration at the time of deployment and recapture. By plotting the latitude against time for a range of sun elevation angles and for each track ( Fig. S1 in the Supplement), the sun elevation angle that was manually selected (1) minimized the amplification of the latitudinal error close to the equinoxes, (2) resulted in matching latitudes at both sides of the equinox, (3) resulted in positions that fitted the latitude of the colony at the beginning and the end of the track and (4) fitted the shape and position of the oceans and continents when plotting the positions on a map (Fig. S2). The resulting sun elevation angle varied from −4.5 to −2.5° (mean angle −3.3°). The method also included rooftop calibration of geolocator models, with the purpose to select model specific thresholds that would result in ap proximately the same sun elevation angles among geo locator models. The mk-series geolocators from the British Antarctic Survey and Biotrack were as signed a threshold of 1 unit, while Intigeo geolocators from Migrate Technology were assigned a threshold of 11 units.
Although longitudes can still be determined reliably around the equinoxes, estimation of latitudes is inherently imprecise during this period, because day length is similar around the globe . Therefore, locations around equinoxes were excluded (8 Sep−20 Oct, 20 Feb−3 Apr; Bråthen et al. 2021). Additionally, continuous daylight during the polar summer (or continuous night during polar winter) does not allow geolocation-based tracking using light-level sensors. To fill these gaps and reduce biases along the trajectories, missing locations were re-estimated by interpolation between known locations using an algorithm that was specifically developed for SEATRACK (Fauchald et al. 2019), based on a method originally proposed by Technitis et al. (2015). In short, this algorithm is based on the determination of so-called space-time prisms, which are 3-dimensional volumes defined by the coordinates (x,y) and time (z). The space-time prism delineates all the potential paths that can be followed by an individual moving from point A to point B, given 3 parameters: the distance from A to B, the time budget available, and the maximum rate of movement (Miller 1991). When projected onto a 2-dimensional plane, the space-time prism becomes the potential point area (hereafter Ppa; Technitis et al. 2015). Although the 3dimensional representation of the space-time prism is useful to understand its concept (Neutens et al. 2007), it is naturally more convenient to work with only 2 dimensions when dealing with discrete time steps, as is the case in tracking studies, where locations are obtained at specific time intervals. Computing the Ppa in this context is straightforward (Technitis et al. 2015), given that the 3 above-mentioned parameters are known. Let us consider a startpoint (A) and start time (t i−1 ), and an endpoint (B) and end time (t i+1 ). Knowing the maximum rate of movement and the time t i at which a new location (N i ) is to be created, one can determine the circle defining the maximum range (rg i−1 ) from point A to the new location and the circle defining the maximum range (rg i+1 ) to point B, centered on B. The Ppa corresponds to the area of overlap between those 2 circles of maximum range, i.e. the area delimiting all locations that are reachable from both A and B, given the time budget and maximum movement rate. This process can be repeated any number of times, depending on the number of new locations that need to be generated. The new locations are generated in a random order (i.e. not chronological), thus creating a sort of correlated random walk respecting the constraints set by the relative position of A and B, the time budget, and the maximum movement rate. Here, we used a dynamic value for the maximum movement rate parameter, based on the distribution of observed movement rates as a function of time elapsed between 2 locations from the dataset. To do so we calculated, based on each individual track, the movement rates for random combinations of known locations separated by varying time-intervals. We used the 75 th percentile from that distribution as the maximum movement rate (Fig. S3). The 75 th percentile was computed by quantile regression, using the function 'rq' from package quantreg (Koenker 2020). Finally, the algorithm uses additional information to constrain the new positions obtained: (1) im mersion data to determine attendance at the colony and force a new location to remain close to the colony during the breeding season, (2) land masks (land filters) to constrain positions over the ocean, (3) longitudes (ob tained from the geolocator data, as longi tude can still be estimated during the equinoxes), and (4) light levels to determine whether the new position was north of the latitudinal limit of the polar day in summer or night in winter (i.e. continuous day/night recorded by the loggers).
In all further analyses, an annual track refers to the non-breeding period extending from the colony departure in fall to the return to the colony area the following spring. Departure from the colony and return to the colony were identified using Lavielle partitioning algorithm ('ts.LaviellePart' function from R package adehabitatLT v.0.3.25; Calenge 2006, Barra quand & Benhamou 2008) over a 5 d running maximum of the saltwater immersion data indicating a transition between land use (mostly dry) and continuous pelagic behaviour (mostly wet). Departure and arrival dates were adjusted according to visual inspection of the individual's locations right after the behavioural transition from land use to pelagic in fall, and right before the transition from pelagic to land use in spring. In spring, foraging trips after the first visit to the colony area were excluded, as individuals start to display a central place foraging behaviour, including long pre-laying trips as far as Iceland (Bogdanova et al. 2011).

Consistency in intertrack distances
To estimate the consistency in non-breeding movement over the entire annual tracks either among or within individuals, we used an approach based on the nearest neighbour analysis (for similar methods see Guilford et al. 2011, van Bemmelen et al. 2017. For each location of a focal track, we calculated the orthodromic distance to the nearest location on (1) a randomly selected track from an other individual to estimate among-individual consistency and (2) a track from the same individual but from another year to estimate within-individual consistency. Prebreeding movements were ex cluded, and positions were considered to be fixed at the colony after the first visit in the colony area in spring. This nearest neighbour distance was calculated over a large time window (60 d) to assess spatial consistency in movement. The 60 d time window was se lected after running a sensitivity analysis using different time windows varying from 1 to 120 d (with 10 d intervals) to assess when apparent variation related to timing differences fades (Fig. S4). The time window we selected allows spatial comparison without overlaps between fall and spring migrations that would be created if using a larger time window (i.e. >100 d). We repeated the method described above with a 1 d time window (i.e. daily comparison) to compare dissimilarities among tracks generated by both spatial and timing differences. We bootstrapped the resulting distances 10 000 times to calculate the mean intertrack distance per day among and within individuals. The results were log-transformed before further analysis to meet the assumptions of homo scedasticity and normality of residuals. We fitted linear mixed-effects models (LMER, R package lme4 v.1.1-23; Bates et al. 2020) to determine if among-and within-individual distances differ, using the mean intertrack distance for each annual track, with individual and year as random factors. Sex had no effects on individuals with known sex (LMER; using the 60 d time window: β = 0.05, SE = 0.05, df = 107.1, p = 0.318; using the 1 d time-window: β = 0.04, SE = 0.03, df = 116.7, p = 0.191) and was therefore discarded from the final models on all individuals. Following the method used by van Bemmelen et al. (2017), the mean intertrack distance calculated over the large time window (60 d) was mapped to illustrate both among-and within-individual spatial consistency in movement during the non-breeding period.

Variability in migration schedule
To illustrate the variability in migration schedule, we extracted parameters associated to important migration phenology events: (1) the departure from the colony, (2) the start of the southward migration movement, (3) the crossing of the Arctic Circle in fall corresponding to the end of the post-breeding staging in the northern regions (Barents and Greenland Seas), (4) the start of the northward migration movement and (5) the arrival at the colony defined as the first visit to the colony area in spring. We estimated the variability in the timing of each phenological parameter at the population level. To do so, we calculated the difference between the latest date in a given year and the dates from all other individuals in the same year and then averaged these differences over the 11 yr of the study. We also estimated the within-individual variability in migration timing for each phenological parameter for individuals tracked multiple years (n = 33 individuals, 104 annual tracks). This was done by calculating the difference between the latest date for a given individual and the dates from all the other years this individual was tracked and then averaged over all individuals. Mean dates are reported with their standard deviation. Finally, we estimated the individual repeatability r (intra-class correlation) for the timing of each phenological parameter. This was done using the function 'rpt' of the R package rptR v.0.9.22 (Stoffel et al. 2017) and using only individuals with multiple years of tracking (see Table S3). We used 'year' as random effect only for the 'colony departure' to account for important interannual variation in this parameter. The repeatability estimate r is not a measure of absolute consistency but the proportion of total variance accounted for by differences between groups (i.e. in our case among individuals; Nakagawa & Schielzeth 2010). All analyses were carried out in R version 4.0.2 (R Core Team 2020).

Non-breeding staging areas
We used Hidden Markov Models (HMMs; Zucchini et al. 2017) to examine at-sea behaviour of individuals and identify the sequence of discrete behavioural states that best fitted the non-breeding tracks of individuals. HMMs were fitted to all individuals at once using R package moveHMM v.1.7 (Michelot et al. 2016) with gamma and von Mises distributions to describe the frequency of step length and turning angle distributions, respectively. Different initial parameter values were tested to ensure numerical maximization of the likelihood through the iteration process. A 4state model was selected after examination of the pseudo-residuals (Michelot et al. 2016) and because it better fitted the geolocation data than simpler 2-or 3-state models based on AICs and initial inspection of distribution of movement parameters (see Table S4 and Fig. S5). States 1 and 2 were de fined by short steps (< 85 km) with either frequent shifts in direction (angle concentration of 0; State 1) or moderately directional movement (angle concentration of 1; State 2) and attributed to periods of staging in more intensively utilized areas (see Table S4). States 3 and 4 were defined by long steps (>150 km) with moderate shifts in direction (angle concentration of 0.5; State 3) or highly directional movement (angle concentration of 7; State 4) that characterized transient and commuting behaviours during travelling periods (see Table S4). In further analyses, the first 2 states were combined into a broader 'staging state' and the last 2 states into a 'travelling state' to de termine stationary and travelling positions, respectively. We were only interested in the stationary positions (staging state), which we mapped separately for the fall migration, winter, and spring migration periods to illustrate the main staging areas during each period (see Fig. 1). Therefore, the 'staging areas' re fer to more intensively utilized areas throughout the non-breeding period. Finally, we calculated the 80 and 50% utilization distribution kernels (UDs) over the stationary positions projected using a Lambert Azimuthal Equal Area coordinate system and the R package adehabitatHR v.0.4.18 (Calenge 2006) with a smoothing factor (h) of 200 km and grid cells of 50 × 50 km. We used the kernels to identify important staging areas for our study population and to illustrate the correspondence between these staging areas and the spatial consistency in the intertrack distance of individuals.

Individual tracking and migratory routes
The mean annual departure date from the colony varied from 27 August to 25 September (5 September ± 9 d [SD]). Individuals staged within the Arctic in the Barents and Greenland Seas (mean annual staging from 38 to 64 d, 52 ± 7 d) before migrating southwest along a corridor between East Greenland and Iceland (Fig. 1). The mean annual date of onset of fall migration ranged from 10 to 29 October (19 October ± 5 d). All individuals spent the winter in the North Atlantic Ocean, with the main wintering area extending from the Grand Banks of Newfoundland to the mid-Atlantic ridge. The winter distribution of the population was therefore largely pelagic, but alternative staging areas were also used along the continental shelves of Northeast America and Western Europe ( Fig. 1 and see UDs of Fig. 3). The northeast movement of spring migration (mean annual starting date from 31 March to 9 April, 4 April ± 3 d) was spread over a larger front than in fall, with routes passing both north and south of Iceland (Fig. 1). The mean annual date of arrival in the colony area varied from 10 to 20 April (15 April ± 5 d).

Non-breeding movement consistency in space
The nearest neighbour analysis conducted over a large time window (60 d) indicated high spatial consistency during the non-breeding period in the studied population, both within and among individuals ( Fig. 2A, see also examples of repeated tracks from several individuals in Fig. S6). The mean within-individual intertrack distance was 261 km (95% CI: 180−376 km), significantly lower than the mean intertrack distance of 545 km found among individuals (95% CI: 377− 793 km; LMER, β = 0.71, SE = 0.04, 95% CI: 0.62; 0.80, df = 218.9, p < 0.001). The mean within-individual intertrack distance was consistent regardless of the number of years individuals were tracked (LM, β = 12.2, SE = 23.7, t 31 = 0.5, p > 0.609), indicating that site fidelity persists over longer tracking periods (range 2 to 7 yr). Mapping intertrack distances (Fig. 3) showed areas of high spatial consistency among individuals in the western part of the North Atlantic Ocean as well as in the Greenland Sea (east coast of Greenland and Svalbard) and in the Barents Sea (between Svalbard and Novaya Zemlya, Fig. 3A). These sectors of high spatial consistency correspond to the main staging areas identified with the Hidden Markov Models (Fig. 1) and illustrated with 80 and 50% utilization distribution kernels (Fig. 3A). Similarly, areas associated with high consistency within individuals were also associated with staging areas in the Northwest Atlantic Ocean and the Greenland and Barents Seas but also with areas along the coasts of Northeast America and, to a lesser extent, Western Europe (Iberian Peninsula and British Isles, Fig. 3B). Across the non-breeding season, the kittiwake distribution was mostly pelagic with 73% of the locations (32 327 out of 44 071 locations in total) in areas deeper than 500 m. High spatial consistency within individuals was also found in coastal areas and shallow waters of the Barents Sea (see Fig. 3), but the use of these areas was limited. The utilization of deeper waters was even more evident during the winter period (mid-November to mid-March), with only 12% of the locations found in areas of 0−500 m depth (2816 out of 24 034 winter locations). Overall, spatial consistency within individuals remained high over deeper waters (Fig. S7), indicating that site fidelity was also common in pelagic areas.

Non-breeding movement consistency in time
The nearest neighbour distance analysis conducted over a short time window (daily comparison to include variations associated with timing effects) also showed higher consistency in mean intertrack distances ( Fig. 2B; LMER, β = 0.39, SE = 0.02, 95% CI: 0.35; 0.44, df = 204.9, p < 0.001) within (847 km, 95% CI: 694−1026) than among individuals (1257 km, 95% CI: 1034−1546). Using a short time window in the analysis revealed peaks in timing variability during fall (October to mid-November) and spring (mid-March to mid-April) associated with migratory stages (Fig. 2B). Individual repeatability estimation of pheno logical parameters during migration stages (Fig. 4B) showed low individual repeatability in timing of post-breeding colony departure (r = 0.20, 95% CI: 0.02; 0.39, p < 0.001) driven by a relatively high intra-individual variability (see Fig. 4A and Table S4), high repeatability in the onset of the fall migration movement (r = 0.80, 95% CI: 0.65; 0.88, p < 0.001), and moderate repeatability in the timing of the crossing of the Arctic Circle after the post-breeding staging in the Barents and Greenland Seas (r = 0.54, 95% CI: 0.33; 0.69, p < 0.001). In spring, the very low repeatability estimates in the timing of the onset of the migration movement (r = 0.36, 95% CI: 0.12; 0.55, p = 0.002) and the arrival in the colony area (r = 0.28, 95% CI: 0.04; 0.49, p = 0.016) were driven by high consistency in the timing of these events among individuals (see Fig. 4B and Table S3). Overall, the timing of phenological events, both within and among individuals, showed higher variability during the fall period than during the spring period ( Fig. 4A and Table S2), indicating the spring migration occurred during a shorter time interval.

DISCUSSION
Kittiwakes showed high spatial consistency in nonbreeding movements at both macro (>1000 km) and meso scales (100−1000 km), suggesting that individuals benefit from predictability in resources at these coarse scales. The significantly higher consistency within individuals than among individuals indicated specialization in space-use strategy at the individual level. With a deviation in routes among years of less than 300 km on average, the inter-annual spatial consistency of individuals indicated important site fidelity at the individual level, especially when considering the coarse resolution of geolocator measurements and the wide ocean-scale distribution and long migrations of the studied population. In contrast, we found low repeatability in the timing of important phenological events in general, which suggests that those intrinsic timing mechanisms known to be repeatable across years, such as the individual's biological clock, may not have been a strong driver of the movement consistency of individuals over the years of the study (but see Section 4.3). Thus, repeatability is likely governed by memory of the location of high-quality patches relatively unconstrained by intrinsic factors. This individual spatial consistency was stable over the entire non-breeding sea-

Fidelity to deep-water areas
High fidelity to staging areas outside the breeding season is not uncommon in seabirds and is often associated with more predictable areas along continental shelves and less with oceanic habitats (Weimers kirch 2007). For instance, highly productive up welling habitats along shelf edges in the Canary and Benguela current systems were associated with site fidelity during the non-breeding period in black-browed alba trosses Thalassarche melanophrys, Cory's shearwaters Calonectris borealis and long-tailed skuas Ster co rarius longicaudus (Phillips et al. 2005, Dias et al. 2011, van Bemmelen et al. 2017. Similarly, kittiwakes tracked in this study used the shelf edge of the Grand Banks, an area known for its high biological productivity (Heywood et al. 1994, Maillet et al. 2005, Frederiksen et al. 2012. However, many individuals also wintered in deep waters, from the edge of the Grand Banks plateau to the mid-Atlantic Ridge. Despite this highly pelagic distribution, individuals showed area fidelity, suggesting this deep-water area can provide habitats with enough coarse-scale predictability in resources to stimulate area fidelity. In comparison, black-legged kittiwakes tracked in the North Pacific also showed some degree of individual fidelity to pelagic areas (~25% of locations within 400 km grid squares; Orben et al. 2015), but to a lesser extent than what we found in the North Atlantic (~75% of locations within 400 km of mean nearest neighbour distance). However, a similar degree of spatial consistency (median nearest neighbour distance 250 to 400 km) to what we found with kittiwakes was described in long-tailed skuas for the same oceanic area west of the mid-Atlantic ridge, which they used as a staging area during migration (van Bemmelen et al. 2017).
This oceanic region is crossed by the subpolar front extending from the Newfoundland Rise to the mid-Atlantic Ridge and characterized by a strong horizontal temperature gradient, eddies, and nutrient mixing and retention inducing biological productivity enhancement (Heywood et al. 1994, Scales et al. 2014, Hátún et al. 2016. There is evidence that this mid-ocean frontal zone is an important staging site and a diversity hotspot for multiple seabirds (e.g. Guilford et al. 2009, Egevang et al. 2010, Weimerskirch et al. 2015, as well as other marine predators such as sharks (Queiroz et al. 2016), tunas (Walli et al. 2009), chelonioid turtles (Eckert 2006), and cetaceans (Doksaeter et al. 2008, Skov et al. 2008. This large overlap in distribution has stressed the potential vulnerability of marine vertebrate populations relying on this area to large-scale changes in environmental conditions affecting resource predictability and availability (Frederiksen et al. 2012), such as the weakening of the subpolar gyre and warming of the North Atlantic (Häkkinen & Rhines 2004, Descamps et al. 2013, Fluhr et al. 2017, Hátún et al. 2017. Populations showing important site fidelity are expected to be more sensitive to extrinsic factors affecting resource predictability and might thus be particularly impacted by such carry-over effects (Phillips et al. 2017).

Movement strategy in space
To a certain extent, marine predators are expected to adjust their space-use strategy in response to environmental changes affecting resource predictability of foraging patches over time (Davoren et al. 2003, Wakefield et al. 2015. This ability to respond to environmental variability by using a 'win-stay, lose-shift' strategy directly depends on the flexibility and plasticity in behavioural traits intrinsic to the population and individuals (Canale & Henry 2010). Area shifting was uncommon in our study, and we found an overall high spatial consistency, even for individuals tracked up to 7 years. This long-term area fidelity could be expected to arise if individuals benefited from consistency in resource availability in known staging areas over the course of the study, thus preventing the need for important shifts in distribution at meso or macro scale. Alternatively, long-term site fidelity can be reinforced by site familiarity and the benefits of acquiring information specific to an area, leading to individuals favouring an 'always-stay' strategy (Irons 1998, Wakefield et al. 2015. Advantages of returning to a known area include increased foraging efficiency through knowledge about food location and availability, movement efficiency by using prevailing wind corridors, dominance during competitive interactions as well as avoiding potential risks of visiting unfamiliar places, such as higher risks of getting stranded (Piper 2011). In contrast, some seabird populations show high flexibility in their non-breeding distribution with individuals possibly having several preferred migratory strategies (Dias et al. 2011, van Bemmelen et al. 2017. For example, some long-tailed skuas shifted their winter distribution between years at the ocean scale, between the Benguela current and the Falkland current (van Bemmelen et al. 2017). They followed a specific route for each alternative wintering site and kept the same route over the years, indicating this shifting behaviour was not accidental but based on past experience. Similarly, Cory's shearwaters used wintering areas independently of prevailing wind currents encountered on route, suggesting the choice of using one or another alternative site was deliberate and predetermined (Dell'Ariccia et al. 2018). We also found that some kittiwakes showed flexibility in their migratory decisions by shifting their staging areas or by displaying a more exploratory behaviour with different degrees of itinerancy, although this occurred at a more modest scale than in long-tailed skuas and Cory's shearwaters. If a decrease in resource availability can stimulate area-shifting, the cause of the flexibility ob served in some individuals can emerge from a diversity of factors both extrinsic or intrinsic, and thus remains unclear. For instance, the breeding status is likely to vary across years and is known to introduce inter-annual variability in individuals' space-use strategies (Phillips et al. 2017). In some species, failed breeders engaged in longer (Bogdanova et al. 2011) or shorter migration (Phillips et al. 2005) compared to successful breeders, indicating that the choice of wintering sites can be conditiondependent to the breeding investment.

Plasticity in timing of movement
The dynamics of large-scale climatic and oceanographic systems, such as the subpolar gyre, generate strong interannual variations (up to several weeks) in the timing of biological productivity (i.e. phytoplankton bloom) in the North Atlantic (Gaard et al. 1998, Henson et al. 2009), with cascading effects on higher trophic levels (Henson et al. 2009, Eliasen et al. 2011. In response to these interannual fluctuations, marine predators following foraging patches learned on past travels are expected to show higher movement consistency in space than in time. This pattern is what we observed with individuals showing high spatial consistency, suggesting that individual specialization is driven more by space utilization in relation to extrinsic factors with higher predictability in space than in time. Similarly, when investigating phenological parameters associated with migratory movement, we observed overall low repeatability in timing across years. Notably, we found low repeatability in the timing of the colony departure, a phenological parameter often related to breeding status and reproductive investment (Bogdanova et al. 2011, Whelan et al. 2020. In another study, food supplementation enabled fed kittiwakes to initiate departure from the colony earlier than unfed kittiwakes, indicating that individuals experiencing lower breeding costs are able to transition into the non-breeding season in better condition and initiate migration earlier (Whelan et al. 2020). Interestingly, high repeatability only occurred for the timing of the onset of the fall migration. Considering the large range of departure dates in the population, this high within-individual repeatability suggests that intrinsic timing mechanisms played an important role in the phenology of the fall migration movement. This has been demonstrated with blue whales Balaenoptera musculus, where movement phenology was driven not by proximate environmental factors, but rather by the average timing of these factors across years at specific staging sites (Abrahms et al. 2019). This highlights the importance of memory and individual biological clock as intrinsic factors regulating the timing of individual movement in this long-lived species.
Additionally, the low individual repeatability in phenology of the spring migration was driven by low variation among individuals. This synchronicity in phenology in spring, coupled with high variation in timing of phenological parameters among individuals in fall, indicated that spring migration might be more time constrained. This is often the case for polar and subpolar migrants, as the optimal arrival date at the breeding site is a compromise between costs and benefits of an early arrival with individuals benefiting from e.g. a better nesting site or higher mating success while confronted with the risk of being exposed to harsh environmental conditions in early spring (Kokko 1999). These constraints are usually relaxed or absent in fall, leading to more variability in the timing of the fall migration (Nilsson et al. 2013), but carry-over effects arising from the breeding status and investment can also contribute to among-individual variations in post-breeding migration patterns (Harrison et al. 2011).
In conclusion, black-legged kittiwakes were remarkably consistent in their overwinter locations among years despite wintering primarily in pelagic regions of the North Atlantic. The consistency was stronger within than among individuals, implying that individuals were using memory to return to profitable locations year upon year. Consistency was also stronger in space than in time, suggesting that it was driven by consistent resource pulses that vary in time more so than space, and that intrinsic drivers (photoperiod, hormones, condition) that are known to be repeatable within individuals were somewhat less important. Of course, intrinsic and extrinsic factors often interact with one another (i.e. foraging success can influence body condition which can influence hormones), and in reality, both types of factors will be integrated by the individual to decide when to stay or go. Consistent use of key foraging locations is likely associated with foraging success, risk of starvation, and a strong driver of overwintering survival, linking variation in timing of resources to individual fitness and, ultimately, population trends. We are thankful to the Norwegian Institute for Nature Research, the Norwegian Polar Institute as well as the French Polar Institute and Alfred Wegener Institute for Polar and Marine Research (AWIPEV) for their logistical support in the field and to the Centre de la Science de la Biodiversité du Québec for their help with travel logistics. We are very grateful to Antonio A. Cuba Martinez for his precious help with R coding and to the Krykkjefjellet field teams without whom this study would not have been possible, including Alexandre Corbeau, Céline Clément Chastel, Hilde Dørum and Anna Lippold. This study was approved by the Norwegian Food Safety Authority (FOTS ID 2086, 3319, 4169, 6291, 6348, 8482, 15603, 15611, 19970), the Governor of Svalbard and the McGill University Animal Care Committee.