Reproductive constraints influence habitat accessibility, segregation, and preference of sympatric albatross species

Study area and tracking activities

Breeding Laysan and black-footed albatrosses were studied at Tern Island (23.87° N,
166.28° W), French Frigate Shoals, Northwest Hawaiian Islands, during the incubation,
brooding, and chick-rearing periods. The incubation period for both species lasts
approximately 60 days, beginning in mid-late November with the laying of a single
egg and ending with the onset of hatching in late January. The brooding period typically
lasts 2–3 weeks, until chicks are left alone at the nest in mid-February. The chick-rearing
period lasts approximately 4–5 months; after a period of fasting, chicks fledge independently
in June and July 33], 34].

We used a combination of satellite tracking and light-level based geolocation to determine
at-sea locations of adult Hawaiian albatrosses throughout the reproductive period.
Satellite tags were used for short-term deployments during incubation and brooding,
whereas geolocation tags were used to obtain foraging positions during chick-rearing,
and to supplement data from the incubation period. While geolocation tag deployments
also spanned the brooding period, the single position derived per day by light-level
based geolocation did not have sufficient resolution to accurately capture movements
during brooding when average trip durations were 2–3 days (Table 1; satellite tags are accurate to 10 km 53]; geolocation tags are accurate to ~200 km 54]). Tracking activities were conducted during five consecutive breeding seasons, from
2002–03 through 2006–07; satellite tracks were obtained during the 2002–03 to 2005–06
seasons and geolocation data were obtained during the 2003–04 to 2006–07 seasons.
Sex of tracked individuals was determined from blood samples 55]; for six individuals for which we did not obtain blood samples, sex was determined
by comparison of culmen lengths 28].

Table 1. Summary characteristics (Mean?±?SD) of Laysan and black-footed albatross foraging
trips. To reduce the influence of individuals tracked for multiple foraging trips,
a single trip per individual was randomly-selected for each reproductive stage to
include in the calculation of mean values

Satellite tracking

One hundred and forty-seven adult albatrosses (76 Laysan and 71 black-footed) were
equipped with satellite platform terminal transmitters (30 g Pico-100, Microwave Telemetry,
Columbia, MD; or 35 g SPOT4, Wildlife Computers, Redmond, WA); tags were attached
to dorsal feathers with adhesive tape (tesa®, Hamburg, Germany), and transmissions
were downloaded via the Argos satellite system (Service Argos, Inc., Largo, MD). Satellite
tags were programmed to transmit continuously every 90 s, with the exception of 19
foraging trips, when tags were programmed to use a 6:18 h on:off duty cycle (11 trips),
20:4 h on:off duty cycle (6 trips), or 9:15 h on:off duty cycle (2 trips) to conserve
battery life on anticipated longer trips. In three instances, the battery on the satellite
transmitter failed or the satellite transmitter fell off the bird before completion
of the foraging trip. During 2002–03, 2004–05, and 2005–06, satellite-tracked individuals
were also equipped with archival tags (10 g Lotek LTD 2400, Lotek Wireless, St. John’s,
Newfoundland) attached to a plastic leg band so that temperature recordings (±0.05 °C)
every 12–40 s could be used to characterize foraging activity while at sea 56]. The combined mass of devices deployed on individuals was less than 2 % of total
bird body mass in all cases, below the recommended limit for studies involving albatrosses
57].

Before calculating trip characteristics, satellite locations were first delimited
by observations of departure and arrival times at the breeding colony or by visual
inspection of the tracks. To remove unlikely locations from the data set, we applied
an iterative forward/backward averaging speed filter 58] implemented in Matlab (The MathWorks, Natick, MA) with a maximum speed limit of 80 km h
^?1
(following 29], 59]) to remove unrealistic flight speeds 60]. Satellite tracks were then interpolated to every 10 min using a Bézier curve with
??=?0.3 [? controls the elasticity of the curve; following 61], and subsampled to two locations per day to match the temporal scale of geolocation
positions (see below).

Tracking with geolocation tags

Archival geolocation tags (10 g Lotek LTD 2400, Lotek Wireless, St. John’s, Newfoundland)
were deployed and successfully recovered from 34 Laysan and 26 black-footed albatrosses.
Each tag recorded ambient light intensity and temperature every 480 s or 540 s to
determine a single daily location: longitude based on the establishment of local noon
in comparison to Universal Time, and latitude based on day length for the established
longitude 62], 63]. Temperature sensors on the geolocation tags allowed a refinement of the location
data based on sea surface temperature [SST; 54], 64], as well as providing a record of foraging activity while at sea.

Geolocation tags were deployed for up to one year, only a portion of which overlapped
our on-colony research activities, therefore we had limited information on the presence
or absence of geolocation-tagged individuals at the colony. Due to the nature of the
SST-processing algorithm 54], 64], it was necessary to first delimit individual foraging trips so that on-colony locations
were removed from the dataset. We calculated minimum daily temperatures recorded by
geolocation tags to determine days when albatrosses were likely at sea; evaporative
cooling due to immersion in water would lead to cooler minimum temperatures on days
when albatrosses were off-colony. In order to detect shifts in minimum daily temperatures
indicative of on- or off-colony periods, we used an algorithm developed for detecting
climate regime shifts 65], which implements a sequential version of the partial CUSUM method combined with
the t-test, and is available for download at: http://www.beringclimate.noaa.gov/regimes/. For those individuals for which exact colony attendance patterns were known (25
individuals during incubation), this detection algorithm correctly classified 99.0 %
of days when birds were known to be on-colony and 83.6 % of days when birds were known
to be at-sea (n?=?548 bird days). This method was therefore effective in delimiting foraging departures,
but some at-sea locations may have been excluded. Estimated departure and arrival
times for geolocation-derived foraging trips were calculated based on Argos-derived
transit rates during the incubation period 32], and distance to the colony from the first and last off-colony location.

To remove unlikely light-based positions, we applied the same speed filter as above,
with a more conservative speed limit of 50 km h
^?1
to account for greater error in geolocation position estimates 54], and interpolated to two positions per day using a Hermite spline [following 61].

Foraging movements

Descriptive characteristics of all Argos- and geolocation-derived tracks were calculated
to compare foraging behavior during each reproductive stage. Maximum distance traveled
from the colony was calculated using great-circle distances to account for the earth’s
curvature. Destination bearings, as defined by the azimuth to the most distant point,
were calculated to describe overall direction of foraging trips. Trip duration was
defined as either the time elapsed between the observed departure and arrival of the
bird, or the time elapsed between the estimated departure and arrival times.

We used linear mixed-effects models 66], with individual as a random effect, to compare maximum ranges and trip durations
between species and reproductive stages. Mixed-effects models were followed by contrast
analysis with the multcomp package 67] in the program R 68]. To compare destination bearings between species and breeding stages, we used Watson-Williams
tests for circular data in the CircStat toolbox in Matlab 69] after ensuring data followed the von Mises distribution (equivalent to normal distribution
for circular data). Accounting for random effects has not been developed in the circular
statistics modeling framework, therefore we randomly selected a single trip for each
individual during each reproductive stage for this analysis.

Foraging distribution

To determine patterns of interspecific habitat segregation during each reproductive
stage, we used kernel estimation 70] to determine utilization distributions (UD) for each species. Because individuals
contributed varying numbers of foraging trips to the overall data set, a single trip
was randomly selected for each individual during each reproductive stage so that the
influence of individual behaviors on estimation of kernel densities was reduced. Geographic
coordinates of interpolated albatross locations were transformed using a Lambert Cylindrical
Equal Area projection 71], and UDs were computed on a 50-km grid using the R package adehabitat 72]. To allow comparisons between species, the smoothing factor (h) was set to the mean
of the h values calculated from each species, as determined using least-squares cross-validation
73]. We then employed a randomization analysis to test the null hypotheses that there
was no spatial segregation in foraging distributions between species during each reproductive
stage 74]. For each comparison, species was randomly assigned to tracks using the same species
ratio observed, and kernel analysis applied. The area of overlap divided by the area
of the larger of the two UD polygons was used as the test statistic following Breed
et al. 74], for the 25 % (core area), 50 % (focal region), and 95 % (foraging range) UDs 29]. Each test was iterated 500 times, and the P-value was calculated as the proportion of random overlaps smaller than the observed
overlap 74].

Activity patterns at sea

We determined the proportion of time albatrosses spent in flight and the frequency
of landings on the sea surface during daylight and nighttime hours to characterize
foraging activity patterns during each reproductive stage using geolocation tags equipped
with temperature sensors. Albatrosses are surface-feeders and, by necessity, must
land on the sea surface in order to consume prey. Therefore, landing rates are indicative
of the level of feeding effort. Flight costs measured for albatrosses have demonstrated
that the most energetically demanding activities albatrosses engage in at sea are
take-offs and landings 75]; landing rates were also highly correlated with field metabolic rates in a study
of wandering albatrosses 76].

For albatrosses equipped with temperature recorders, we implemented an algorithm (Iknos
toolkit for Matlab; Y. Tremblay, unpublished) designed to identify landings based
on rapid changes in temperature, and stable periods associated with sitting on the
sea surface 56]. We used civil twilight (sun no more than 6° below the horizon) to define daylight
hours, based on temporally-matched tracking locations and NOAA’s solar calculator,
as implemented in the maptools package in R 77] .

To test whether percent time in flight and landing rates differed between high-resolution
temperature records (12–40 s) and low-resolution records (480 s and 540 s), we subsampled
high resolution records to the lower resolution and compared these metrics. While
percent time in flight did not differ between high resolution and subsampled records
(paired t-test: t₁₄₀
?=?1.11, P?=?0.27), landing rates were significantly lower in the subsampled group (paired t-test:
t₁₄₀
?=?20.9, P??0.0001). To compare all temperature records in subsequent analyses (only low resolution
records were available during the chick-rearing period), we used histogram matching
to rescale landing rates from low resolution records 78]. Histogram matching is an image processing technique used to rescale lower resolution
data so that the histogram of the data after transformation matches that of the reference
data 79]. This allowed us to present landing rates comparable to values in the literature
determined from high resolution temperature loggers; significance tests based on histogram-matched
landing rates yielded the same conclusions as tests based on landing rates from all
the data subsampled to the lowest resolution available.

Percent time in flight was arcsine transformed and landing rate was log transformed
prior to analysis to meet normality assumptions. We used linear mixed-effects models
66] with individual as a random effect to compare percent time in flight and landing
rates between species and reproductive stages, followed by contrast analysis with
the multcomp package in R 67].

Habitat preference

We followed the analytical framework of Aarts et al. 80] to model habitat preference of Laysan and black-footed albatrosses during each stage
of breeding. We adopted a case–control design such that each tracking location was
temporally matched with three randomly generated control locations. Because breeding
albatrosses are central place foragers, it is unrealistic to assume that all points
within the study area are equally accessible 80], 81]. Therefore, we adopted a simple null model of usage that assumes that the accessibility
? 81] of a point in space is inversely related to the distance from the colony (d
_c
) 82]. We assumed that locations were not accessible beyond the maximum range observed
for each breeding stage for each species; control locations were then quasi-randomly
selected at a rate proportional to ? within this range for each species-stage combination.
Because this null model may over- or under-predict true accessibility, we also included
d
_c
as a candidate covariate in our habitat preference models 83] as suggested by Aarts et al. 80].

Laysan and black-footed albatrosses breeding at Tern Island rarely made southward
departures from the colony (Fig. 1); therefore we also modeled habitat preference using a more restrictive null model
of usage to reflect the northern bias of tracking locations (see Appendix A in Additional
file 1) and to ensure conclusions were not sensitive to the choice of null usage model.
Final habitat models and response curves of selected covariates were generally similar
irrespective of the null usage model (see Results and Appendix A in Additional file
1), therefore we present only the results based on the simpler null usage model, which
assumed that locations were accessible within the maximum range observed for each
species-stage, regardless of the direction from the colony. This required fewer assumptions
regarding accessibility of habitats, and followed the implementation of Wakefield
et al. 82].

Fig. 1. Tracking locations of Laysan and black-footed albatrosses breeding at Tern Island,
Northwest Hawaiian Islands. Data encompass five breeding seasons, from 2002–03 through
2006–07, for Laysan albatrosses during incubation a brooding b and chick-rearing c and black-footed albatrosses during incubation d brooding e and chick-rearing f

Habitat preference models were implemented using a binomial distribution and logit
link (the inverse of the logistic function) to relate the response (tracking locations,
value of 1 (Fig. 1) and control locations, value of 0 (Fig. 2)) to environmental covariates. Geolocation tracks (incubation and chick-rearing)
were only included in habitat analyses if they included at least five filtered off-colony
locations, to ensure adequate coverage of sampled habitats along each track. We applied
generalized additive mixed models (GAMMs) to allow for the possibility of a nonlinear
response to environmental covariates 84], and to account for non-independence of points within trips and the variable number
of trips contributed by each individual 85].

Fig. 2. Randomly generated control locations selected at a rate proportional to accessibility. Maximum observed ranges of Laysan and black-footed albatrosses during the incubation
a brooding b and chick-rearing (c) periods were used to limit spatial extent of control locations

Environmental covariates

We selected marine habitat variables based on their potential to characterize physical
features that may stimulate or aggregate albatross prey resources. Because Hawaiian
albatrosses are surface-feeders foraging within the top meter of the water column,
they rely on neustonic or vertically migrating prey, as well as on carrion and fisheries
discards 33], 34], 86], 87]. Hawaiian albatrosses consume a diverse array of prey items, with Ommastrephid squid
and flying fish (Exocoetidae) eggs comprising the largest proportion of the Laysan
and black-footed albatross diet, respectively 86]. The diversity of prey types consumed by Hawaiian albatrosses reflects the variety
of marine habitats encountered during their long-ranging movements; for instance,
flying fish eggs are likely unavailable to black-footed albatrosses when foraging
in the temperate waters of the California Current, where squid may be a larger component
of their diet.

Previous studies have demonstrated a relationship between catches of Ommastrephid
squid and SST 88]–90], productivity 91], and position relative to frontal features of the North Pacific Transition Zone (NPTZ)
90]. Areas of surface convergence, such as those found near fronts and eddies, can also
aggregate algae or floating material which many species of flying fish use to attach
their non-buoyant eggs with sticky filaments 92], 93]. Adult flying fish, also prey of Hawaiian albatrosses 86], can be found at the outer edges of rapidly rotating eddies 93], 94].

The habitat relationships of these prey items informed our selection of environmental
parameters with which to investigate albatross habitat use. We used SST and primary
productivity (PP) to characterize the regional thermal and phytoplankton production
regimes; latitudinal gradients in SST and proximity to the Transition Zone Chlorophyll
Front (TZCF) to describe large-scale frontal characteristics; wind stress curl to
describe wind-driven oceanic upwelling or downwelling; and sea surface height anomaly
(SSHa) and eddy kinetic energy (EKE) to investigate overall intensity of eddy activity
95]–97]. As albatrosses use near-surface winds to engage in gliding flight 98], we also characterized wind speed and direction for all locations used in the habitat
analysis. Because Hawaiian albatrosses periodically visit the continental margins
of the western coast of North America and the Aleutian Islands during breeding 28], 29], 32], and bathymetric features can stimulate local production and aggregate prey 99], 100], we also used sea floor depth to characterize marine habitats; bathymetry data was
extracted from NOAA’s ETOPO2 data set (http://www.ngdc.noaa.gov/mgg/global/etopo2.html).

Environmental data were obtained by querying the NOAA OceanWatch Live Access Server
using Matlab and ERDDAP (http://coastwatch.pfel.noaa.gov/erddap/). Where possible, we used satellite products with the finest temporal resolution
available (daily or weekly); we used monthly composites when necessary to compensate
for data gaps due to cloud cover.

We used a blended product of SST derived from both microwave and infrared sensors
carried on multiple platforms at a spatial resolution of 0.1 degrees and as a 5-day
temporal composite 101]. As SSTs in the central North Pacific are distributed in broad latitudinal bands,
latitudinal gradients in SST (dySST) were used to describe frontal structure 102], by computing the local derivative of adjacent 0.1-degree pixels in the north–south
direction.

Monthly PP estimates 103] were derived at a spatial resolution of 0.1 degrees from monthly chlorophyll-a values and photosynthetically available radiation obtained from the Sea-viewing Wide
Field-of-view Sensor (SeaWiFS) on the Orbview-2 satellite, and SST values from the
AVHRR Pathfinder Oceans Project.

Monthly chlorophyll-a values obtained from SeaWiFS were also used to calculate distance to the TZCF. Following
Polovina et al. 46] and Bograd et al. 50], the TZCF was defined as the 0.2 mg m
^?3
chlorophyll-a contour.

Daily composites of wind velocity fields were gathered by the SeaWinds scatterometer
aboard NASA’s QuikSCAT satellite at a resolution of 0.25 degrees and a reference height
of 10 m, and used to assign daily values of wind speed and direction to each location.

The NOAA CoastWatch West Coast Node also processes wind vector fields to calculate
wind stress curl from these data. We used wind stress curl as a metric of surface
convergence/divergence 50], as a proxy for local aggregation of neustonic or buoyant prey (e.g., fish eggs,
fish and squid larvae, zooplankton, fishing discards). For each location of interest,
we calculated a mean wind stress curl from the previous 30 days to get a representation
of past physical forcing in the area of interest.

EKE was calculated from surface ocean currents derived from four satellite altimeters
(Jason-1, ENVISAT, ERS-1 and 2, and TOPEX/Poseidon) provided by the AVISO program
at a resolution of 0.25 degrees. For each location of interest, a 10 x 10-degree spatial
mean was removed from each 7-day composite of zonal and meridional surface ocean currents
to generate zonal and meridional velocity anomalies (u’ and v’). The EKE (per unit mass) was then calculated as

(1)

SSHa was derived from AVISO sea level data compared to the mean geoid as measured
from 1993–1995.

Oceanographic data were extracted so that composite data matched closest in time to
tag-derived locations (and matching control locations). The median value of each oceanographic
variable was calculated for grid cells falling within the approximate error of Argos
or geolocation positions. To estimate error of Argos locations, we conducted trials
on seven satellite tags affixed to permanent structures at the University of California
Santa Cruz, Long Marine Lab, for a period of three weeks. Argos locations were then
compared to GPS-derived locations; mean error for all location qualities was 0.06°
(6.7 km) from a total of 2,585 location fixes. To be conservative, satellite tag-derived
locations (and matching control locations) were extracted within a 0.15° (16.7 km)
longitude by 0.15° (16.7 km) latitude grid centered on each at-sea location. For geolocation-derived
positions, data were extracted within a 1° (111 km) longitude by 2° (222 km) latitude
grid centered on each location [the approximate error of the geolocation method; 54].

Model fitting and selection

We modeled habitat preference for each species and breeding stage considering the
following candidate covariates: SST (°C), PP (mg C m
^-2.
day
^?1
), dySST (°C km
^?1
), minimum distance to the TZCF (dTZCF; km), wind stress curl (Pa m
^?1
), SSHa (cm), EKE (cm
²
 s
^?2
), wind speed (m s
^?1
), wind direction (°), sea floor depth (m), and d
_c
(km). EKE and PP were log-transformed to achieve an even spread of covariate values,
thus avoiding undue leverage of a few high values 104]. To ensure that there was not strong collinearity among parameters, we used generalized
linear models with a binomial distribution and logit link to relate the response variable
(presence/control) to all covariates, and calculated tolerance values (inverse of
the variance inflation factor) for each species-stage; all tolerance values were greater
than 0.1, the approximate guide suggested by Quinn and Keough 105].

GAMMs were implemented within the gamm4 package in R 106]. All candidate covariates were fitted using cubic regression splines; a cyclic spline
was used to model wind direction. We first investigated the potential effects of transmitter
duty-cycle and tag type (geolocation versus satellite tracking), as well as year and
sex effects, by comparing intercept-only GAMMs for each species and reproductive stage
to models including each of the above terms. In all cases, Akaike’s information criterion
(AIC) values were not improved by the addition of these terms, therefore they were
not subsequently included in models of habitat preference.

To arrive at an inferential model for each species and breeding stage, we used AIC
as a guide, and relied on cluster-level cross-validation for final model selection
80]. This approach was taken due to the fact that model selection based on AIC alone
can lead to overparameterized models when data are spatially or temporally autocorrelated,
as is expected with both tracking and oceanographic data 80]. Therefore, we randomly selected two-thirds of the individuals tracked during each
breeding stage for each species to use in the first stage of model selection. We started
with intercept-only models and used forward model selection with AIC to arrive at
a suite of candidate models for each species-stage. We then fit data from the remaining
individuals to the sequence of models obtained from forward selection. The model with
the lowest AIC from the cross-validation step was retained as the inferential model.
To reduce the chance of over-fitting, we then replaced each spline in turn by a linear
term and selected a final inferential model for each species-stage based on AIC values.
Candidate GAMMs were fit using maximum likelihood to allow comparison of models with
different fixed effect structures 107]; final inferential models were fit using restricted maximum likelihood estimation
108]. We used the proportion of the total deviance explained by the model to assess goodness
of fit 107], and the percent deviance contributed by each coefficient to assess how much variability
in the response could be explained by each main effect 109].