Wintering North Pacific black-legged kittiwakes balance spatial flexibility and consistency

Background Marine environments are inherently dynamic, yet marine predators are often long-lived and employ strategies where consistency, individual specialization, routine migrations, and spatial memory are key components to their foraging and life-history strategies. Intrinsic determinates of animal movements are linked to physiological and life-history traits (e.g. sex, colony, experience), while extrinsic influences occur as the result of an animal’s interactions with either other animals or the environment (e.g. prey availability, weather, competition). Knowledge of the factors affecting animal movements is critical to understand energetic bottlenecks and population dynamics. Here, we attempt to understand the interaction of some of these factors on the winter distributions of a surface-feeding seabird in the North Pacific. Between 2008 and 2011, we tracked 99 black-legged kittiwakes breeding at St. Paul and St. George in the Pribilof Islands, Alaska using geolocation loggers. We tested for colony and sex differences in winter distributions, and individual spatial fidelity over two consecutive winters of 17 individuals. Then we linked tracking data to associated environmental conditions as proxies of prey availability (e.g. sea surface temperature, mesoscale eddies, chlorophyll a, and wind) to understand their influence on kittiwake space use at an ocean basin scale. Results Black-legged kittiwakes from both Pribilof Islands primarily wintered in pelagic sub-arctic waters, however, distributions spanned seven ecoregions of the North Pacific. There was a high degree of similarity in area use of birds from the two closely situated colonies and between sexes. Birds tracked for two consecutive years showed higher fidelity to wintering areas than occurred at random. Annual changes were apparent, as distributions were further north in 2009/10 than 2008/09 or 2010/11. This occurred because 70 % of birds remained in the Bering Sea in the fall of 2009, which corresponded with lower October sea surface temperatures than the other two years. Conclusions Although individuals returned to wintering areas in consecutive years, our results suggest that under current conditions individual black-legged kittiwakes have a high capacity to alter winter distributions.


Background
Marine environments are inherently dynamic, yet marine predators are often long-lived and employ strategies where consistency, individual specialization, routine migrations, and spatial memory are key components to their foraging and life-history strategies [1][2][3]. At small scales, prey species such as squid, fish, and krill are patchily distributed-both vertically and horizontally in the water column [4,5]. Surface foraging seabirds are adapted to find food at small scales through social communication, olfactory and visual cues [6,7]. During the non-breeding period many of these species travel thousands of kilometers to upwelling regions and frontal zones where prey are more predictable [8,9]. This results in individuals choosing migratory destinations at an ocean basin scale, presumably without knowledge of conditions at their destinations [10].
Where and how animals move across a landscape is driven by a myriad of factors than can be simplified into intrinsic or extrinsic factors [11,12]. Intrinsic determinates of animal movements are linked to physiological and life-history traits of species (e.g. reproductive status, past experience, physiological needs and capacity, navigation abilities), while extrinsic influences occur as the result of an animal's interactions with either other animals or the environment (e.g. competition, predation, prey availability). While disentangling intrinsic and extrinsic factors and understanding how they may interact is difficult [13], intrinsic factors appear to be particularly important. For example, for seabirds, age, sex, and breeding colony may influence migrations [14][15][16][17][18]. As a consequence, these differences imply individuals within populations will be differentially affected by extrinsic factors such as prey availability. Furthermore, there is undoubtedly an influence of past experience and learning that shapes some portion of migratory paths, wintering areas, and prey preferences in long-lived seabirds [19][20][21][22][23][24]. Yet, over a lifetime individuals are likely to be flexible in their foraging choices as oceanic habitats are dynamic and conditions change from one season to the next [25,26]. Understanding how intrinsic and extrinsic factors interact to influence animal movements is paramount for an understanding of population trends and spatial capacity, particularly in highly mobile long-lived species.
Black-legged kittiwakes show spatial and dietary variability in foraging habitats both within and between colonies [27,28], and throughout the annual cycle, including the winter non-breeding period [29][30][31][32][33]. Breeding individuals are known to display marked fidelity to locations and foraging in concert with tidal cycles [34], yet little is known about the importance of fidelity during the non-breeding period (but see [35]). In this study, we use the black-legged kittiwake as a model species to understand the importance of flexibility and consistency during their winter migrations to the dynamic pelagic environment. First we address the effects of breeding colony location, sex, and individual experience on the winter distributions of the black-legged kittiwake breeding at the Pribilof Islands. Then, we link three winters of kittiwake tracking data to associated environmental conditions and oceanic habitats to understand how annual habitat conditions influence space use at an ocean basin scale.

Methods
To study black-legged kittiwake (hereafter kittiwake) wintering ecology, geolocation loggers (1. 8 Table 1). All birds were captured off active nests using either a noose pole or foot snare supplemented with a CO 2 powered net gun containing custom nets (Super Talon Animal Catcher, Advanced Weapons Technology, California) at recapture. At deployment and recapture birds were weighed and measured, nest contents were recorded and blood samples for sexing via DNA were taken [36]. Overall, 86.6 % of birds were resighted and 77 % of loggers were recovered (Table 1). Seven loggers failed during data recovery and 13 loggers failed before recapture. Over the three study winters, 17 birds carried loggers for 2 winters. A total of 113 complete winter trips were recorded from 99 kittiwakes. It is possible that though the loggers used in this study were small (~0.6 % body mass or less) they still may have had deleterious effects on the individuals carrying them [37,38], however no negative effects have been thus far reported for kittiwakes carrying geolocation loggers in terms of body mass, breeding participation or survival [29,31,32].
Data processing, spatial analyses and statistical tests were conducted using MATLAB (v2014a, The Mathworks, Natick, MA) and R v3.1.1 (R Development Core Team, 2014). Significance was set to p < 0.05. Birds were not individually marked in 2008 so resights rely on birds attending the same nests in subsequent years. Like recaptures, resights presented are cumulative, so birds from 2008 had three seasons of following effort, while birds from 2010 only had one year of recapture and resighting effort. Our objective was to recatch rather than resight birds, which often required hiding from birds rather than reading alphanumeric bands. Most birds were deployed with GPS and GLS loggers for a summer foraging study [27,28], were subsequently recaptured and GLS loggers redeployed for overwinter. Thus deployment numbers also include 7 birds that were not recaptured during the summer breeding season (GPS tags fell off when tail feathers were molted); subsequently only 2 of these birds were recaught or resighted. Omitting these birds (and those whose nesting ledges fell) gives an overall resight of 91 % a Includes 1 bird deployed for the summer with a GPS and a GLS logger; the GPS was not recovered and the bird was not resighted in subsequent seasons b 3 birds were deployed in on a cliff face that collapsed overwinter and never resighted. These birds may have relocated to other areas of the colony where, due to the size of the cliffs, we were unable to resight them c Includes 2 birds deployed for the summer with a GPS and a GLS logger; the GPS was not recovered and birds were not resighted in subsequent seasons d Includes 2 birds deployed with both a GPS and GLS logger; the GPS was not recovered. 1 bird was recaught with a GLS logger in 2009, the other was not resighted e Includes 2 birds deployed with both a GPS and GLS logger; the GPS was not recovered. 1 bird was recaught with a GLS logger in 2010, the other was not resighted f High Bluffs, where 2 GLS loggers were deployed in 2010, was only visited twice in 2011 and 1 bird was never seen

Geolocation processing
Geolocations were estimated from September 1 thru May 30 with the colony as a fixed start and end location using the 'TripEstimation' package in R to implement a Bayesian model that incorporates a land mask and flight speed into the prior distribution and calculates locations from light levels using the 'template-fit' method [39][40][41].
We used a mean speed of 33 km hr −1 [42], however tag derived activity data (time spent dry, [32]) indicated that birds only spent 15 ± 9 % of each 24 hr period flying, therefore, as a conservative estimate we chose to assume that birds spent 33 % (two standard deviations from the mean) of their time in flight, restricting an individual's range over 12 hours, on average, to 130km. Thus a mean speed of 10.89 km hr −1 and speed variance equal to half this were entered into the model to parameterize a lognormal distribution. The most probable track was then obtained through a Kalman filter with six Markov Chain Monte Carlo (MCMC) simulations of 1000 iterations each, after a 500-iteration burn in period. Tracks were visually compared and found to be similar to those calculated using smoothing and filtering methods that result in errors of 186 ± 114 km [43]. Subsequently, all analysis were constrained to October thru February to avoid greater errors inherent in locations estimated around the equinox period.

Distributions
How birds share space over large scales relates to the size of individual ranges, the distribution of these ranges, and the density of birds. Individual range size was calculated as the number of 45 × 45 km grid squares occupied by the last chain of the MCMC estimation (1000 iterative tracks). To assess how variable each group (colony, sex, year) was in the areas they used as well as a quantitative assessment of sample sizes, the cumulative number of grid cells occupied was calculated with the addition of one track from each individual randomly selected for 10,000 iterations [44]. Groups with less variability or higher amounts of area shared between individuals, have shallower curves than groups with high individual variability in area use [44]. To assess the amount of shared area, an index of similarity was calculated as the ratio of the number of shared grid cells (45 km) occupied by the last chain of the MCMC estimation to the total grid cells used by both groups, where identical groups would be equal to 1. To test these indexes against what similarity might randomly occur, birds were randomly assigned to groups and the index of similarity was calculated 10,000 times with the same number of individuals in each group as the original dataset. P-values were calculated as the percentage of iterations that resulted in an similarity index smaller than observed [32,45]. Finally, to assess relative density of birds and overall distributions, densities were calculated from the last chain of the MCMC estimation for each individual over a 5 × 5 km grid (to aid in visualization). This method incorporates the uncertainty in the location estimates and negates the reasons to use a method such as a kernel density estimate [40,46].

Individual Spatial Fidelity
Area fidelity can occur at varying spatial scales. To identify the spatial scales in which kittiwakes showed fidelity to wintering areas we calculated a monthly index of similarity, a measure of the number of shared grid cells, between all estimated locations for repeat trips from individual birds (n = 17) at a series of grid sizes ranging from 10-400 km. Then to quantify a random measure of overlap we calculated the percentage of overlap between 59 random pairs, between individuals of the same sex and colony from consecutive years.

Annual difference in large marine ecosystem use
To better understand annual changes in large-scale distributions and habitat use of kittiwakes we assessed the use of oceanic biogeographical provinces of the North Pacific [9]. For each province, at the scale of 45 km grid cells, we calculated the percentage of the biogeographic province occupied by kittiwakes (again incorporating location error by using the last chain of the MCMC estimates) (# grid cells occupied/total grid cells in province), average density of occupied grid cells, and the proportion of the overall bird distribution in each province (# grid cells occupied inside region/total # grid cells occupied). We also assessed how individuals used these regions by calculating the number of ecoregions used by individual birds and by sexes.

Annual habitat conditions
Habitat variables were extracted along the best-fit tracks at a 1°grid scale. Sea surface temperatures (SST) were extracted as an eight-day blended product from data hosted by NOAA's Environmental Research Division (http://oceanwatch.pfeg.noaa.gov/thredds/catalog.html). Bathymetry was extracted from the ETOPO1 dataset [47], and bathymetric slope calculated. As eddies are known to condense prey for surface foraging predators including kittiwakes [28,48], we extracted sea surface height (SSH) and surface currents used to calculate eddy kinetic energy (EKE) were extracted from the Navy Layered Ocean Model (1/32°, http://www7320.nrlssc. navy.mil/global_nlom/) using the nctoolbox (https:// github.com/nctoolbox). We calculated the distance to eddy edge using mesoscale eddy trajectories [49]. Distance to productive seamounts or knolls, defined as those within 1500m of the sea surface, was calculated for each bird location, as these features are known to enhance biological productivity [50,51]. Monthly composites of MODIS-Aqua chlorophyll a, spanning 2007-2012, were constructed using DINEOF 3.0 to interpolate regions where clouds obscured satellite data [52]. Locations were matched to monthly composites, with locations from the first week of each month assigned to the previous month.
Wind speed was extracted from the RNCEP Reanalysis II data sets for surface values [53,54].

Residency time
We used residency time to understand the influence of habitat variability on individual movements. We chose this parameter rather than using a metric of foraging based on wet-dry data [32], as kittiwakes can employ both active and sit-and-wait foraging styles [55]. Residency time can be defined as the cumulative amount of time an individual animal spends within a circle of constant radius over a period of time [56]. In our case, we chose a radius of 45 km, as this is rougly equivelant to the mean daily distance kittiwakes traveled over the whole non-breeding period (September-May), and a time constraint of 1 month to avoid false positives if individuals crossed over their own path months later on the return trip. High residency locations, indicating periods of intense search effort, were then chosen as the upper quartile of each individual's residency times [57].

Habitat selection models
Linear-mixed models were used to relate residency time to oceanic habitat characteristics in the top four biogeographic regions frequented by kittiwakes. Variance inflation factors were calculated, but all were <3, so all explanatory variables were retained [58]. To meet the conditions of normality for linear models residency time and chlorophyll a were log transformed and distance to eddy and bathymetric slope were square-root transformed. Best-fit models for each ecosystem were constructed using a reductive approach and identified from Akaike Information Criterion (AIC) scores based on restricted maximum likelihood estimates [58]. Individual was included as a random effect. Marginal R 2 values were used to describe the variance explained by the fixed effects and conditional R 2 values, the combined fixed and random effects [59,60].

Results
Kittiwakes from the Pribilof Islands wintered predominantly across the deep oceanic waters of the central and western subarctic North Pacific all three years (Fig. 1). Area use of individual birds did not significantly differ between years, sexes, or colonies (1,557,000 ± 358,000 km 2 , p > 0.05; Table 2). As a group, females covered more area than males (Figs. 2 and 3). There was no difference in the size of cumulative area covered by birds from the two study colonies up to 40 individuals (n = 40,  (Fig. 3). The amount of overlap between groups was not different than random for colonies ( (Fig. 4).

Scale of Individual Spatial Fidelity
Birds that were tracked for two winters showed a tendency to return to the same regions and, for example, had on average 9 % overlapping grid squares for 100 × 100 km grid squares. During the mid-winter period (October -February) birds showed site fidelity at all scales, as the amount of shared grid cells was higher for individual bird repeat trips than for randomly paired tracks (Fig. 3). Some individuals repeated very unique routes, for instance, one bird traveled to the northern Bering Sea in the fall and then foraged along the Emperor Sea Mounts before heading back to the vicinity of the Pribilof Islands (Fig. 5).

Annual difference in large marine ecosystem use
Wintering Pribilof kittiwakes used seven biogeographical regions of the North Pacific (Table 3) (Table 3). The highest densities of tracked kittiwakes also occurred in the PSAW (Table 3). Annual differences were apparent in regards to when kittiwakes occupied each of the main ecoregions; specifically in October and November of 2009/10 there was a notable increase in the continued use of the Bering Sea (Fig. 3).
The change occurred because 70 % of birds remained in the Bering Sea in the fall, which corresponded with lower October sea surface temperatures than the other two years (Table 4). In both 2008/09 and 2009/10, Dec-Feb, over 30% of bird locations were in the NPPF; while in 2010/11 the majority of locations during these months were farther north in the PSAW (Fig. 6). On average individual kittiwakes used 4.9 ± 1 ecoregions and this did not differ between sex, or year, however only female kittiwakes traveled to the California Current and birds from St Paul used fewer ecoregions (4.6 ± 1) than birds from St George (5.3 ± 0.08) (Fig. 7).

Habitat conditions and selection
Throughout the winter, kittiwakes experienced a broad range of habitat conditions. Sea surface temperatures were on average 6.  (Table 3). Habitat variables were not a good predictor of residency in each ecoregion, with the best models only explaining <15% of the variation (Table 5).

Discussion
Kittiwakes from the Pribilof Islands underwent extensive pelagic migrations to diverse subarctic biogeographical regions in North Pacific. Both sex and past experience appeared to influence where individuals wintered. However, individuals displayed a large amount of flexibility and annual changes in distributions were larger than differences observed between sexes or colonies. As generalist predators kittiwakes are likely adapted to a high amount of environmental variability.

Interplay of intrinsic influences
The reasons for colony specific wintering areas are not always clear. In some cases these colony differences may arise as the result of differing local conditions altering the phenology of migrations or simply be the result of a range expansion from geographically separated colonies as birds are still tied to these at an annual or biennial time scale [15,30,46,61,62]. Alternatively there are examples where birds from different colonies winter in the same region [61,63]. The Pribilof colonies are only 70 km apart and the timing of breeding is highly synchronous [64]; thus it is not surprising that distributions of kittiwakes are similar. Yet, they are almost entirely distinct from wintering areas frequented by kittiwakes originating from Prince William Sound [33]. Therefore, there is likely additional colony associated variability in wintering areas across the North Pacific basin similar to that seen in kittiwakes in the North Atlantic [30].  We found a slight effect of sex on wintering distributions; though area use was largely similar, females appear to be more dispersive and no males traveled to the California Current System. Sex differences in distributions are often linked to differences in timing relative to nest defense or breeding roles and niche partitioning in sexually dimorphic seabirds [65][66][67]. Black-legged kittiwake males are larger than females in body and bill size [32] and this may lead to differences in energetic requirements or foraging preferences that may influence  [27,68,69], and that breeding outcome (or elevated stress levels) may carry over to affect the non-breeding distributions of one but not the other sex [29,70].
Fidelity during periods when individuals are constrained by central-place foraging, for instance in breeding seals and seabirds, appears to be relatively high [40,[71][72][73].
Much less is known about spatial fidelity when marine predators undergo migrations, often at the scale of ocean basins. Both flexibility and fidelity have been observed in migrating shearwaters, suggesting that individuals are able to explore and utilize multiple suitable wintering areas during a lifetime [24,25]. Like migrating Atlantic puffins [20], black-legged kittiwakes displayed a tendency to return to areas that they frequented the year before. Both puffins and kittiwakes  Marine biogeographical ecoregions are those defined by Longhurst [9] show characteristics of a dispersive migration, as routes often move gradually away from the breeding colony, rather than traveling to a distinct destination (see Fig. 4). However for kittiwakes, the amount of fidelity individuals showed was relatively small (~25 % of locations within 400 km grid squares) compared to the fidelity quantified for Atlantic puffins (median nearest neighbor distance <5 degrees) [20]. Kittiwakes are generalist predators that can only access the top 1 m of the water column, and they are much better fliers than puffins. Perhaps it is the combination of these life-history strategies that limits the amount of spatial fidelity that is advantageous, as prey resources can be patchy in pelagic environments [8]. Like puffins, it is unlikely that kittiwakes socially learn wintering areas as there is both individual fidelity and a high amount of variability within the population, instead it seems more likely that individuals rely on initial exploration, spatial memory and current conditions to inform travel paths. Our dataset is dominated by comparisons to 2009/10 when population level distributions were distinctly different from the other two years. This contrast likely led to less fidelity than might occur under similar conditions. More research is needed to understand if individuals are returning to areas where they previously experienced good foraging conditions and how fidelity changes when conditions change.

Oceanographic habitats and annual conditions
During all three of our study years the highest densities of wintering kittiwakes occurred in the central subarctic Pacific at the confluence of the Western Subarctic Gyre (PSAW) and the Eastern Subarctic Gyre (PSAE). This region, characterized by the subarctic current may support relatively high zooplankton biomass in the summer [74,75]. Kittiwakes from the Pribilofs used the entire PSAW. The PSAW has higher primary productivity than the Gulf of Alaska [76], more intense eddy activity [77], and supports a greater diversity of myctophid species [78]; potentially providing a more predictable winter habitat than the other ecoregions. Annual changes in kittiwake distributions occurred at a much higher magnitude, measured by the amount of overlap, than differences between intrinsic groups (colony and sex), and this reflects the contrasting conditions that occurred during the three-study years. Even restricted to ecoregions our models relating residency time and environmental variables were not able to explain much of this relationship. It may be, that as generalists, kittiwakes are foraging on a diverse suite of prey (with variable preferences for oceanographic conditions) thus limiting the ability of these models to identify strong environmental associations. Or, as we are limited to the resolution of geolocation derived predator data, meaningful habitatkittiwake associations maybe operating at smaller scales (both spatial and temporal), then we are able to quantify [28,79]. Regardless, overall our results show that  [80]. During this winter the western subarctic experienced anomalously high sea surface temperatures, while the central subarctic had anomalously cool sea surface temperatures [81]. The distribution of wintering kittiwakes in 2009/10 had a restricted longitudinal range and birds stayed much longer in the Bering Sea, potentially facilitated by colder sea surface temperatures (see Fig. 4 and Table 6), than during 2010/11, classified as strong La Niña year or 2008/09 classified as neutral conditions [82]. This shift to more northerly distributions in 2009/10 was also noted for wintering kittiwakes originating from the Shoup Bay colony in Prince William Sound [33]. This similar response suggests that during 2009/10 conditions suitable for wintering kittiwakes shifted northward across the North Pacific.
The winter of 2010/11, switched to one of strong La Niña conditions [83]. In the Gulf of Alaska, zooplankton biomass, survival estimates for age-1 pollock, and catch rates of juvenile pink salmon were all low [84,85]. In 2010/11, kittiwakes spent much less time in the Polar Front region, however birds that did venture there did not experience conditions that were notably different from either of the other years, except for a  All models include a temporal correlation term (corCAR1(form =~date|id)). Summary statistics of each full model (SST + depth + d2ed + slopetrn + EKE + ssh + d2hill + wind + chla) are presented first, followed by the best-fit model for each ecoregion. Akaike Information Criterion (AIC) were used to identify the best-fit model and marginal R 2 (R 2 (m)) and conditional R 2 (R 2 (c)) are presented. Abbreviations for the environmental variables used in the table are: sea-surface temperature (SST), distance to mesoscale eddy center (d2ed), sea-surface height (SSH), eddy kinetic energy (EKE), distance to productive seamounts and knolls (d2hill), monthly chlorophyll a (chla), bathymetric slope (slope) and bathymetry (bathy) marked decrease in EKE. Eddies and surface currents are known to condense and facilitate prey capture for surface foraging seabirds [28,79], thus is may be that this difference is linked to lower use. Relative to 2008/09 the more northerly distributions in 2009/10, when birds used the Bering Sea more, and in 2010/11 when kittiwake range in the subarctic was decreased, could be closer to what kittiwake wintering distributions may be like in the future. Climate change is predicted to shift the North Pacific Transition Zone farther north, causing the subarctic zone south of the Aleutians to shrink in size [86]. This shift in habitats is likely to be challenging for Hawaiian albatrosses as it moves preferred foraging areas farther Yearly means ± SD are calculated from individual bird means in each ecoregion from breeding colonies [87,88]. For kittiwakes, these changes will shrink the area of available wintering habitat, initially to a greater extent than what will open up to the north (due to the presence of land). Furthermore, in the North Pacific, Russian stocks of pink salmon (Oncorhynchus gorbuscha) currently follow a predictable alternation between large and small cohorts [89]. During large cohort years, prey availability is suppressed and across the North Pacific blacklegged kittiwakes respond through later hatch dates, lower laying success, smaller clutch size, and lower overall reproductive productivity-these effects appear to diminish in colonies in the Gulf of Alaska [89]. A shrinking subarctic may exacerbate competition between pink salmon, other salmonids, and seabirds. This could then lead to density dependent regulation of kittiwake populations, particularly if winter day length limits how far north kittiwakes can remain during the winter [90]. Though black-legged kittiwakes wintering in the pelagic North Pacific are able to retain spatial memories from one year to the next they also appear to have a high capacity to shift distributions in relationship to annual conditions. It remains unknown how long spatial memories may persist and if these long-lived individuals are simply returning to areas visited in the more distant past. It is also unknown if individuals are targeting previously profitable foraging areas. We were unable to find strong associations between kittiwake residency time and oceanographic habitats. This may reflect the scale of our analysis or be a characteristic that highlights the capacity of blacklegged kittiwakes to use multiple wintering ecoregions, presumably foraging on different prey species across both space and time.

Availability of supporting data
Raw geolocation data are available from the North Pacific Research Board [91].