Laysan albatross exhibit complex behavioral plasticity in the subtropical and subarctic North Pacific Ocean

: Animals that regularly traverse habitat extremes between the subtropics and subarc-tic are expected to exhibit foraging behaviors that respond to changes in dynamic ocean habitats, and these behaviors may facilitate adaptations to novel and changing climates. During the chick-provisioning stage, Laysan albatross Phoebastria immutabilis parents regularly undertake short-and long-distance foraging trips throughout the vast central North Pacific Ocean. We examined GPS tracking data among chick-provisioning albatrosses in Hawai‘i to characterize habitats during short-and long-distance trips. The study period encompassed a marine heatwave (2014) and the cooling period after an extreme El Niño event (2016), enabling us to examine foraging habitats under novel and changing climates. First passage time and generalized additive mixed models indicated that during 183 short and 110 long trips (n = 32 birds), wind-assisted flight efficiency, proximity to productive areas, and moonlit-searching were important in both subtropical and sub-arctic habitats. Laysan albatross took foraging trips that had similar lengths and durations in 2014 and 2016 and visited similar areas, indicating that their foraging range did not expand in response to climatic variability. A strategy that uses similar foraging areas across years combined with reliance on environmental processes that enhance flight efficiency (wind) and that enable searching behaviors (moonlight) indicate that Laysan albatross exhibit complex behavioral plasticity that allows them to utilize subtropical and subarctic habitats affected by dynamic climate variability. This strategy may benefit their ability to respond to oceanographic and climatic change, including expanding warm water regions and changing atmospheric conditions influenced by global warming.


INTRODUCTION
Animals that regularly traverse habitat extremes between the subtropical and subarctic marine environment could be expected to exhibit plasticity among behaviors and foraging strategies to respond to dynamic local conditions. In the subtropical central Pacific Ocean, albatrosses nest in Hawai'i but regu-larly travel to the subarctic gyre to forage (Kappes et al. 2010). This route traverses a > 20°C change in sea surface temperature (SST), several mesoscale atmospheric and oceanic features, and presumably areas with varying prey availability. The North Pacific Ocean exhibits basin-scale and mesoscale oceanographic processes and features that are predicted to shift with climate change, including expanded warm water regions (Le Borgne et al. 2011), expanded high-pressure systems (Sun et al. 2019), intensified climatic events (Cai et al. 2018), and shifts in atmospheric circulation patterns (Vecchi et al. 2006). Knowledge of species behaviors and how they use environmental conditions throughout variable habitats could help reveal the capacity for species to adapt to a changing climate and dynamic environments (e.g. Gilmour et al. 2018). In this study, we examined how a pelagic seabird, the Laysan albatross Phoebastria immutabilis (Hawaiian name: Mōlī), forages during the chick-provisioning stage in variable environmental conditions.
Laysan albatross parents that are rearing chicks undertake foraging trips in subtropical and subarctic waters to search for squid and other prey (Gould et al. 1997, Conners et al. 2018. Both parents return to their breeding colony at variable intervals until their chick fledges at 5.5 mo old (Awkerman et al. 2020). During the chick-provisioning stage, Laysan albatross and other related procellariform seabirds balance the costs of parental care and parental body condition by alternating short, near-colony foraging trips with longer, more distant trips to access preya strategy that is common among subtropical species (Baduini & Hyrenbach 2003). Temporal and energetic constraints require that chick-provisioning Laysan albatross parents use energetically efficient flight and environmental gradients (e.g. wind and waves) to minimize commuting costs while searching for food during this demanding life-history stage. Yet, the importance of oceanographic and environmental features to foraging albatross appear to vary among colonies throughout the North Pacific Ocean due to high intra-individual variation in movements (Gutowsky et al. 2015) and access to varying wind speeds, SST, and chlorophyll a (chl a) concentrations (Henry et al. 2021).
Coupling oceanographic features important to predators and prey can inform how far-ranging Laysan albatross navigate the dynamic ocean environment, especially throughout extreme habitats. Basin-wide and mesoscale environmental features can aid efficient travel and also transport prey (e.g. water density; Arkhipkin et al. 2015). For example, wind aids albatross flight (Suryan et al. 2008) and it also increases vertical mixing of the water column, which increases productivity and contributes to squid abundance (Nishikawa et al. 2015). Similarly, the position and seasonality of frontal zones and major currents define squid life history (Boyle & Rodhouse 2005) and are also important features of seabird foraging habitats (Scales et al. 2014). Identifica-tion of important habitat characteristics and the scales at which these features are used by albatross are therefore essential to understand how species can efficiently traverse the ocean, search for food, and exhibit plasticity in response to changing ocean conditions. Behavioral plasticity in foraging strategies is especially important for Laysan albatross as the species expands its range at sea and colonizes new nesting locations in the Pacific (e.g. Young et al. 2009b, Henry et al. 2021. Such behavioral plasticity is also important for individuals experiencing humanassisted social attraction or translocations to establish new breeding areas (Miskelly et al. 2009, Vander -Werf et al. 2019. In this study, we combined GPS tracking data with remotely sensed and model-derived oceanographic data to determine at-sea habitat characteristics that influence the foraging behaviors among chickprovisioning Laysan albatross. Previously, tracking studies of Laysan albatross have been limited to the chick-brooding and chick-guard stages  or had temporally low resolutions (sampling resolution: 2 locations d −1 ; Young et al. 2009a, Kappes et al. 2015. This is the first high-resolution study from the chick-provisioning stage that leads up to fledging, when chick flight-feather growth is fastest, and energy demands are greatest (chickprovisioning stage defined by age range of 17−165 d; Rice & Kenyon 1962). We conducted tracking during 2 breeding seasons that coincided with 2 climatic events: an anomalous marine heatwave and the cooling period after an extreme El Niño event -events that we suggest might be representative of future ocean conditions as our climate changes.

Field methods
We tracked chick-provisioning Laysan albatross from 3 colonies in the main Hawaiian Islands: 2 colonies in northeastern Kaua'i (NE Kaua'i 1 in 2014 and NE Kaua'i 2 in 2016) and 1 colony in northwestern O'ahu (Ka'ena Point in 2014). The colonies on Kaua'i were re-colonized during the 1970s, and the O'ahu colony was established in 1992 (Young et al. 2009b). Seabird nesting locations within island groups are often pooled for analyses. In this study, the Ka'ena Point colony was 130 km from the Kaua'i colony, but substantial inter-colony differences in area-restricted search (ARS; see Section 2.3) were not expected at this small scale. Chicks were 38−71 d old (mean ± SD: 59.1 ± 7.9 d, n = 33) when tags were deployed on parents during the post-guard stage of chick rearing (during late March to mid-April in both study years). We captured birds by hand at marked nests only if they were observed feeding a chick or were previously banded and a known parent. Nonbanded birds were fitted with a size 7B hard-metal leg band for identification (U.S. Geological Survey Bird Banding Lab). The sex of individuals was determined from blood or feather samples via DNA sexing (Fridolfsson & Ellegren 1999). Regular monitoring throughout the breeding season provided Laysan albatross chick hatch dates and ages. In 2014 (NE Kaua'i 1 and Ka'ena Point), e-obs Bird Battery 1Alight GPS/accelerometer 30 g tags (e-obs GmbH) were attached to several back feathers using cloth tape (Tesa 4651). In 2016 (NE Kaua'i 2), we taped a plastic template to birds' back feathers to which we attached e-obs Bird Solar GPS/accelerometer 30 g tags with plastic cable ties. GPS tags were 1.3% of the mean body mass of albatrosses in this study (mean ± SD mass: 2347 ± 274 g, n = 22 birds), which is less than the maximum of 3% of body mass recommended for tracking tag attachments on birds (U.S. Geological Survey Bird Banding Laboratory, https:// www.usgs.gov/labs/bird-banding-laboratory/science/ auxiliary-marking-authorizations). Tagged birds were returned to their nests within 10−20 min of initial capture. Tags were presumed to fall off back feathers when the cloth tape decomposed or when back feathers molted during the post-breeding period (Awkerman et al. 2020). Bird Solar tags included a solar panel that recharged the battery and Bird Battery tags relied exclusively on an internal battery. Tags were programmed to collect locations every 15 min in 2014 and every 5 min in 2016 (15 min when battery levels were low). Tags transmitted location data to a receiver that we installed in each colony to acquire and archive tracking data, thus absolving the need to recapture the birds to recover the tag.

Tracking data processing
We processed and analyzed all GPS data using the program R (version 4.0.2; R Core Team 2020). Tracks were filtered to remove locations within 2 km of the colony. A speed filter was applied to remove erroneous locations for which speed between 2 locations generated impossibly fast transit speeds (144 km h −1 ; 0.3% [1114 / 412 589] of raw data points were removed) based on the known biology and flight behaviors of albatrosses . Each bird's track was separated into individual foraging trips (function 'MakeTrips', R package 'trakR', version 0.0.9, Fleishman et al. 2020) that were ≥ 1 h in duration and that occurred at least 10 km (calculated via great circle distance) from the colony (Kappes et al. 2010). Because of differences in sampling rates between years, data were standardized by inter polating to 15 min intervals (function 'redisltraj', R package 'ade-habitatLT', version 0.3.25, Calenge 2006). Trips were manually examined for completeness (identified departure and return to colony) and only complete trips were used in descriptive statistics (trip length and duration) and in comparative analyses of these metrics between sexes and years. To examine the directions that Laysan albatross departed from and returned to the colony, outbound, inbound, and mid-trip locations were identified (function 'InOutPoints', R package 'trakR', Fleishman et al. 2020), and the mean longitude per outbound and inbound segment was calculated for each foraging trip. Preliminary examination of trip durations indicated a strongly bi-modal distribution ( Fig. S1 in the Supplement at www.int-res. com/articles/suppl/m697p125_supp.pdf); therefore, we examined short (<100 h) and long trips (≥100 h) separately in all subsequent analyses.

First passage time
ARS is often used to identify locations where animals spend greater proportions of time and enact tortuous paths that are often interpreted as searching or foraging behaviors (Fauchald & Tveraa 2003). To discern parts of trips where birds exhibited ARS behavior, first passage time (FPT) analyses were used (function 'fpt', R package 'adehabitatLT', Calenge 2006) on all trips > 4 h in duration, following Pinaud (2008). FPT is the time required to cross a circle of a given radius and this metric can provide information about foraging areas and/or behavior (e.g. scale ;Fauchald & Tveraa 2003), such that small FPT values indicate short searching times. FPT does not incorporate behavioral variability (e.g. does not differentiate between 'sit and wait' and foraging in flight strategies; Conners et al. 2015) but does enable identification of habitats in relation to small-scale and mesoscale environmental processes (e.g. Suryan et al. 2006). Non-overlapping ARS habitat (hereinafter 'ARS zones') was identified along each trip through an iterative process following Suryan et al. (2006). Briefly, FPT was calculated for all trips at 8 different radii separated by 30 km intervals (radii: 20, 50, 80, 110, 140, 170, 200, 230 km). These 8 radii encom-passed (1) the median radii (55 km) of FPT from all trips identified from plots of the maximum log variance of FPT vs. radii when FPT was calculated for 5−1000 km at 5 km intervals; and (2) the range of radii reported as important scales in other studies of procellariform seabirds (range: 29−121 km; Suryan et al. 2006, Pinaud 2008, Kappes et al. 2010. For each radius, the maximum FPT value was identified for each trip, and all locations within a circle of the specified radius (e.g. 20 km) were removed surrounding that location; this was considered the 'ARS zone'. Subsequent ARS zones within the trip were identified by using the nexthighest FPT value from the remaining locations, and locations within that 20 km circle were removed. This process was repeated until no additional points could be assessed without overlapping existing ARS zones, and then the process was conducted for each of the remaining radii. This process resulted in independent, non-overlapping ARS zones that included the 8 different radius scenarios for each trip. A scaled FPT value (FPT value for the ARS zone divided by the area of the circle for the radius of interest) was calculated for each ARS zone, following Kappes et al. (2010). These scaled FPT enabled comparisons between trip types and among birds.

Habitat variables
For each ARS zone, the median value of each environmental variable (Table 1) was calculated from publicly available rasters. SST, SST anomaly, chl a concentrations, and sea level pressure were obtained via the NOAA Coast Watch ERDDAP service (function 'rxtracto_3D', R package 'rerddapXtracto', version 1.0.1, Mendelssohn 2020). SST was obtained from the dataset 'Group for High Resolution SST (GHRSST) Level 4 MUR Global Foundation Sea Surface Temperature Analysis (v4.1)', with spatial and temporal resolutions of 0.01° and 1 d, respectively (https://doi.org/10.5067/GHGMR-4FJ04; JPL MUR MEaSUREs Project 2015). SST anomalies were calculated as daily SST values minus a 30 yr climatological mean and were obtained from the dataset 'SST, Daily Optimum Interpolation (OI), AVHRR Only, Version 2.1, Final, Global, 0.25°, 1981−present', with spatial and temporal resolutions of 0.25° and 1 d, respectively (Huang et al. 2021). Chl a concentrations were from the dataset 'High Resolution chlorophyll-a concentration' from MODIS/Aqua (https://coastwatch. pfeg.noaa.gov/infog/MW_sstd_las.html; NOAA Coast -Watch Program and NASA Goddard Space Flight Center 2020), with spatial and temporal resolutions of 0.025° and 1 mo, respectively. Sea level pressure was obtained from the Fleet Numerical Meteorology and Oceanography Center (https://www.usno. navy.mil/FNMOC/) 'Sea Level Pressure 360X180' dataset, with spatial and temporal resolutions of 1.0° and 6 h, respectively. Depth was obtained from the NOAA National Centers for Environmental In formation ETOPO1 database (Amante & Eakins 2009), with a spatial resolution of 0.0166° (function 'get-NOAA.bathy', R package 'marmap', version 1.0.6, Pante & Simon-Bouhet 2013). Wind speed was calculated from hourly uand v-vectors of wind at a spatial resolution of 0.25° and with a reference height of 10 m from the European Centre for Medium-Range Weather Forecasts ERA5 Reanalysis (Hersbach et al. 2018) that were downloaded from the Copernicus Climate Change Service (https://cds.climate. copernicus.eu/). Wind speed was calculated with the equation: . To describe wind directions relative to birds' headings, the relative angle between the bird heading and wind along the trajectory was calculated with custom code (Methods S1). Filament activity was described by the AVISO product backward-in-time finite-size Lyapunov exponent (FSLE) derived from the dataset 'SSALTO/Duacs delayedtime global ocean absolute geostrophic currents' ('DUACS2018 DT MADT UV products') and downloaded from Aviso Satellite Altimetry Data via the Centre National d'Études Spatiales (https://www. aviso.altimetry.fr/es/data.html) at spatial and temporal resolutions of 0.04° and 1 d, respectively. Lagrangian coherent structures (LCS) were identified as |FSLE| values > 0.1 d −1 , following Tew Kai et al. (2009). Particles with maximal FSLE values tend to form continuous lines that guide transport processes throughout the world ocean and can thus identify areas with high dispersion rates and the occurrence of stretching fluid parcels, which may aid travel and foraging for marine predators (Tew Kai et al. 2009, Abrahms et al. 2018. In addition to the median FSLE value per ARS zone, the proportion of |FSLE| > 0.1 d −1 was calculated across each ARS zone, and this is referred to as the 'proportion of LCS'. The percent moon illumination was calculated with the function 'moonAngle' (R package 'oce', version 1.2.0, Kelley & Richards 2020). Percent moon illumination was calculated for each ARS zone by multiplying moon altitude with either 0 or 1 that represented day (defined by astronomical day phase, calculated from the angle of the sun with the function 'sunAngle', where sun angles > −18° were considered 'day') or night (sun angles <−18°), respectively. This resulted  (Kappes et al. 2010(Kappes et al. , 2015 Squid have water temperature preferences (Rowell et al. 1985, Arkhipkin et al. 2015 Sea level pressure Temporal: 6 h Spatial: 1.0° (~111 km) hPa Related to seasonal manifestation of the North Pacific High and Aleutian Low, which affect frequency and location of storms and wind speed and direction (Peterson et al. 2017), which can affect albatross flight (Suryan et al. 2006) High pressure = low wind = less water mixing = low food (Robinson et al. 2016) Wind in a metric of moon illumination for nighttime when the moon was above the horizon (positive moon altitude value), but zero values during daylight hours or when the moon was below the horizon (negative moon altitude value). We also calculated distances between Laysan albatross ARS zones and oceanographic frontal zones. At the largest scale, the North Pacific Ocean can be divided latitudinally by the 18°C SST isotherm which coincides loosely with the subtropical frontal zone (STFZ; Bograd et al. 2004) and the 0.2 mg m −3 chl a contour, which is frequently used to demarcate the transition zone chlorophyll front (TZCF) and the southern extent of subarctic waters (Polovina et al. 2001). This study occurred in the spring−summer, and we refer to the area between the 18°C isotherm and the 0.2 mg m −3 chl a contour as the North Pacific transition zone (NPTZ). Distances to the STFZ and TZCF were calculated separately by constructing a raster of the isotherm or contour for each tracking day and calculating the great circle distance between the center point of each ARS zone and the closest point along the isotherm or contour.

Statistical analyses
Because all birds had multiple recorded trips, we used linear mixed effects models (function 'lmer', R package 'lme4', version 1.1.29, Bates et al. 2015) with individual bird ID as a random factor to examine whether trip length, trip duration, or mean longitudes of inbound and outbound trips differed between years or sexes. One bird whose sex was unknown was excluded from analyses in which sex was a factor. Because birds were tagged at 3 separate colonies, we also used linear mixed effects models to ask whether trip length and maximum distance traveled from the colony differed among the 3 colonies. Type 3 sum of squares were then used to assess the significance of year, sex, or colony from the linear mixed effects models (function 'Anova', R package 'car', version 3.0.13, Fox & Weisberg 2019). Chisquared tests (function 'chisq.test', R package 'stats', version 4.1.2, R Core Team 2020) were used to assess whether SST anomalies in ARS zones were evenly distributed or whether Laysan albatross disproportionately foraged in areas that had SST anomalies either greater than or less than zero. To assess the degree of similarity in ARS zones between 2014 and 2016, we calculated an overlap index (Bhattacharayya's affinity, BA; Bhattacharyya 1943) for all points from short and long trips separately. This index results in a value ranging between 0 (no overlap) and 1 (complete overlap; function 'kerneloverlap', R package 'adehabitatHR ', version 0.4.19, Calenge 2006). To increase the robustness of our comparisons of spatial overlap, we also tested for significant differences between the observed BA values and a null distribution generated from 1000 randomizations for short and long trips, separately (function 'kernalOverlapBA_p', R package 'trakR'; Fleishman et al. 2020). A significant difference from this test occurs if < 5% of randomized BA values do not exceed the observed BA values. For short trips, we used a grid size with 1 km resolution, a 5 km buffer, and an h-value of 11.64. For long trips, we used a grid size with 10 km resolution, a 25 km buffer, and an h-value of 195.05.
Generalized additive mixed models (GAMMs) were used to assess the relationships between scaled FPT (hereinafter 'FPT s ') and environmental factors for short and long trips, separately. The distributions of FPT s had heavily skewed tails and the distributions of environmental factors were varied and, when plotted against FPT s , had linear and non-linear relationships, so we used GAMMs for location, scale, and shape (GAMLSS; R package 'gamlss', version 5.4.3, Rigby & Stasinopoulos 2005). GAMLSS models location (mean; μ), scale (σ), and shape (ν, τ) parameters for response variables that may include high skewness and kurtosis (e.g. Hernandez et al. 2009).
To build GAMLSS models, we first assessed the best distribution that described the response variable (FPT s ) with the function 'fitDist'. The skew power exponential type 2 ('SEP2'; Azzalini 1986) and sinharcsinh ('SHASH') distributions (Jones & Pewsey 2009) were chosen for short and long trips, respectively, based on relative Akaike information criterion (AIC) values.
Because environmental factors were on different scales and had uneven distributions, they were centered (the mean was subtracted from each value) and scaled (each value was divided by the SD of centered values) with the function 'scale' (R Core Team 2020). All trips with duration > 4 h (Pinaud 2008) were used in GAMLSS models except for 1 bird whose sex was unknown (n = 6 trips) and short trips > 4 h but for which FPT s could not be calculated with the focal radii (n = 7 trips). The variables distance to the 18°C SST iso therm (STFZ) and distance to the 0.2 mg m −3 chl a contour (TZCF) were only used as predictors in GAMLSS models for long trips because the short trips did not extend beyond 27° N and did not traverse either frontal zone. In addition to environmen-tal factors, the fixed categorical factors sex, year, and FPT s radius and the random effect individual bird ID were included in GAMLSS models. Because short trips mostly did not encompass the largest radii of interest, only the 4 smallest radii were included in GAMLSS models for short trips.
We assessed the potential for uninformative parameters with tolerance values (the inverse value of variance inflation factors, a measure of collinearity) using a linear mixed effects model (function 'lmer', R package 'lme4', Bates et al. 2015) on all candidate variables with individual bird ID as a random effect. Variance inflation factors were calculated with the function 'vif.mer' (Frank 2014) and predictor variables whose tolerance values were < 0.1 were removed from subsequent GAMLSS models (Quinn & Keough 2002). Because we conducted separate GAMLSS models for short and long trips, different variables were removed from each full model. For short trips, the predictors FPT s radius and proportion of LCS were removed. For long trips, the predictors SST and distance to the 0.2 mg m −3 chl a contour (TZCF) had tolerance values < 0.1 and were removed.
We then fit separate parametric models for short and long trips with the function 'gamlss'. Parametric models included all remaining predictor variables and individual bird ID as a random intercept. We included a Gaussian correlation structure to account for spatial autocorrelation among ARS zones with the terms latitude and longitude. To assess assumptions for parametric models, we visually checked residual plots for normality and independence, and we also checked that the residual distribution had a mean near 0, variance near 1, coefficient of skewness near 0, and a coefficient of kurtosis near 3 (DeCarlo 1997). Because GAMLSS models may include terms for location, scale, and shape of the response variable as a function of the predictor variables, models were built with increasing complexity, starting with the location parameter (μ). To better meet parametric model assumptions, additional models were constructed by including formulas for the distribution parameters of variance (σ) and shape, represented by skewness (ν) or kurtosis (τ), as needed. Sigma formulae were added to short-and long-trip GAMLSS models and a nu formula was added to the long-trip GAMLSS model. The terms chosen for sigma and nu formulae were based on plots of each predictor variable against the response variable: plots in which the response variable had a high fluctuation at higher ranges of predictor variables were included in the sigma formula, and plots in which the response vari-able was skewed with respect to the predictor variables were included in the nu formula. For short trips, variables included in the sigma formula were depth, FSLE, sea level pressure, SST, and wind speed; for long trips, these variables were distance to the 18°C SST isotherm (STFZ), FSLE, proportion of LCS, and sea level pressure. The terms chosen for the nu formula were chl a, depth, distance to the 18°C SST isotherm (STFZ), FSLE, proportion of LCS, percent moon illumination, sea level pressure, and wind speed. This process resulted in 2 and 3 parametric models that were generated for short and long trips, respectively, that included (1) only μ terms; (2) μ and σ terms; and (3) μ, σ, and ν terms (long trips only). Generalized AIC (function 'GAIC', R package 'gamlss', Rigby & Stasinopoulos 2005) was used to compare relative model fits to choose the final parametric model; within this function, the argument 'k' was defined as: log(sample size).
We then fit reduced separate non-parametric models for short and long trips to derive a more unbiased estimate of important predictors. To do this, smoothing terms (penalized beta-splines) were used with the function 'pb' for the significant μ predictors from the parametric models; non-significant parametric terms were retained in the model. The smoothing function was applied to each term with the default number of knots (n = 20 knots, which is recommended for large sample sizes; Rigby & Stasinopoulos 2005) and the smoothing parameter λ was determined with the method 'GAIC', where the argument 'k' was defined as the log of the sample size for short and long trip models, separately. The default number of knots for percent moon illumination during long trips resulted in an overfit model that produced a marginal effect plot that was not biologically meaningful. So, we restricted the degrees of freedom (df) for percent moon illumination for long trips to 4.28, which corresponded to the effective degrees of freedom (edf) generated for this term by the final nonparametric model for short trips. We report estimated coefficients and p-values for all factors and edf for smoothed factors in final non-parametric models (Zuur et al. 2009) and Cox-Snell pseudo-R 2 values (Smith & McKenna 2013) for all models.
The maximum distance traveled from colonies during long trips was significantly different among colonies (ANOVA: χ 2 = 9.452, p = 0.009) such that trips were farther from Ka'ena Point (mean ± SD of maximum distance: 3173.4 ± 649.3 km, n = 22 trips, 11 birds) than from NE Kaua'i 1 (mean ± SD of maximum distance: 2715.7 ± 511.0 km, n = 32 trips, 11 birds) and NE Kaua'i 2 (mean ± SD of maximum distance: 2661.6 ± 659.0 km, n = 50 trips, 11 birds). The Ka'ena Point colony is only 70 km farther south than the Kaua'i colonies, but these birds traveled 400−500 km farther north into the Aleutian Islands and the Gulf of Alaska, where they used the continental shelf break near Kodiak Island (Fig. S2). The maximum distance traveled during short trips was not different among colonies (ANOVA: χ 2 = 0.255, p = 0.880). Total trip distances were not different among colonies for both short (ANOVA: χ 2 = 0.038, p = 0.981) and long trips (ANOVA: χ 2 = 2.317, p = 0.314).

Foraging habitats
Chick-provisioning Laysan albatross took distinct short and long trips to areas within 12.1−4283.0 km of their nests. Long trips traversed distinct habitats that were delineated by large ranges of productivity, SST, FSLE, wind speeds, and sea level pressures (e.g. Fig. 3); long trips crossed both the STFZ and TZCF (Fig. S3), but short trips remained within the subtropics to 27°N (Fig. 1A) and did not approach either front. Long trips ranged to at least 45°N, north of the TZCF that delineates the southern reaches of the subarctic; some birds traveled into the Gulf of Alaska and near the Aleutian Islands, >1000 km beyond the TZCF. The mean latitude of the STFZ shifted seasonally from 30°N in March to 40°N by July, and the TZCF shifted slightly north and became more variable throughout summer (Fig. S3).
Laysan albatross encountered winds from all directions during long trips, where wind direction changed with latitude, resulting in mostly quartering tailwinds relative to birds' flight paths (Fig. S4). Easterly winds occurred near Hawai'i, and wind direction at locations > 30°N (near the STFZ) shifted to westerly winds. During short trips, albatross primarily encountered winds from the east and northeast, resulting in crosswinds relative to birds' flight paths (Fig. S5).
3.2.1. Environmental predictors of FPT s FPT s ranged from 26.4 min to 11.8 d across the 8 radii (Table 3). Multiple ARS scales were evident within trips (Fig. 4), such that variance in FPT s was greatest for small radii (e.g. 20 km) and variance was lower for larger radii (e.g. 110 km) during short and long trips. FPT s was generally 1.5× longer during short trips compared with long trips for radii at 20, 50, 80, and 110 km, but there was greater variation in FPT s during long trips (Table 3).
3.2.1.1. Short trips. Significant terms from the final GAMLSS model for FPT s during short trips were chl a, depth, percent moon illumination, sea level pressure, SST, wind speed, sex, and year (Table 4). During short trips within ~700 km of the colony, FPT was longest at moderate chl a concentrations, shallow depths, greater amounts of moonlight, high sea level pressure, low and high SST, and at low wind speeds (Fig. 5). Females had longer FPT s than males, and FPT s was longer in 2016 than 2014. 133 Fig. 2. Laysan albatross flew in anticyclonic loops during long trips. Boxplots of mean longitude of outbound and inbound trip segments per foraging trip (n = 102 long trips, n = 201 short trips, 32 birds; data do not include bird whose sex was unknown). Boxplots represent the median (thick vertical line), upper and lower quartiles (leftmost and rightmost limits of the box, respectively), and minimum and maximum values (horizontal lines and dots) of each subset of data. *Mean longitudes of outbound and inbound trip segments were significantly different in linear mixed effects models. There were no statistically significant differences between sexes and years, so data are pooled for plots  The final selected model that described FPT s during short trips was: FPT s ~ FSLE + pb(chl a) + pb(depth) + pb(% moon illumination) + pb(sea level pressure) + pb(SST) + pb(wind speed) + Sex + Year + re(random = ~1|Bird ID) Sigma formula = ~depth + FSLE + sea level pressure + SST + wind speed Correlation = corGaus(form = ~ Latitude + Longitude) where FPT s from non-overlapping ARS zones, re represents a random effect for individual bird, and pb is a penalized beta spline.

Long trips.
Significant terms for the final GAMLSS model for FPT s during long trips were chl a, distance to the 18°C isotherm (STFZ), FSLE, pro-portion of LCS, percent moon illumination, sea level pressure, wind speed, radius, sex, and year (Table 5). There was an overall increase in FPT s with increasing sea level pressure (especially at the greatest sea level pressure values), an increase with greater distance to the 18°C SST isotherm (STFZ), and a cyclic, bell-shaped curve for FPT s during nights with greater moon illumination that reached a maximum at intermediate values of illumination (Fig. 6). There was an overall decrease in FPT s with increased wind speed and decreased proportion of LCS (Fig. 6) (Humphries 2021) centrations ( Fig. 6B). At higher concentrations, FPT s generally decreased with increasing chl a, although the latter part of the trend included a large outlier (Fig. 6B). Additionally, FPT s generally increased with increased search radius, but radii of 200 and 230 km were not significant. FPT s was significantly longer in 2016 compared with 2014, and females had significantly longer FPT s than males (Fig. 6). The final selected model that described FPT s during long trips was: FPT s ~ depth + pb(chl a) + pb(distance to SST isotherm) + pb(FSLE) + pb(proportion of LCS) + pb(% moon illumination) + pb(sea level pressure) + pb(wind speed) + Sex + Year + radius + re(random = ~1|Bird ID) Sigma formula = ~distance to 18°C SST isotherm + FSLE + proportion of LCS + sea level pressure Nu formula = ~chl a + depth + distance to 18°C SST isotherm + FSLE + proportion of LCS + % moon illumination + sea level pressure + wind speed Correlation = corGaus(form = ~Latitude + Longitude) where FPT s is from non-overlapping ARS zones. To better understand how the ARS scale changed with habitat variables, we evaluated an additional GAMLSS model that included the interaction between radius and each smoothed continuous environmental covariate. Most variation in FPT s among radii was caused by wind speed (Fig. S6, Table S1). Regardless of radius, there were overall in creasing relationships with chl a and sea level pressure, and overall decreasing relationships with FSLE (Fig. S6). The model fit was similar for the model with and without radius interactions (Cox-Snell pseudo-R 2 = 0.27 and R 2 = 0.25, respectively), but GAIC results indicated the simpler, more parsimonious model without radius interactions to be the better model (GAIC, non-radius-interaction model: −238 007.0, df = 126.8; radius-interaction model: −235 222.9, df = 462.5).  Laysan albatross occurred in anomalously warmer waters than expected by chance during the marine heatwave in 2014 (long trips: χ 2 = 18 311, df = 1, p < 2 × 10 −16 ; short trips: χ 2 = 876.96, df = 1, p < 2 × 10 −16 ; Fig. 7). Laysan albatross occurred in anomalously cooler waters than expected by chance during long trips in 2016, during the cooling period after the extreme El Niño event (χ 2 = 14437, df = 1, p < 2 × 10 −16 ), but they occurred in relatively warmer waters than expected by chance during short trips (χ 2 = 681.85, df = 1, p < 2 × 10 −16 ; Fig. 8).

DISCUSSION
Laysan albatross traversed subtropical and subarctic environments to find food for self-maintenance and chick provisioning. Across 2 non-consec-utive years, Laysan albatross visited similar locations. Laysan albatross, considered squid specialists, enacted ARS at multiple scales throughout these trips, especially at radii 20−170 km, indicating use of multiple environmental features to find food. Important processes included large-scale environmental features that optimize albatross flight efficiency (wind, sea level pressure) and small-scale features that may help optimize prey searching (SST, FSLE, chl a, moonlight, depth). ARS occurred in anomalously warmer waters during the marine heatwave in 2014. In 2016, ARS occurred in anomalously cooler waters during long trips and in anomalously warmer waters during short trips. Use of similar locations between years could indicate that Laysan albatross have the capacity to adapt to variable environmental conditions. Such adaptation may benefit Laysan albatross cope with increasing environmental heterogeneity as the climate changes.

Environmental processes that optimize flight efficiency increase habitat accessibility
Large-scale phenomena like wind and sea level pressure fields, and proximity to frontal zones and filaments, may aid flight-and foraging efficiency. Consequently, albatrosses may increase foraging habitat accessibility (habitat that is detectable, identifiable, reachable, and available to albatross throughout the breeding-season foraging range) to enable chickprovisioning adults to exploit food resources 1000s of kilometers from their nest. For example, wind facilitates energetically efficient flight in albatrosses (Weimerskirch et al. 2000), and Laysan albatross benefit from high wind speeds via fast transit rates (Thorne et al. 2016), which could result in the inverse relationship between FPT s and wind speed observed in Laysan albatross in this study. High wind speeds facilitate a dynamic soaring flight pattern, where albatrosses gain energy from horizontal wind movement in 4 phases (gliding along waves, turning into headwinds, ascending while gaining kinetic energy, gliding while descending back to sea surface; Sachs 2005). Wind may also differentially enable efficient flight between sexes, resulting in segregated habitats based on wind speed (Shaffer et al. 2001). Male Laysan albatross were 7.4% heavier on average than females in this study , and typically, males have higher wing-loading, meaning that they are adapted to make more efficient use of greater wind speeds to fly compared with females (e.g. Shaffer et al. 2001, Wakefield et al. 2009). These differ-139 Fig. 6. (A) Marginal effect plots for the partial effects of chl a, distance to the 18°C sea surface temperature (SST) isotherm (subtropical frontal zone, STFZ), finite-size Lyapunov exponents (FSLE), proportion of Lagrangian coherent structures (LCS), percent moon illumination, sea level pressure (SLP), wind speed, radius, sex, and year on scaled first passage time during 110 long trips by 32 chick-provisioning Laysan albatross. Continuous variables were centered and scaled before analyses and are unitless. Black lines: smoother estimated with penalized beta splines (pb) of the continuous variables in the model; grey shading: 95% confidence intervals. Distributions of individual scaled values represented by rug plots (dark grey vertical bars). (B) Marginal effect plot for the partial effect of chl a with a modified y-axis to illustrate the trend more clearly in lower values; thick black line: smoother; thin black lines: 95% confidence intervals ences may help explain the shorter FPT s in males compared with females. The inverse relationship between wind speed and sea level pressure also helps to describe spatial patterns of Laysan albatross foraging trips. Dynamically soaring albatrosses are expected to avoid highpressure systems, where wind speeds are low (Weimerskirch et al. 2000), but also may be 'trapped' (such that they cannot fly) in areas with very low wind speeds (Suryan et al. 2006). Laysan albatross FPT s appeared to reflect this paradoxical pattern; individuals spent more time (longer FPT s ) in areas with higher sea level pressure (potentially trapped), but they also partially avoided the high-pressure regions by undertaking anticyclonic (clockwise) outbound trajectories during long foraging trips. This stereotypical pattern would allow foraging albatross to avoid the seasonal North Pacific high-pressure system located to the north of Hawai'i and enable albatross (via more efficient transit) to reach the NPTZ and higher latitudes, which are characterized by greater productivity and greater wind speeds (Adams & Flora 2010).
Frontal zones and LCS at the boundaries of water masses may also provide predictable foraging opportunities for diverse predator assemblages at large scales (Scales et al. 2014). Thermal properties like SST can also distinguish areas where water masses and aggregated biota converge (Bograd et al. 2004). During short trips in this study, FPT s was longer at low and high SST, and Laysan albatross avoided intermediate values of SST. Laysan albatross could focus search efforts in areas that aggregate prey with either convergent, relatively warm surface water or divergent, upwelled, relatively cool water. However, FPT s increased with greater distance to the STFZ, indicating that albatross spent more time engaged in ARS in subarctic waters with more primary productivity. The STFZ, despite having greater proportions of LCS, appeared less important to far-ranging Laysan albatross. These patterns also indicate that while mesoscale fronts did help predict Laysan albatross foraging habitats (Kappes et al. 2010, Thorne et al. 2015, the absolute values of FSLE and distance to the STFZ may simply represent the latitude at which albatross enact ARS during long trips. Persistent features like zones with greater wind speeds and convergent water masses at northerly latitudes may combine to provide predictable features with environmental characteristics that increase habitat accessibility and profitability. For example, the number of LCS in the central North Pacific are greatest at 35° N (Abernathey & Haller 2018), but few ARS zones were observed at this latitude. In fact, the greatest proportion of ARS zones occurred north of the TZCF at ARS radii of 20 and 50 km, and the greatest number of larger ARS zones with radii ≥ 80 km occurred between fronts, in the NPTZ. These ARS patterns align with previous observations that described increased foraging within the NPTZ (Hyrenbach et al. 2002, Kappes et al. 2010, where wind speeds are greater and likely beneficial for albatrosses (Suryan et al. 2008). However, there may be an upper limit to the scale of utility of mesoscale features to foraging Laysan albatross: in our study, only radii up to 170 km were determined to be important for describ-ing FPT s and ARS, limiting the importance of larger, basin-wide features for describing smaller-scale behaviors at sea.

Laysan albatross use small-scale environmental processes to optimize ARS
Mesoscale features could help Laysan albatross identify and reach optimal foraging regions, and at smaller, sub-mesoscales, moonlight, chl a, and depth may combine to help Laysan albatross locate food. For example, seabirds likely use a combination of vision, sound, and olfaction to locate prey (Nevitt et al. 1995, Bolin et al. 2009. Laysan albatross forage both during the day and at night, when moonlight likely enhances visual acuity (Conners et al. , 2018. Laysan albatross FPT s was greatest during periods with intermediate values (during long trips) and greatest values (during short trips) of moon illumination. During the lunar phases that correspond with intermediate and sub-maximum illumination levels, nocturnally migrating squid and other prey may approach the surface more frequently than during maximal moon illumination, when they might be more exposed to predators  (Benoit-Bird et al. 2009a,b). The response of FPT s to moonlight could therefore represent opportunistic foraging behaviors that occur approximately monthly; albatross likely forage regardless of lunar phase (e.g. Pinet et al. 2011), but if their foraging trip coincides with intermediate and high amounts of moonlight, they may benefit from increased light that facilitates visual foraging, and from potentially enhanced prey availability closer to the sea surface. Localized productivity could also be associated with Laysan albatross foraging habitats. Longer FPT s coincided with intermediate concentrations of chl a and, during short trips, with shallower depths. Marine productivity (indexed by chl a concentration) is important for describing habitat associations among ommastrephid squids (Boyle & Rodhouse 2005, Ichii et al. 2011 and gonatid squids in the subarctic gyre (Kearney et al. 2013), and for characterizing areas used at sea during short and long trips for several procellariiform seabirds (Baduini & Hyrenbach 2003), including Laysan albatross (Henry et al. 2021). Locations with persistent productivity (see Suryan et al. 2012) could help explain why the relationship between FPT s and chl a was maximal at intermediate values, especially among short trips (our Fig. 5; Kappes et al. 2015). Although dynamics affecting surface chlorophyll concentrations are complex, intermediate chl a concentrations could reflect areas with elevated relative productivity and presence of grazing zooplankton, which constantly recycle chlorophyll at and below the sea surface (see Strom et al. 2001). Dimethyl sulfide produced by zooplankton grazing on phytoplankton could generate olfactory foraging cues for seabirds that extend up to 4 km away (Nevitt et al. 1995). Static bathymetric features also provide predictable foraging habitats at sea because shallow features like seamounts aggregate biota (De Forest & Drazen 2009). Aggregations of predators may attract other seabirds by sight and sound up to several kilometers away (Bolin et al. 2009, Thiebault et al. 2016). Thus, actively seeking out these habitat cues could represent a reliable strategy employed by Laysan albatross during both short and long trips.

Interpretations of FPT benefit from consideration of scales
Because FPT is a metric of time, it could be interpreted in 2 ways: Did an albatross spend a long time (high FPT) in an area because it was searching for prey, or because the environmental conditions were unfavorable for flight? This apparent dichotomy can be limiting and confusing, and it is important to assess FPT within the context of expected scales of behaviors and other factors like environment and physiology. Consideration of environmental contexts, especially the scale at which environmental processes occur, and the scale at which processes may influence flight or searching behaviors, is important when interpreting FPT values. For example, FPT s was inversely related to wind speed (relatively coherent at 100s to 1000s of km), which we interpreted as wind aiding albatross flight. By placing this result within the context of the known flight physiology of albatrosses (Weimerskirch et al. 2000), a high FPT value would indicate unfavorable flight conditions if only wind conditions were considered. Conversely, high FPT s values were associated with intermediate chl a concentrations. Chlorophyll can vary at scales of 10s of km, but chl a does not affect albatross flight physiology. An interpretation of high FPT in relation to chlorophyll could indicate optimizing search behaviors. Although these are contrasting interpretations, statistical approaches like GAMMs that assess FPT in relation to environmental conditions can account for additive factors and help clarify interpretation. Despite the complexity of these interpretations, however, FPT information is valuable and enables examination of large-scale patterns and environmental features. Whereas analyses that identify discrete behaviors (e.g. expectation-maximization binary clustering, EMbC; hidden Markov movement models) work well for high-resolution data with temporal resolutions ≤ 5 min (Mendez et al. 2017, Conners et al. 2021, FPT is meaningful for analyses where behaviors related to larger-scale processes, and where tracking data sampled at coarser scales (e.g. 15 min) are of interest.
Consideration of the scales at which ARS activity occur are important to the interpretation of albatross' interactions with prey and their environment. The scales at which albatross identify and interact with environmental features are typically > 20 km (Suryan et al. 2006, Kappes et al. 2010. Finer spatial scales (100 m to 10 km) that occur over seconds to 10s of minutes can provide valuable insight about smallscale behaviors, like drift-foraging or resting periods . However, inclusion of periods when albatross sit on the water could inflate the variance in FPT at very small scales, and these periods are sometimes removed prior to analyses to reduce this error . Drift periods (which could encompass rest or drift-forage behaviors) were previously identified using GPS data collected at very short intervals (10 s; Weimerskirch et al. 2007. Our sample intervals (15 min) did not allow for identification of these fine-scale behaviors or distinction between resting and drift-foraging. Therefore, we did not remove rest periods prior to analyses. Because albatross-environment interactions occur at an order of magnitude greater than drift behaviors (20−200 km vs. 100 m to 10 km), the retention of drift periods in our analyses likely did not affect the variance in FPT s at radii in which we were most interested. Small-scale drift movements could indicate use of local visual and olfactory cues associated with seamounts ). Availability of GPS data at higher temporal resolutions would have enabled a more comprehensive examination of the importance of moonlight, chl a, and depth within our smallest ARS zones (20 km), but even at a scale of 20 km, these apparent patterns enabled identification of important ARS habitats.

Changing climates
Laysan albatross benefit from increased foraging habitat accessibility through a combination of processes that optimize flight efficiency at large scales and that aid searching behaviors at smaller scales. The predictability of some of these features may increase site fidelity and the ease with which they can repeatedly detect and reach important foraging areas. With climate change affecting the entire Pacific, mesoscale features including sea level pressure, wind, and SST, and sub-mesoscale features including chl a concentrations and fronts are predicted to change in ways that may adversely affect Laysan albatross' use of these features while foraging. A weakening of the Walker Circulation (Vecchi et al. 2006) and a weakening of the Aleutian Low (Sun et al. 2019) are expected to result in stronger and more expansive seasonal high-pressure systems in the North Pacific that could decrease the area at sea available to Laysan albatross that benefit from moderate to strong winds to maximize flight efficiency. An accompanying change in wind patterns would further reduce efficient access to large areas at sea and potentially force increases in the lengths and durations of foraging trips (e.g. Suryan et al. 2008). Throughout large areas of the central North Pacific Ocean, climate change may result in the expansion of warm water (Collins et al. 2010, Polovina et al. 2011, and important subarctic foraging areas will be located farther from nesting sites (Thorne et al. 2015) or could overlap with fisheries (e.g. Wren et al. 2019). Among fishes and certain albatross prey, some species could shift poleward and species composition could homogenize across larger gradients (Magurran et al. 2015). For example, optimal SST for squid spawning habitat (21−25°C; Ichii et al. 2009) could move poleward with expanded warm waters, farther from Laysan albatross nesting colonies. Changes in regional prey composition could also result in less spatial overlap with prey (e.g. Sadykova et al. 2020) and potentially could result in less optimal trophic niche shifts (e.g. Bond & Lavers 2014).
Two climate events occurred during this study that presage altered foraging conditions for albatross in response to climate change. An anomalous marine heatwave with SST up to > 4°C warmer than average occurred throughout the eastern North Pacific during 2014 (Bond et al. 2015). The heatwave extended as far west as 160° W and as far south as 30° N (Peterson et al. 2017) during the 2014 tracking period (April−June), covering one-quarter of Laysan albatross foraging areas (Fig. 7). In this region, there were mixed heatwave effects on the marine foodweb, including immediate and lagged changes (up to 5 yr; Suryan et al. 2021). In 2015−2016, a strong El Niño resulted in SST up to 2.5°C warmer than average, weak trade winds, and low nutrient and high oxygen concentrations in the upper 200 m throughout the eastern tropical Pacific (Stramma et al. 2016). The 2016 tracking period corresponded with the end of this event-a transition period during which multivariate ENSO index (MEI) values decreased from 1.3 to −0.5 and represented cooling SST (https://psl. noaa.gov/enso/mei/). Although seabirds may benefit from reduced thermocline depth and in creased water mixing during cooling SST periods (Ribic et al. 1992), months-long lagged relationships between SST and chlorophyll could have longer-term consequences on foodwebs (Kislik et al. 2017). FPT was shorter in 2014, indicating short search times, where Laysan albatross concentrated ARS in relatively warmer waters. In 2016, albatross spent more time engaged in scale-specific ARS that may have resulted from lower prey availability or changes in wind conditions (e.g. Thorne et al. 2016). Although Laysan albatross are considered squid specialists, foraging on epipelagic squid at night and opportunistically on micronekton (e.g. floating dead squid; neuston) during the day (Conners et al. , 2018, individuals may have some capacity to adapt to changes in prey availability and composition. In response to marine heatwaves, specialist seabird species expand foraging ranges and may reduce resting periods (Osborne et al. 2020), whereas generalist seabird species experience less physiological stress than specialists (Tate et al. 2021). Laysan albatross took trips that had similar lengths and durations in 2014 and 2016 and visited the same areas, indicating that their foraging range did not expand in response to climatic changes in SST. However, changes in prey composition may have other effects: some prey species experienced decreased growth and lipid accumulation during the marine heatwave (von Biela et al. 2019), which could affect predators throughout the foodweb (Piatt et al. 2020). Although diet was not assessed in this study, limited observations from the Laysan albatross that we studied on Kaua'i had a relatively high chick-fledging rate in both 2014 (0.85 [17/20] chicks fledged; Kaua'i Albatross Network and USFWS unpubl. data) and 2016 (0.79 [169/198] chicks fledged). These fledging rates are greater than the mean fledging rate for Kaua'i overall (0.74 ± 0.15 SD during 2015−2021), indicating that foraging conditions were sufficient despite anomalous SST conditions. Fledging success on Kaua'i was also high in 2015 (0.88 [217/247] chicks fledged), indicating no lagged effects on fledging were observed during the year following the marine heat wave. Although reproductive success was high, increased variability in SST (Cai et al. 2018) may force Laysan albatross to adjust to greater habitat variability within their vast foraging range. Foraging strategies that enable Laysan albatross to cover a large region may integrate enough habitat variability to take advantage of small-scale pockets of productivity, enhancing behavioral plasticity.

Conclusions
While far-ranging albatross traverse 1000s of kilometers across dynamic ocean habitats, they experience varying scales of environmental processes that modify habitat accessibility and likely affect foraging efficiency. As individuals search this environment for food, albatross respond and modify their searching behaviors and engage in ARS. Environmental processes and ocean features help Laysan albatross locate foraging areas and find prey within the patchy subtropical and subarctic central North Pacific Ocean. The potential for Laysan albatross to find and track these features at different scales indicate that Laysan albatross are flexible foragers during an energetically and temporally constrained breeding stage. Such behavioral plasticity may benefit Laysan albatross as they colonize new sites and explore new at-sea habi-tats (Henry et al. 2021). However, several of the features on which albatross rely are predicted to change in the future as the planet warms (Collins et al. 2010, Polovina et al. 2011, Sun et al. 2019. Reduced wind speeds and a stronger and more expansive North Pacific high-pressure system could affect albatross flight efficiency, changing how and where they access foraging habitats. Increased SST could also affect prey assemblages that are sensitive to changes in SST (e.g. Magurran et al. 2015, Piatt et al. 2020, Donahue et al. 2021. Variable external factors, combined with persistent use of similar lo cations between years, indicate there may be a limit to how quickly and how well Laysan albatross can adapt to future environmental changes. Laysan albatross currently demonstrate incredible adaptations to changes in environmentsalternating foraging trips between subtropical and subarctic habitats; but the degree to which they can adapt to future changes remains unknown.