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.