The intensity of horizontal and vertical search in a diving forager: the harbour seal

Data sampling

Data on the horizontal and vertical movements of harbour seals were collected from
14 animals belonging to a resident population in the Porsangerfjord, northern Norway
(70-71° N, 25-26° E, Fig. 5). The animals were equipped with GPS phone tags (SMRU Instrumentation, University
of St. Andrews, UK) glued on their fur in the shoulder area right after the moulting
period (July) to ensure the maximum longevity of the sampling (the tags would fall
off during the next moult). Six animals were caught in September-October 2009, five
in September 2010, one in August 2011 and two in September 2012. The animals’ handling
procedures (capture, transport, tagging, and release) were approved by the Norwegian
Animal Research Authority and are described in Ramasco et al. 36].

Fig. 5. The study area and its bathymetry. The study area (main plot and rectangle in lower
right plot) comprises an entire fjord system in Northern Norway: the Porsangerfjord.
The fjord’s bathymetry (grey scale) ranges between very shallow areas in the inner
western part of the fjord to deeper areas in the outer parts

The tags were set to recover the animal’s GPS position at 20 min intervals. Due to
the changing availability of the satellites or the tags being at times underwater,
registrations were occasionally delayed, resulting in irregular time series. Time,
retrieved from an onboard clock, and depth, measured through a pressure sensor, were
recorded regularly at 4 sec intervals and stored in the form of time-depth profiles
of 11 inflection points equally spaced in time. The maximum depth of each dive was
also recorded. A conductivity sensor detected at any time if the animal was underwater
or at surface. If the tag was dry for longer than 10 min, a haul-out start was registered,
which ended when the tag was wet for more than 40 sec. Data were temporarily stored
in the tag memory and later relayed through the GSM network.

Horizontal movement data and foraging indices

GPS data were filtered to retain only good quality positions (maximum error of 50 m,
49]). The irregular time series of GPS positions were cut into separate bouts if no position
was available for 24 h at sea or 48 h hauled out. Only the bouts of duration??3 h
were used in the analysis. Switching state-space models were fitted to the irregular
locations for each individual 37]. Two states (or movement types, MT) were allowed, assumed to correspond to fast directional
movements (extensive search, MT?=?0) or slow and tortuous movements (intensive search,
MT?=?1). From the model, horizontal speed (HS) was predicted at regular 20 min intervals.
Residence time (RT) was calculated from the predicted regularized locations as described
in Barraquand Benhamou 6]. The index corresponds to the time an animal spends within a circle of a given radius
(r) centred on each point along the trajectory. More precisely, RT is equal to the
time elapsed from the moment the animal enters the circle to the moment it leaves
it for longer than a given time threshold (t). RT values within a radius (r) distance
from haul-out site were excluded from these models since biased by the time of residence
at the haul-out site.

Vertical movement data and foraging indices

Errors in the registration of the seals’ vertical movements could arise due to missed
surfacing registrations caused by failure of the pressure or conductivity sensors
to detect respectively shallow depths or dry conditions. If one or more surfacings
were missed, multiple dives were compressed into one single 11-points profile causing
implausibly long dive durations at the limit of the tag’s registration capacity (25 min)
and largely above expected maximum dive durations for the species (ca. 10 min, 50]). Potentially incorrect dive records were defined as the ones that, excluding the
first descending and last ascending phases, showed one or more depth readings in the
upper 25 % of the dive without surfacing or when the animal stayed in the upper 25 %
of the maximum depth for more than 50 % of the dive duration. These thresholds were
chosen assuming that an animal is neither likely to decide to dive again after reaching
the upper part (upper 25 %) of the water column without surfacing, nor prone to spend
over half the time submerged right below the surface without emerging before or after
a deeper dive. Ninety percent of the dive records with durations at the upper edge
of the distribution, but not detected by the method previously explained, were less
than 5.6 m deep. This indicated that, most likely, for dives shallower than that threshold,
surfacings were often missed by the sensors. Dives shallower than 5.6 m were therefore
excluded from further analysis.

Bottom time (BT) was computed as the time in the bottom phase of each dive. We used
the definition of bottom phase as in Austin et al. 17] as the time spent in the lower 15 % of the dive’s maximum depth. BT was then standardized
across depths (stBT) by transforming it into a % of maximum potential BT (maxBT) for
a given dive depth and duration. For dive i:

(1)

(2)

Minimum travel time was defined as the time the animal would use to reach the depth
of the bottom phase (15 % of maximum depth) from surface at maximum vertical speed
(set as the 0.95 quantile of the individual’s distribution of vertical speeds, mean
1.97 m/s, range 1.75 – 2.16 m/s across individuals). For comparison, the Time at Depth
index (TAD, 29]) was also calculated and the correlation between BT, stBT, and TAD investigated [see
Additional file 2]. Since in general several dives occur between two successive locations (i.e. a trajectory segment), dive characteristics were averaged for each trajectory segment.

Covariates

Among the factors potentially affecting the relationship between horizontal and vertical
foraging indices, we considered the following variables, potentially affecting the
searching intensity of the animals while foraging: dive depth, trip direction, and
predatory tactic (benthic or pelagic). Moreover, we considered the presence of resting
behaviour while diving as a potential confounding signal for the detection of foraging
using indices based on the allocation of time in space (see Table 1 for a description of the covariates).

To differentiate between periods of benthic and pelagic diving behaviour, we calculated
for each dive??5.6 m the distance between its maximum depth and the depth of the
sea bottom, expressed as the depth of the water column at mid tide (modelled bathymetry,
grid cell 100 × 100 m, data from the Norwegian Mapping Authority 51]). Each dive was located on the bathymetric map by assuming constant swim speed on
a straight line between two successive GPS locations. To account for the variance
in the estimated distance to the sea bottom, which is due to the combined errors in
dive location, bathymetric predictions and tidal state, we fitted a mixture of n normal distribution functions (1??= n??=5) to the frequency distribution of bottom distances and modelled the probability
of each dive to belong to any of these distributions. From the best model (n?=?4, Fig. 6) we assumed the distribution having its mean closest to zero to be the distribution
of bottom distances for benthic dives. The mean and upper (95 %) quantile of that
distribution were found to be respectively –2.8 m and 2.2 m, therefore all dives for
which bottom distance was??= 2.2 were considered benthic. The negative average distance
from the sea bottom for benthic dives suggests that the tagged individuals were diving
on average in deeper waters than predicted. This bias, as mentioned above, was partially
due to errors in locating dives and in the bathymetry predictions. However, given
bottom depths were calculated for mid tide, the negative mean distance from the bottom
can partially indicate that seals were diving more often during high tide than low
tide (tidal range ±1 m).

Fig. 6. Detection of benthic and pelagic diving. The histogram shows the empirical distribution
of the distance between maximum dive depth and sea bottom, while the lines show the
fitted mixture of n normal distributions (best model for n?=?4). The fitted distribution with mean closer to 0 was assumed to be the distribution
of benthic dives (mean – 3.15 m). The dashed line shows the threshold (2.2 m) used
for distinguishing between benthic (= 2.2 m) and pelagic dives (2.2 m)

A categorical variable (Direction) was created to reflect the major movement direction
with respect to haul-out sites in each foraging trip. The different categories distinguished
between movements directed away from or towards the haul-out sites, movements close
to a haul-out, movements of transition between haul-out sites, or movements with no
particular directionality. To classify the different parts of a trip, first the difference
in distance to haul-out site was calculated for each position in a movement trajectory,
then a running average was computed (window width 6 h) to smooth the fine scale variation
in direction changes (Fig. 7a), and then each trip divided into sections whenever average distance difference changed
sign (Fig. 7B). ‘Outward’ (and ‘inward’) trip sections were defined as the continuous succession
of positions with a persistent increasing (decreasing) distance from the last (to
the next) haul-out site and starting (ending) within 2 km from the haul-out site [see
Additional file 6 for the choice of window width and threshold distance parameters]. Sections entirely
within the threshold distance from haul-out sites (2 km) were classified as ‘within
range’. Segments starting within the range of one haul-out site and ending within
the next were classified as ‘transiting’. The remaining segments, showing relatively
stationary movement behaviour away from haul-out sites, were by exclusion classified
as ‘other’. The approach used was able to detect multiple return trips between two
successive haul-out events (Fig. 7). We assumed the main trip directions to reflect the satiation state of the animals,
with a higher degree of satiation during ‘inward’ trip sections rather than ‘outward’
and ‘others’. This assumption was based on the biology of the species, a short range
forager and generally performing foraging trips relatively close to haul-out sites
35], allowing to assume that the individuals will on average return to the haul-out site
for resting purposes after foraging 34]. The categories ‘within range’ and ‘transiting’ (respectively 21 and 4 % of points)
were not considered to be related to a particular satiation state or to foraging behaviour
in general and were therefore not included in the analysis.

Fig. 7. Trip direction. The categorical variable (Direction) was defined as the major movement
direction with respect to the previous and next haul-out sites (black vertical bands).
The distance from the last and to the next haul-out sites (respectively solid and
dashed line) was calculated for each trajectory location (a). The difference between successive distances was then computed and averaged across
a running window of 6 h (b). Trip segment limits were defined when average difference distance crossed zero
(grey line in b). ‘Outward’ (medium dark grey bands) and ‘inward’ (light grey bands) trip segments
were identified as the segments respectively starting and ending within 2 km from
the last or next haul-out site (grey line in a). The Direction category ‘Other’ (dark grey bands) was identified by exclusion (see
Methods). The approach used was able to identify multiple return trips to the vicinity
of haul-out sites between two haul-out periods. Data shown are for individual pv30-01-09
during a selected period of time

Finally we detected periods of resting while diving as described in Ramasco et al.
36]. We calculated the proportion of vertical ascent to descent speed and termed this
proportion dive skewness (SK). We assumed that periods of resting while diving would
be indicated by series of consecutive skewed dives and that changes in SK would occur
abruptly at the shifts between resting and other behaviours. We then used a multiple
changepoint method 52] to detect breaks in continuous series of log(SK) based on shifts in the mean. The
segments obtained were then classified by fitting a mixture of normal distributions
to the frequency distribution of mean log(SK) for each segment (see 36] for details). We then summarized the information in a categorical variable (RestingD)
indicating if the majority of the time per trajectory segment was spent resting or
in activity.

Model fitting

Trajectory segments (duration 20 min) including haul-out, surfacing or shallow diving
(5.6. m) behaviour for more than 50 % of the time were excluded from the dataset,
together with segments occurring within a RT radius distance from haul-out site or
belonging to the ‘transiting’ and ‘within range’ Direction categories (59 % of data).
A certain degree of autocorrelation was still assumed to be present in the reduced
dataset (N?=?73 629), therefore resampling of s random subsets of n data points each was used to repetitively fit the models at each stage of model selection
in order to reduce the effect of autocorrelation on parameter estimation.

Linear mixed effects models were fitted with individual as a random effect. Sex was
not accounted for in the models, since no differences were found in vertical and horizontal
indices across sexes (stBT: females CI?=?0.21 – 0.94, males CI?=?0.21 – 0.91; HS:
females CI?=?0.01 – 1.13 m/s, males CI?=?0.01 – 1.20 m/s). To render the relationship
between the vertical and the horizontal foraging indices linear and positive in case
of increase in searching intensity in the two spaces, HS and RT were transformed into
–HS and –1000/RT respectively (corr (–HS, –1000/RT)?=?0.87), where the factor 1000
was used to scale their values to similar ranges for comparison of the model parameters
[see Additional file 3]. For model selection, first the appropriate random variance structure was investigated
by comparing three full models with respectively no random variance, random intercept
() and random intercept and slope for the hFI ():

(3)

(4)

(5)

with ?X_ij being the matrix of covariates and their parameters, and e_ij the error for the ith individual and jth point. The three models, all fitted using
reduced maximum likelihood estimation (REML), were compared by likelihood ratio tests
(significant p-values??0.01, s?=?100, n?=?7000) and the best structure was considered to be the one selected most often across
the s repetitions.

Using the chosen random structure, fixed effects were then added by forward model
selection, from a minimum model including only the hFI,

(6)

up to a full model including all covariates and 2-ways interactions (fitted using
maximum likelihood estimation, ML),

(7)

The interaction between RestingD and hFI was not tested since resting dives occur
almost exclusively when the animal is stationary (small values of HS, large values
of RT and MT?=?1). The variables were sequentially added (based on the Akaike’s Information
Criterion scores) and kept in the model if a likelihood ratio test was found significant
with a p-value threshold of 0.05 (a p-value threshold of 0.01 was also used for comparison, see Fig. 2). Models were fitted 30 times (n?=?7000) and the frequency and order of selection of each covariate used to choose
the best set of fixed effects. Covariates included in the models at least 1/3 of the
times were chosen. The final model was fitted (using REML) and parameter errors estimated
by bootstrapping (100 repetitions, n?=?7000). Model validation was performed by visually assessing the presence in the
residuals of non linear patterns or violation of the assumptions of homogeneity and
normality.

To assess the influence of the temporal resolution of the data we resampled the trajectories
every p points, with p?=?3, 9, and 15, simulating decreasing temporal resolutions of respectively one point
every 1, 3, and 5 hours (20 min * p). Numerical covariates were re-estimated either by averaging the values every p trajectory segments (for the dive variables) or re-estimating the variables from
the new trajectories (for HS, RT and Direction). To avoid fitting new switching state-space
models at lower resolution, due to the high computational effort required for these
models, movement type was estimated by assigning to the new trajectory segments the
most frequent of the two states computed at the highest resolution (hence the choice
of a set of uneven p values in order to always have a majority of either state). We then performed forward
model selection and parameter estimation as previously described.

All data processing and analyses were performed in R 3.1.1 53]. State space models were run using the bsam package 54]. RT was computed using the adehabitatLT package 55]. Mixed models were fitted using the nlme package 56].