Introduction
Light drives the spatial distribution of lower-trophic organisms,
influencing their availability to predators across diel cycles
(Hamilton, 1971; Milinski, 1984; Thiebault et al. , 2016). Diel
variations in distribution patterns occurs across taxa and systems, and
are documented in terrestrial vertebrates (Slavenko et al. ,
2022), invertebrates (Chang and Hanazato, 2004; McMunn and Hernandez,
2018), amphibians (Rocha et al. , 2015) and aquatic species
(Haney, 1988; Cardinale et al. , 2003; Regular, Hedd and
Montevecchi, 2011; Mehner, 2012; Isaksson et al. , 2019). During
twilight transitions, light may still be sufficient for visually
oriented predators to exploit temporal windows of enhanced foraging
conditions while prey shifts between distribution patterns (Pyke, 1984;
Nishimura, 1992; Mella et al. , 2018). Thus, diel variations in
prey distribution modify their availability within a foraging site,
causing the cost of prey capture to vary greatly for predators
(Hamilton, 1971; Kacelnik, 1979; Kacelnik and Houston, 1984; Milinski,
1984; Wilson et al. , 1993; Benoit-bird and Au, 2003; Simset al. , 2005; Bollens et al. , 2011; Regular, Hedd and
Montevecchi, 2011; Thiebault et al. , 2016). Uncovering
predator-prey relationships under natural diel cycles is crucial for
understanding how species divide into ecological niches and adapt to
environmental constraints. These insights are essential for protecting
biodiversity, as they reveal the mechanisms that sustain ecosystem
stability and resilience.
Small pelagic species are fundamental in aquatic food webs, and often
have distinct diel distribution patterns. Many species utilize diel
vertical migration by ascending to the surface at night and descending
at dawn to spend daylight hours at deeper depths or closer to the seabed
(Cardinale et al. , 2003; Hays, 2003; Axenrot et al. , 2004;
Klevjer et al. , 2016; Kaartvedt, Christiansen and Titelman,
2023). In addition, many small pelagic species form aggregations e.g. in
schools/shoals/swarms (Delcourt and Poncin, 2012) during daylight and
disperse to feed at night (Magurran, 1990; Nilsson et al. , 2003;
Hensor et al. , 2005; Bertrand et al. , 2006). For prey,
aggregation behaviour reduce an individual’s chance of being predated
through the dilution effect (Hamilton, 1971). For predators,
aggregations may be easier to detect than dispersed prey but
intimidating and confusing upon approach, leading to reduced efficiency
and catch success rate in foraging attempts. Thus many predators have
adapted foraging strategies aiming to disrupt aggregations for isolating
more easily catchable individual prey (Lett et al. , 2014).
Interestingly, changes in spatial distribution of prey by depth and
aggregation can interact, potentially leading to great impact on
predators foraging behaviour (e.g. efficiency, effort, prey type
selection and hunting tactics), while light levels must be suitable for
prey detection. Thus for predators, identifying times when prey are both
abundant and vulnerable is crucial (Waggitt et al. , 2018), yet
how predators adapt to these fluctuating opportunities remains
underexplored due to challenges in collecting high spatiotemporal
resolution data on both predator behaviour and prey distributions
(Regular et al. , 2010; Regular, Hedd and Montevecchi, 2011). We
here address this gap using high-resolution hydroacoustic data collected
by an autonomous unmanned surface vehicle (USV), matched with the
foraging activity (monitored by animal trackers) of two sympatric
seabird species reliant on pelagic prey. By capturing the spatiotemporal
overlap between prey’s diel distributions and predator’s foraging
patterns, we test how prey distribution and light conditions influence
the predators’ foraging behaviour.
We hypothesize that predators’ foraging effort should increase when dive
efficiency (i.e. time cost/gain of dives and resting) becomes higher,
aligning with shallow depth distributions and low levels of aggregations
of prey during twilight. Using two closely related predators that
partition their foraging niche based on flight and diving cost
efficiencies, Uria aalge (common guillemot/murre) and Alca
torda (razorbill), we illustrate how dive capacity (Carbone and
Houston, 1996) and sensitivity to light levels (Smith and Clarke, 2012)
shape temporal foraging strategies and adaptations. We investigate diel
prey distribution patterns (i.e. depth distribution and number of
aggregations; hereby PDPs), to test whether differences in dive ability
between the alcids translates to differences in nyctohemeral foraging
effort. Finally, we assess how PDPs and light levels influences niche
partitioning of the predators through differing behavioural responses
(number of dives, bout length, dive depth) and dive efficiency
(bottom-dive cycle duration ratio, BDCr).
Methods
1 Study system
The study was performed in the western part of the central Baltic Sea at
and around the island of Stora Karlsö (Appendix A1, Fig. A1), a breeding
site for 26,000 pairs of U. aalge , and 12,000 pairs of A.
torda . Both species forage predominantly on Sprattus sprattus (Sprat) and Clupea harengus (herring) (Kadin et al. , 2016;
Engwall, Waldenström and Hentati-sundberg, 2022). The role of the
rapidly increasing Gasterosteus aculeatus ( three-spined
stickleback) population (Olin et al. , 2022) as food source for
these seabirds has not yet been evaluated, though it occasionally has
been recorded as chick feed (Kadin et al. , 2012) and is thus
likely to be part of adult diet.
2 Data collection and
preparation
All time data was set to Central European Time (CET), excluding daylight
saving, to ensure that noon and midnight fell at solar positions of 180°
and 0° azimuth, respectively.
2.1 Predator telemetry
2.1.1 Capture and handling
Over 10 years (2010 + 2015-2023), 54 U. aalge and 16 A.
torda were equipped with time depth recorders (TDRs) attached with
cable ties to plastic leg rings (Evans et al. , 2013; Isakssonet al. , 2019). Birds were captured on their nests, both on
natural cliffs and on an artificial breeding ledge (Hentati-sundberget al. , 2025), using a noose-pole or by hand. TDRs used were
either CEFAS G5 DST (accuracy ±1.5m), mass ≈ 3g, log interval 1-12s, or
Lotek Wireless LAT1500 (accuracy ±1.0m), mass ≈ 5g, log interval 4-12s.
Their deployment period varied, ranging from 24h to over-winter data
sampling, though only data from May-July was selected for this study.
Birds were equipped mainly during chick rearing (late June/early July),
except in 2022 when birds were captured during incubation (May/June).
For pairs tracked from 2010-2019, the breeding stage was determined at
deployment and retrieval (i.e. incubation or chick rearing). For
individuals from 2020-2022, breeding stage was assessed every day from
surveillance video, and divided into incubation, chick-rearing,
post-fledge (after chicks fledged) or non-breeding (individuals that did
not lay an egg or lost their egg/chick during the study). One nest was
too far from the camera to be identified, upon which the stage was set
to the mean of the nearby pairs.
2.1.2 Foraging behaviour and
efficiency
All dive sequences were visually checked using the R package diveMove
(Luque, 2022) and calibrated with zero offset correction adjusted for
each tracking. Dives shallower than 2m or shorter than 3s were excluded
(i.e. to remove cleaning dives and inaccurate surface values). Dives
were classified as V-shaped when horizontal movements at the bottom
lasted less than 2s. V-shaped dives constituted 18% of dives inA. torda and 14% dives in U. aalge , and were excluded
from analysis on foraging efficiency and dive depth due to their lack of
horizontal bottom movement generally considered active foraging (See
Appendix A1). Post-dive surface durations longer than 300 seconds
(Chimienti et al. , 2017) were used to distinguish ‘bouts’ (i.e. a
set of sequential dives). Foraging behaviour was summarised through 4
parameters: dive depth, bout length, number of dives per 10° azimuth
(See 2.3 Light levels ) and dive efficiency. Dive efficiency
represents the proportion of a dive cycle allocated to foragingt , based on the total dive duration T and the surface
duration after the dive S , referred to as bottom-dive cycle ratio
(BDCr), i.e. \(\frac{t}{T+S}\), according to Carbone and
Houston(Carbone and Houston, 1996). For efficiency estimations, the last
dive of each bout were excluded as it represents activities beyond the
recovery from the dive (e.g. flying to a new patch).
2.2 Prey distributions
Prey distribution was assessed from hydroacoustic data collected via a
remotely operated unmanned surface vehicle (USV), Sailbuoy (Offshore
Sensing AS, Bergen, Norway), equipped with a Simrad Wide Band
Transceiver WBT-mini echosounder with an ES-200CDK transducer from
Kongsberg Maritime (frequency sweep 185-255 kHz, see Carlsen et
al. , 2024). Data were collected from April to July in 2019-2023 (42-81
days per year), according to 10 min cycles (7min of recording – 3min
pause). This almost entirely covered the seabirds foraging area in space
(~60x80km2, see Appendix A2).
Prey depth was extracted from Echoview Software Pty Ltd v. 13. After the
estimation of calibration values, Nautical Area Scattering Coefficients
(NASC; m2 nmi−2) (Simmonds and
MacLennan, 2005) were estimated through echo-integration over the
7-minute intervals per 4m-depth layer (Appendix A3 and Carlsen et al.
2024). Prey depth was then extracted as the mean depth of the median
layer (i.e. at which the cumulative backscatter accounted for 50% of
total) to inform on prey vertical migration and change in prey
availability, based on 717,525 NASC values, divided over 36,292 water
columns summaries.
Prey aggregations (i.e. schools/shoals) were detected through image
processing of the raw echograms in Python 3.11 using echoedge which
relies on echopype (Lee et al. , 2024) 0.8.1 for reading and
pre-processing of raw data (See appendix A4). Aggregations were
extracted through thresholding methods with a minimum size fixed to 2x2m
and excluding aggregations within 1m from the seabed. As drone velocity
affected the number of aggregations detected, we (i) removed
aggregations observed at extreme speeds (<0.17 or
drone velocity as a covariate in all models. The final dataset had
40,570 observations of single aggregations, summarized into 7,089 bins
of 10° azimuth (~40min time-periods/ 1.3km with average
USV speed).
3 Statistical analyses
All statistical analyses were performed in R (version 4.2.2, R Core
Team). Some outliers were removed but this never represented
information.
3.1 Model variables
Analyses were divided into three groups: light dependent for (i)
predators and (ii) prey behaviour (hereby predator~light
and prey~light models), and (iii) prey dependent, i.e.
predator behaviour based on predictions from light dependent prey
behaviour (hereby predator~prey models).
3.1.1 Dependent variables
The predators’ behaviour response variables (i.e. efficiency-, number-
and depth of dives, and bout length) were selected based on their
ability to represent different aspects of foraging behaviour, while
being generalizable across species. Number of dives were summed per 10°
azimuth within a day/individual. Bottom-dive-cycle ratio (BDCr) informs
on the predators efficiency of presumably foraging dives (Carbone and
Houston, 1996).
For modelling prey distributions, the response variables were prey depth
and number of prey aggregations per 10° azimuth as they are likely to
affect diving predators (see 2.2 Prey distributions ).
3.1.2 Independent variables
Light variables were calculated based on the geographical position of
the USV for prey data, and the position of the breeding colony for
predators foraging data (Appendix A6). Light dependent models all
included the fixed effects ‘azimuth’, ‘zenith’, ‘week’, and interactions
between the three. In addition, the prey~light model for
aggregations included USV ‘velocity’ as a covariate to correct for the
speed of the drone (see Prey aggregations above).
Predator~light and prey~light models
included the temporal variable ‘week’, whilst
predator~light models in addition included the factor
breeding stage (‘Stage’, i.e. incubation, chick rearing, post-fledging
and non-breeding) in interaction with light variables.
The PDPs in predator~prey models were based on the ‘fish
depth’ and ‘number of aggregations’ as predicted from the
prey~light models, with interactions between the PDPs,
and between those and ‘week’. For predator~light and
~prey models where single dives/dive cycles were the unit
of observations (i.e. dive efficiency and dive depth), the number of the
dive in a bout (‘Diven’), and the total number of dives in the bout
(‘bout length’) were included as fixed effects, to reduce temporal
autocorrelation. All predator models had fixed effects for ‘sex’, and
random effects of identity of the individual (‘BirdID’), logging
interval of the TDR (‘TDRlogint’) and ‘Year’. However, we did not
include year effects to explain prey aggregation or depth, as those
models were used to predict over the entire study period (2010,
2015-2023) even when no prey data were collected. See table 1 for full
model structures (i.e. before model selection).
3.2 Distribution families
Model distribution families were gamma (link: log or identity) or
Tweedie (power 1.1-1.99), with only natural logarithmic bout length and
prey depth fitted with a Gaussian distribution. We visually inspected
residual distributions (Appendix A8) to assess assumptions and model
fit, and heteroscedasticity was inspected by plotting residuals against
model fit.
3.3 Modelling
Generalized additive mixed effect models (Wood, 2024) were used in
hypothesis testing (GAMMS; R package mgcv, see Wood, 2017). To account
for variance, correlations, and cohort effects within predator and prey
behaviours, we performed model selection to prevent overfitting while
ensuring biologically meaningful groupings in the data were accounted
for. We used backward stepwise selection from full models, based on
generalized cross-validation (GCV) and deviance explained (DE),
following the principle of parsimony and model simplicity. All
parameters were checked for correlation and variance inflation to ensure
they could appropriately be included in the same model. See full and
final models, including distribution families and smooths in Appendix
Table A7.
The light dependent behaviour models for both predator and prey allowed
for investigation into the effects of light levels independently of prey
distribution patterns (Fig 2). All effects mentioned in results and
discussion were statistically significant (p<0.05) unless
otherwise stated. Parameter effects are supported by deviance explained
(DE) and partial deviance (pd) in % figure texts. We here
focus solely on the effects of solar azimuth and PDPs, but see model
summaries in Appendix A9 and figures of interactions in Appendix A10 and
A11. When prey distributions were predicted on the predators movement
(i.e. for predator~prey models), the velocity of the USV
was set to the mean of observed velocity, 0.55m/s (See Prey
aggregations above). Otherwise, the predictions were based on the focal
values of predators dives with no modification made to the PDP models
before.