Using Predicted Patterns of 3D Prey Distribution to Map King Penguin Foraging Habitat

Introduction

The Southern Ocean is one of the most pristine ecosystems on the planet. At present, 7.2% of the ocean surface is protected, and this may expand to 11.2% if proposed marine protected areas (MPAs) are sanctioned (Hindell et al., 2020). As part of the Global Ocean Alliance push for ‘30 by 30’, there is an aspiration to afford protection to more. Understanding the distribution of predators at sea is part of the process of developing spatial management plans. Hindell et al. (2020) proposed areas of ecological significance (AESs) based on predator tracks, but these AESs do not consider prey distributions directly. Determining how marine predators find their prey and hence better understanding their foraging ecology and role in ocean food webs is important for marine conservation.

Foraging habitat choices by central-place foragers during the breeding season are determined by their colony locations and dive capabilities, relative to the three-dimensional spatial distributions of prey (that may themselves show diel or seasonal variability due to migration) and competition with conspecifics and other predator species (Fauchald, 2009; Boyd et al., 2015). In turn, the self-organising behaviour of prey into, e.g., schools, swarms, patches and layers for protection from predation (Brierley and Cox, 2015; Rieucau et al., 2015) leads to prey aggregating at multiple spatial scales (10s of m for schools to 100s of km for layers). The ability of predators to track prey aggregations across these different scales depends on their search strategies (Sims et al., 2008), that typically involve a mixture of memory, social facilitation and area-restricted search (ARS) (e.g., Davoren et al., 2003; Wakefield et al., 2013). The strength of predator–prey relationships are therefore often scale-dependent, since large, persistent and predictable prey aggregations are more likely to be found than small, ephemeral patches (Fauchald, 2009). Small patches may be nested within larger aggregations, adding further complexity (Benoit-Bird et al., 2013). Spatial distribution of air-breathing predators at sea will also vary as a function of prey depth and predator diving capability, so predator distribution may also be more closely linked to accessibility of prey rather than abundance (Boyd et al., 2015). This can lead in some instances to apparent avoidance of dense prey patches that are too far from the colony (Bertram et al., 2017), or too deep in the water column (Boyd et al., 2015). Predators may also not show strong associations with prey when food is superabundant (Bertram et al., 2017), distributed uniformly (Logerwell et al., 1998), or very scarce (Smout et al., 2013; Campbell et al., 2019). Studies of predator distribution relative to that of prey therefore need to be conducted in the context of the predators’ foraging behaviour and the spatial and temporal scales at which prey aggregate (Benoit-Bird et al., 2013).

In the open ocean, fish and zooplankton form super-aggregates that are typically 10s of m in height and that can extend for 100s of km. These structures, known as deep scattering layers (DSLs; often observed between 200 and 1000 m, but have been observed outside this range) since they are readily detectable using echosounders, may offer an almost ubiquitous prey field, potentially providing deeper diving predators (e.g., elephant seals; Robinson et al., 2012) with a reliable and predictable food source. Mesopelagic fish are an important component of DSLs whose global biomass is estimated to be between 1,800 and 16,000 million tonnes (Irigoien et al., 2014; Proud et al., 2019), making them a major component of ocean ecosystems (Saunders et al., 2015). Myctophids are a small (2–20 cm) but abundant forage fish present in most oceans (Gjøsaeter and Kawaguchi, 1980; Catul et al., 2011) and form an important component of mesopelagic biomass in the Southern Ocean (Hulley, 1981; Lubimova et al., 1987). They comprise the main prey of numerous diving predators such as southern elephant seals (Mirounga leonina) and king penguins (Aptenodytes patagonicus) (e.g., Adams and Klages, 1987; Cherel et al., 2002; Bradshaw et al., 2003). However, since we lack observational data on the underwater distribution of potential prey across seascapes (conventional acoustic surveys are limited to observations along transects, but see Makris (2006) for areal acoustic observations), the times and depths at which prey become accessible to predators remains unclear.

Tracking seabirds in the horizontal dimension enables exploration of drivers of distribution in the context of widely available physical variability (from remote sensing). Typically, tagged seabirds are tracked from their colonies or counted from vessels, and their density is either related to the distribution of prey (from vessel-based acoustic/net/environmental sampling, e.g., Logerwell et al., 1998; Hentati-Sundberg et al., 2018; Campbell et al., 2019), or associated with remotely sensed oceanographic variables that influence prey densities like sea surface temperature (SST), bathymetry, chlorophyll concentration and the locations of oceanographic fronts (e.g., Bost et al., 2011; Scheffer et al., 2012; Scales et al., 2014; but see Grémillet et al., 2008; Sherley et al., 2017). As the climate changes, foraging habitats are expected to shift geographically in line with these associations, and this could have consequences on marine predator populations (Sherley et al., 2017; Cristofari et al., 2018).

King penguins are distributed widely across the Southern Ocean and breed on subantarctic islands such as Kerguelen, Crozet, Macquarie and South Georgia (Woehler, 1993; Ratcliffe and Trathan, 2011). These islands, and much of the penguin foraging grounds, are protected by MPAs, e.g., the South Georgia and the South Sandwich Islands (SGSSI) MPA (2012). The diet of king penguins is predominantly comprised of myctophids, particularly Krefftichthys anderssoni at South Georgia, whose predominance in the diet is associated with higher king penguin breeding success (Perissinotto and McQuaid, 1992; Olsson and North, 1997). King penguin foraging distributions have been associated with warm-core eddies and strong temperature gradients in the Antarctic Polar Frontal Zone (Scheffer et al., 2010). For example, the effects of a one degree latitude shift in SST anomalies on the location of the Polar Front could extend the travel distance of king penguins breeding at Crozet Islands by 130 km, which is a substantial proportion of their expected foraging range (typically between 300 and 600 km) (Bost et al., 2015). But since king penguins are central place foragers, that have fixed breeding localities (islands) and a scarcity of alternative breeding sites, their ability to change foraging distribution in response to spatial shifts in prey distributions may be limited (Bost et al., 2009; Sherley et al., 2014). Constraints on foraging include distance that an individual can travel from its colony that may vary with stage of the breeding season (Scheffer et al., 2016), diving capability (<400 m; Charrassin et al., 2002) and the fact that king penguins are visual foragers, which means foraging is predominantly carried out in daylight (Bost et al., 2002); this limits their ability to feed on myctophids within DSLs that migrate to the surface at night (Hazen and Johnston, 2010). These constraints lead us to expect that king penguins would preferentially target fine-scale (10s of m) dense fish schools rather than larger scale prey structures (10s–100s of km) such as DSLs that may be inaccessible most of the time, conditional on the schools being at relatively shallow depth during the day and within the stage-specific foraging range of the colony.

In this study, we attempt to relate the distribution of king penguins during guard stage tracked with GPS, time-depth recorders (TDRs) and accelerometers from South Georgia to that of their prey, inferred from diet studies and acoustic surveys conducted during the Antarctic Circumnavigation Expedition (ACE, December 2016 to March 2017). Despite our best intentions, we were not able to obtain coincident observations of predators and prey, simply because the penguins did not forage close to the vessel: hence we built models to predict prey distribution at penguin dive locations. The specific objectives of this study were to (i) compare fine-scale (10s of m) vertical distributions of fish schools and layers (determined from low frequency echosounder observations) and king penguin prey encounters (analysis of accelerometery data) observed in the same interfrontal zone (located between the Antarctic Polar Front and the Southern Antarctic Circumpolar Current Front) and during the same period (February and March 2017) off South Georgia; (ii) combine in situ echosounder observations recorded during ACE with environmental data to build predictive models of local-scale (10s of km) fish echo intensity and vertical distribution for the Southern Ocean, and (iii) model the horizontal habitat suitability of king penguins based on predicted prey distributions to answer the question: what characteristics of prey do penguins prefer to travel to? Note that this is different to modelling prey encounter dives versus non prey encounter dives which asks a different question (i.e., what characteristics of the prey are associated with penguins finding prey).

Materials and Methods

Data were collected as part of the ACE (20th December 2016 to 18th March 2017, Figure 1) aboard the research vessel ‘Akademik Tryoshnikov’. Echosounder data were collected continuously throughout the voyage and king penguin foraging behaviour data were collected simultaneously off South Georgia using data loggers. Unfortunately, the ACE vessel was not equipped to carry out net tows and therefore no biological sampling was conducted. However, DNA was extracted from penguin faecal samples collected during the study period to determine penguin diet and provide information on prey consumed.

King Penguins

King penguins breeding at Hound Bay (54°23’S, 36°15’W), South Georgia, were studied between 20th February and the 22nd March 2017. This period coincided with the passage of the ACE voyage past South Georgia (2nd to the 8th March 2017). Only adult chick-rearing king penguins in guard stage were included in this study.

Animal Handling and Instrumentation

A total of 18 individual king penguins were captured as they were leaving on foraging trips. They were instrumented with miniature data-loggers attached to their lower backs with waterproof Tesa^® tape and quick-drying two-part epoxy, and marked on their front with permanent black hair dye to aid in reidentification and device recovery (Scheffer et al., 2010). At-sea movements were followed using GPS (Pathtrack nanofix GEO, weight: 32 g; size: 54 × 23 × 17 mm), and dive profiles were recorded using TDRs (G5, CTL, Lowestoft, United Kingdom; weight: 6.5 g; size: 12 × 35.5 mm) recording at 1 Hz. Prey strikes were measured using a tri-axial accelerometer (X8M-3, Gulf Coast Data Concepts, Waveland, MS, United States, weight: 34 g; size: 51 × 25 × 16 mm) recording surge, heave and sway at 25 Hz. Birds were recaptured after one foraging trip and the loggers were retrieved. Battery capacity limited the observation period for both TDR and accelerometer tags. Typically, the accelerometer tags ran out of battery 4–5 days before the penguins returned from foraging at sea.

Foraging Behaviour

Trip characteristics were extracted from GPS data using the ‘trip’ package in R (Sumner et al., 2009) following the application of a speed filter of 14 km h^–1 (e.g., Cotté et al., 2007) (which is the fastest observed swim speed of a penguin, see Kooyman and Davis, 1987) to remove any erroneous position estimates (<0.1% of all positions). For each bird, trip duration (days), path length (km), maximum distance from colony (km), and mean and maximum speed (ms^–1) were calculated. TDR data were processed using IGOR Pro (WaveMetrics, Version 7, OR, United States) to calculate maximum dive depth; dive duration; number of undulations (also called ‘wiggles’, defined as a depth change rate < 0.25 ms^–1); bottom phase duration; and post dive duration (see Ropert-Coudert et al., 2007, for parameter definitions). Depth undulations (wiggles) are indicative of either prey search, prey encounter, predator interaction or other movement, and hence are not a reliable source for detecting prey capture attempts. Shallow travelling dives ≤4 m were excluded (Charrassin and Bost, 2001; Charrassin et al., 2002). For five individuals, the GPS battery expired before completion of the foraging trip so the location of subsequent dives could not be determined.

Prey Encounters

Prey capture attempts were identified using the accelerometry data (e.g., Ropert-Coudert et al., 2006; Zimmer et al., 2011). Vectorial Dynamic Body Acceleration (VDBA), a commonly used proxy of whole body activity (Gleiss et al., 2011), was calculated as:

where A_x, A_y and A_z are the dynamic accelerations of the surge, heave and sway axes respectively, calculated by subtracting the static acceleration of each axis (a smooth over 1 s) from the total acceleration (Shepard et al., 2008). There was an upper inflection point at 2.5 m s^–2 in the frequency of VDBA across all birds, and this was used as a threshold for diagnosis of prey capture attempts (Supplementary Figure 1), i.e., when VDBA (at 1 Hz) became higher than 2.5 m s^–2 during a dive, we inferred that there was a prey encounter and pursuit (Sánchez et al., 2018). This meant that multiple prey encounters (e.g., multiple captures of prey) might be detected as a single prey capture attempt. This was deemed acceptable since our analysis is reliant on the identification of prey encounter dives and prey encounter depths, not the absolute number of prey items (which would require additional sensors, e.g., cameras). For each foraging dive (dives > 4 m) with a prey capture attempt, the depth at which the first VDBA point exceeded 2.5 ms^–2 was considered as the prey encounter depth (m). The first 4 m of the descent phase were excluded from the acceleration analysis to be consistent with the dive data and because birds perform strong wing beats at the beginning of a dive to overcome the effects of buoyancy and inertia (Zimmer et al., 2011).

Diet Study

During the same period as the biologging study, 47 faecal samples were collected from other (i.e., not the instrumented birds) king penguins (Le Guen et al., 2020). Ideally, we would have collected samples from the instrumented birds, but due to limited resources and sampling logistics, we were only able to collect samples from birds at the colony whilst the instrumented birds were at sea. Since all the instrumented birds were rearing chicks, only samples taken from chick-rearing birds were considered in this study (16 out of 47 samples collected between 22nd February and 11th March 2017). It is likely that samples were taken from different individuals as they were collected from different nests across the colony (c. 100–150 penguins in total), but it is possible that faeces from the same bird were sampled more than once since adults swap nests regularly. Contamination was minimised during collection by maintaining strict sampling protocols, e.g., sterilising apparatus with ethanol between sampling and only taking samples from the surface of the faecal matter, i.e., not making contact with the underlying soil (for more information, see Le Guen et al., 2020). These data provide vital information on colony-specific diets, which could potentially change seasonally and/or annually in response to changes in prey distributions (e.g., Saunders et al., 2013). Prey taxa for each bird were identified using DNA sequencing (for DNA extraction and sequencing protocols see Supplementary Material).

Echosounder Observations

Echosounder observations were collected continuously along the ACE track (Figure 1) using a Simrad (Bergen, Norway) EK80 echosounder operating a single-beam (10–20 kHz) Wärtsilä ELAC Nautik GmbH (Keil, Germany) LSE179 transducer. Prior to the survey, at-sea trials were conducted to find the optimum settings (i.e., high signal-to-noise ratio) for detection of mesopelagic fish communities: frequency = 12.5 kHz, ping interval = 8 s, power = 150 W, and pulse duration = 16 ms. At this frequency (wavelength c. 12 cm), the majority of the returned signal is expected to originate from large organisms (fish, squid, etc., not zooplankton; Simmonds and MacLennan, 2005), and organisms that possess gas bladders (Proud et al., 2019), e.g., K. anderssoni. An echosounder calibration was conducted at South Georgia using the standard target method (Demer et al., 2015).

Echosounder Data Processing

Raw echosounder data were imported into Echoview version 8 (Myriax, Hobart, Tasmania) and the calibration parameters were applied. Transient noise, impulse spikes and attenuated pings were removed. Data were gridded at a resolution of 5 m depth by 10 km along track, and integrated values (Mean Volume Backscattering Strength, MVBS, dB re 1m^–1) were exported to provide data for modelling. Data were also integrated and exported at fine-scale (i.e., 1 ping × 1 m sample) to facilitate comparisons with fine-scale king penguin prey encounters. MVBS values were used to calculate the nautical area scattering coefficient (NASC, average echo intensity per nautical mile squared; m² nmi^–2) and weighted mean depth (WMD, mean depth weighted by linear values of MVBS) values for two depth zones: (1) the penguin foraging zone, between the surface and the depth at which 95% of prey encounters occurred, and (2) the mesopelagic zone, this – by common definition – extended between 200 and 1000 m, and was expected to contain king penguin prey (e.g., myctophids), some of which would have been deeper than the dive range of penguins. At the 12.5 kHz echosounder frequency used, the local WMD value for the mesopelagic zone is essentially equivalent to the mean depth of the principal (strongest) DSL, since during ACE a single broad DSL was observed through most of the expedition. Following the same logic, NASC values calculated between 200 and 1000 m were effectively a proxy for DSL NASC. NASC and WMD values calculated between 200 and 1000 m will be referred to as DSL NASC and DSL depth hereafter.

Predicting Nautical Area Scattering Coefficient and Weighted Mean Depth Using Neural Networks

To model prey availability within the foraging domain of the king penguins at South Georgia (Figure 2), neural networks were built to predict NASC and WMD along the entirety of the ACE voyage track (see Figure 1). Using the entire expedition (as opposed to using a subset of the data close to South Georgia) increased the size of the training set by a factor of 25 and enabled daily and weekly variability (recorded away from the South Georgia foraging site) in prey characteristics to be captured. Environmental variables used as inputs for the neural networks (Table 1) were selected based on previous studies and known relationships between NASC, WMD and the environment (e.g., Klevjer et al., 2016; Aksnes et al., 2017; Proud et al., 2017). To explain diel variation in the depth distribution of biomass resulting from diel vertical migration, sun angle was included as an input variable (see Table 1). All input variables were at the same or higher spatial resolution as the summarised track echosounder observations (10 km), except for sea-level anomaly, which was gridded at a resolution of 0.25° × 0.25°. Values of the input variables were assigned positions along the vessel track using linear interpolation.

Each 10 km position along the ACE cruise track was classified as belonging to one of four interfrontal zones (Figure 1) evident in vertical potential temperature profiles. Daily mean potential temperature data at 25 m depth intervals from 100 m to 300 m depths were extracted from the Global coupled FOAM quarter degree model run by the United Kingdom Met office¹ along the ACE track for each day of the cruise. Each ship location was then associated with a interfrontal zone based on the potential temperatures at 200, 300, and 500 m (Figure 1) evident in previous work (e.g., Orsi et al., 1995; Belkin and Gordon, 1996; Boehme et al., 2008). All locations associated with a water depth of less than 1000 m were classified as being on the shelf.

Neural network inputs (Table 1) and outputs (NASC and WMD calculated over the penguin foraging zone and mesopelagic zone) were scaled to between 0.0001 and 0.9999 (weighting the contribution of each input equally). Temporal inputs were predominantly at daily resolution apart from PP (see Table 1), which was monthly and hence values for February and March were assigned to penguin dive locations over the study period. R-packages ‘caret’ (Kuhn, 2008, 2020) and ‘nnet’ (Venables and Ripley, 2002) were used to train and test the neural networks (10-fold cross validation; 80%/20% split for train/test data, respectively). A standard 3-layer setup (input, hidden layer, output) was used. The number of nodes in the hidden layer (6–30 nodes) and decay parameter (regularization parameter to avoid over-fitting: 0.001–0.1) varied, and optimum models were selected based on validation set mean root-mean-square error (RMSE) values (mean over the 10 folds).

Penguin Habitat Selection Modelling

We did not have direct echosounder observations at penguin dive locations, so took a modelling approach to infer the likely prey-field characteristics at dive sites on specific days during the study period. This was made possible by including ACE day (see Table 1) in the NASC and WMD models to capture daily and weekly variability; note that we used ACE day rather than day of year to avoid forcing the model to weight day 1 (January 1st) and day 365 (December 31st) very differently, when in fact these 2 days should be weighted similarly since they are only 1 day apart. Daytime king penguin dive data that could be associated with valid GPS locations were used to study king penguin habitat use. Bird 14 was excluded from this analysis due to anomalous behaviour (see trip with far eastern loop in Figure 1), which was not consistent with other studies for birds in this stage (Scheffer et al., 2012), and also because its TDR logger ran out of battery (note absence of eastern loop in Figure 2). Habitat selection models (e.g., Aarts et al., 2008) were used to define the pelagic foraging habitat suitability of the king penguins. The outer limit of available habitat was set conservatively (to limit spatial extent of sampling domain, see Hazen et al., 2021) as a circle with a radius of 1.1× the maximum foraging range from the Hound Bay king penguin colony (Wakefield et al., 2017). Dives where the seabed was <600 m (the expected maximum depth of DSLs in this region) were excluded from this analysis and hence habitat suitability in this context refers to pelagic dives as opposed to shallow inshore/benthic dives (where we did not have acoustic observations). Dive locations were assigned values of the environmental variables (see Table 1) at the same scale as the modelled NASC and WMD observations. Environmental variables were averaged over an area demarked by a circle (diameter = 10 km) centred on each dive location. Predicted NASC and WMD, distance to colony and seabed depth were extracted for each dive (presence) location and day (Figure 2) as well as for three randomly chosen pseudo-absence locations (three points drawn from the entire accessible geographic area) for every dive location to represent habitat availability (Aarts et al., 2008). Binomial generalised additive models (GAM; logit link function) with presence/pseudo-absence location as the response, and the prey indices (predicted NASC and WMD), seabed depth and distance to colony as explanatory variables were fitted using the R package ‘mgcv’ (Wood, 2006). Models were fitted using interactions between related covariates (using tensor smooths; Wood, 2006), e.g., an interaction between seabed depth and distance to colony is expected to influence habitat suitability since birds only utilise shallow water when close to the colony; note that some degree of correlation is expected between seabed and distance to colony. GAMs were parameterised using full tensor product smooth terms: (te1) distance to colony * seabed depth; (te2) NASC_{50–260 m} * WMD_{50–260 m}, representing the penguin ‘prey’ component, and (te3) NASC_{200–1000 m} * WMD_{200–1000 m} representing the ‘DSL’ component (i.e., addressing the question does the depth and intensity of DSLs influence diving behaviour?).

For each model, the random pseudo-absence selection was repeated 30 times, resulting in 30 model iterations. Model selection was conducted using Akaike Information Criterion (AIC), cross validation and Area Under the Receiver Operating Curves (AUC), which are conservative and permit violation of model assumptions relating to non-independence of errors that are inherent in tracking data (spatial and temporal autocorrelation). Three candidate models were compared using AUC: (1) the null model (te1), (2) the ‘prey’ model (te1+te2), and (3) the ‘prey+DSL’ model (te1+te2+te3). Selected GAMs were then averaged over the 30 iterations using the ‘MuMIn’ R package (Barton, 2020) and used to predict habitat use probabilities around the Hound Bay colony. Since predictions of NASC and WMD will have an associated error (RMSE from neural networks), uncertainty associated with the habitat model will be biased low.

Results

In total, 18 birds were tracked, and faecal samples were collected from 16 non-instrumented birds. Bird 14 was excluded from dive analysis due to anomalous foraging behaviour that was not consistent with previous studies, and no accelerometer data were collected from bird 16 due to a tag malfunction.

King Penguin Foraging Trip Characteristics

All king penguin foraging trips headed north from South Georgia in the direction of, but never reaching, the Antarctic Polar Front. The birds remained in the APF-SACC interfrontal zone during guard stage (consistent with other studies; see Scheffer et al., 2012). Foraging trips took place between 20th February and 22nd March and had an average duration of 10.2 days (range = 4.8 to 17.8 days; 183.8 days total), path length of 820.7 km (range = 293.9 to 1837.2 km) and mean maximum distance from colony of 317.2 km (101.7 to 672.8 km). Specific details of individual trips are given in Supplementary Table 1.

A total of 37,275 dives deeper than 4 m were identified during 165.3 days of total tracking time across all birds between February 20th and March 20th (using TDR tags; see Figure 2 and Supplementary Table 2). On average, birds dove 228 times per day and each dive included 8.2 undulations/wiggles (i.e., potential prey encounters). The maximum dive depth recorded was 352 m and the maximum dive duration was 9 min 10 s. Birds dived mostly during the day (c. 73% by dive duration), with depths increasing progressively after sunrise and reduced around sunset (Supplementary Figure 2) in keeping with both diel vertical migration (DVM) of prey and increased visibility with depth as sun angle increases (Bost et al., 2002; Pütz and Cherel, 2005).

King Penguin Prey Encounters and Pelagic Vertical Habitat Use

A total of 45,958 (43,063 with valid GPS locations) prey encounters deeper than 4 m were detected during 99.3 days of total tracking time across all birds between February 20th and March 16th (using accelerometer tags; see Figure 2 and Supplementary Table 3). On average there were 462.8 prey encounters per day and hence around 2 prey encounters per dive (i.e., 462.8 prey encounters per day/228 dives per day). The average prey-encounter depth was 107.1 m (lower quartile = 67.1 m; upper quartile = 174.8 m), and 95% of encounters occurred at a depth ≤260 m, which was used to define the lower boundary of the king penguin foraging zone (see sections “Predicting Fish Distributions Using Neural Networks”, “King Penguin Habitat”, and “King Penguin Habitat Suitability”): very few prey encounters occurred deeper, so including acoustic backscatter from beyond this depth as ‘prey’ would have resulted in inclusion of features (such as DSLs occurring deeper than 260 m) that were not typically accessed by king penguins. The vast majority of the birds’ foraging effort was off shelf (seabed > 1000 m; 38,394 out of 43,063 prey encounters were offshore). We thus compared the vertical distribution of prey encounters offshore to offshore echosounder observations (seabed > 1000 m; c. 465 km of vessel track) made in the same interfrontal zone (APF-SACC; see Figure 1) between 4th March and the 18th March 2017. Echosounder observations recorded in the top 50 m were excluded due to surface noise. We compared day (29,877 prey encounters; see Figure 3A) and night (8,517 prey encounters; see Figure 3B) separately to avoid complications associated with prey DVM. MVBS was plotted for each depth interval (0–600 m by 10 m intervals) and MVBS level (−80 to −40 dB re 1m^–1 by 3 dB re 1m^–1 intervals) e.g., the cell in the bottom left-hand corner of Figures 3A,D is the MVBS value calculated by first thresholding the data by the MVBS level (−80 to −77 dB re 1m^–1) and then integrating the remaining observed echo intensity over the depth interval (590 to 600 m).

Acoustically Detected Prey Field

DSLs, which were ubiquitous throughout the APF-SACC interfrontal zone, were the source of most of the recorded echo energy (78.3% during the day and 92.4% during the night; note that the top 50 m was excluded) but even at night these DSLs appear to be deeper than the dive capability of the king penguins (Figure 3). Fish schools (isolated echogram features, 10s of m in width and height) occurred ephemerally in the epipelagic zone, observed at both low echo intensity (prey type 1: MVBS level = −70 to −65 dB re 1m^–1; see Figure 3A) and high echo intensity (prey type 2: MVBS level = −50 to −40 dB re 1m^–1; see Figure 3A). The differences in MVBS level imply that prey types 1 and 2 differ in numerical fish density, target strength or both. A third prey type, which we termed ‘prey patches’ (Béhagle et al., 2017) (prey type 3: MVBS level = −70 to −60 dB re 1m^–1; see Figure 3A), comprised clusters of linked schools (see Supplementary Figure 3) that, due to echosounder beam geometry, would appear as layers if they were observed deeper in the water-column (the echosounder beam broadens with depth and hence sampled volume increases with depth, reducing horizontal resolution and blurring the boundaries between discrete prey structures giving them a continuous appearance). The vertical distributions of prey types 2 and 3 shifted markedly between day and night (see Figures 3A,D) and hence were likely composed of species that performed DVM. The modal depth in daytime penguin prey encounters matched with the modal depth of prey type 2 (high energy fish schools) NASC (Figure 3). Prey type 1 was persistent throughout the daily cycle and matched with peaks in penguin prey encounters both during the day and night (Figure 3).

Diet Analysis

Only 9 of the 16 chick-rearing penguin faecal samples returned useful DNA data (likely due to degradation of the DNA in the gut during digestion). Of these, eight contained DNA of eight species of fish (Table 2) but only two contained euphausiid crustacean DNA (Euphausia superba and Thysanoessa macrura). Fish also eat krill and so it is possible that presence of krill DNA in penguin faecal samples is a consequence of secondary ingestion. Higher occurrences of krill were found in the other breeding stages (3 out of 16 samples for incubating birds, and 5 out of 16 samples for non-breeding birds). No prey species were encountered in the procedural blanks, thus it is likely that no contamination occurred during the DNA extraction phase.

Consideration of the known DVM behaviour, biology, and depth distribution of these fish species (Table 2) links them to the three prey types observed (see section “Acoustically Detected Prey Field” and Figure 3). Both adult icefish (Channichthyidae spp.) and painted notothen (Lepidonotothen larseni) are associated with the seabed so probably do not make up a large percentage of the penguins’ diets (90% of prey encounters are pelagic; these fish may be taken opportunistically at the start or end of trips). However, larvae and juveniles of these species occur within the epipelagic zone and since they are without gas-filled swimbladders are relatively weak scatterers, which matches with prey type 1 (icefish also do not perform DVM, which is consistent with the observed DVM pattern). The other two prey types (2 and 3) likely correspond to rhombic lanternfish (K. anderssoni), since they are known to perform DVM, have gas-filled swimbladders (i.e., high intensity acoustic targets) and occur in both the epipelagic and mesopelagic zone.

Predicting Fish Distributions Using Neural Networks

NASC and WMD were calculated for the penguin foraging zone (50–260 m; defined in Figure 3) and the mesopelagic zone (200–1000 m), resulting in four model outputs (NASC and WMD for both depth zones) and hence four neural networks. Models were trained using 80% of the ACE echosounder data (n = 1,636 samples), and model selection and evaluation was based on validation set (n = 164 × 10 folds) and test set (n = 409 samples) RMSE, respectively. Optimum networks for each of the outputs were as follows (see Figure 4): (A) NASC_{50–260 m}, nodes = 22, decay = 0.001 (test RMSE = 295.6 m² nmi^–2; R² = 0.66); (B) NASC_{200–1000 m}, nodes = 30, decay = 0.001 (test RMSE = 289.9 m² nmi^–2; R² = 0.81); (C) WMD_{50–260 m}, nodes = 28, decay = 0.001 (test RMSE = 18.4 m; R² = 0.73), and (D) WMD_{200–1000 m}, nodes = 28, decay = 0.001 (test RMSE = 47.8 m; R² = 0.75).

King Penguin Habitat

Predicted gridded NASC and WMD daily values (sun angle was set to the local daily maximum value) were averaged over the king penguin foraging period and plotted (Figure 5; standard deviation plotted in Supplementary Figure 4). Mean areal values of NASC_{50–260 m} and NASC_{200–1000 m} over the penguin foraging domain were c. 322 and 335 m² nmi^–2, respectively. Mean areal values of WMD_{50–260 m} and WMD_{200–1000 m} over the penguin foraging domain were c. 78 and 336 m, respectively. King penguin dive locations were associated with shallow prey (Figure 5C) that had high NASC values (Figure 5A) in their foraging zone (top 260 m), and with relatively shallow (Figure 5D) and high intensity DSLs (Figure 5B) in the 200–1000 m depth range. NASC and WMD were also predicted for the region surrounding the Kerguelen Islands (see Supplementary Figure 5), which is also home to colonies of breeding king penguins. Spatial patterns of prey were very different between the two regions: in the Kerguelen region, mean NASC_{200–1000 m} was much higher (c. 1599 m² nmi^–2), mean NASC_{50–260 m} was much lower (c. 132 m² nmi^–2), and DSLs were substantially deeper (c. 487 m; a depth far beyond the diving capabilities of king penguins).

King Penguin Habitat Suitability

Excluding Bird 14 (see section “Materials and Methods”) left a total of 19,504 dives (out of 23,646 daytime dives and 37,254 in total) available for the habitat selection models. The ‘prey+DSL’ model had significantly better support than the null and prey models (Figure 6C) and a substantially lower AIC value (mean over the 30 model runs was 18,472 compared with mean values of 21,226 and 30,717 for the prey and null models, respectively, indicating that the more complex prey-based models were not being overfit). The ‘prey+DSL’ model was used to predict mean habitat suitability over the period (Figure 6D); note that NASC values are integrated over the penguin foraging zone and the mesopelagic zone in these models and hence are much larger than the values in Figure 3, which are integrated over 10 m depth intervals.

King penguin habitat suitability declined with distance from the colony, and since we used the distance to colony * seabed depth interaction (te1), shallow seabed depths tended to be used close to the colony but were avoided as distance from it increased. Habitat suitability increased in the presence of high intensity and relatively shallow DSLs (Figure 6B). Three distinct features were evident in the partial plot for the prey component of the GAM (te2; Figure 6A) that were related to the three prey types defined in Figure 3 (based on WMD and NASC) and Table 2 (based on diet studies).

Discussion

Our results indicate that king penguins in guard stage at South Georgia do not feed directly on the ubiquitous DSLs around their colony, as these were found to be outside their diving range during the study period. Instead, they likely target smaller, ephemeral schools occurring above the DSLs. This offers strong support for the idea that availability and prey type (or quality), rather than abundance per se influences predator–prey interactions at fine (10s of m) vertical scales (e.g., Benoit-Bird et al., 2013; Boyd et al., 2015). The link between targeted fish schools and the inaccessible DSLs, which is evident at a 10 km horizontal scale, may be related to ontogeny (e.g., juveniles forming clusters of shallower aggregates) or other food web links (e.g., changes in depth distribution related to species, size, DVM, or environmental conditions, e.g., light). However, these findings are based on a single site, breeding stage and year and warrant further investigation.

Study Caveats, Considerations, and Recommendations for Future Work

The approach taken here provides the foundations of a framework to tackle an important problem in marine ecology – linking at-sea diving predator distributions with their prey across multiple spatial and temporal scales – using a combination of acoustics, diet studies, animal tracking and modelling. Such a comprehensive and multi-faceted approach naturally incurs limitations, which will ultimately need to be addressed in future studies. For example, there was poor spatial overlap between penguin dive locations and the vessel track. Use of a single echosounder frequency did not enable target identification. However, since the frequency used was relatively low (12.5 kHz), most of the received backscattering intensity could be attributed to fish. Due to surface noise, we were not able to make observations in the top 50 m and hence some schools at the surface may have been missed. This may explain why the shallow zone neural networks systematically underestimated NASC and WMD (see Figure 4). In future studies, structured and statistically robust (see Jolly and Hampton, 1990) multifrequency/broadband acoustic surveys with coincident trawling in known (or predicted) penguin foraging grounds should be conducted to enable prey species identification and provide information on prey size distribution.

We compared penguin prey encounters with predicted NASC values to study predator–prey interactions. However, not all predicted encounters will be a result of prey interactions since prey and predator encounter discrimination is not possible using accelerometery data alone, but the number of predation events is expected to be relatively low. We also used a neural network to predict prey characteristics because relationships between acoustic target strength and the environment are known to be non-linear and complex (e.g., see Jech et al., 2015). Neural networks provided very good predictive power but are hard to interpret, precluding any discussion related to the importance of the correlated model inputs.

Spatially resolved ecosystem based management requires information on predator and prey distributions and on the interactions between them. We used model-predicted physical and biological variables to model prey NASC and WMD, which were then used to model predator habitat suitability. Instead of this multi-model approach, a single model based on the environmental covariates could have been used to predict habitat suitability. Whilst this less complex approach would ultimately be more statistically robust, it would not provide the means to determine the importance of the vertical distribution of different prey structures (e.g., scattering layers and schools) in predicting habitat suitability (i.e., providing information on predator–prey interactions). In the absence of coincident prey observations, future modelling work should test the utility of our multi-model approach by comparing predictions with output from simpler models that only consider environmental observations.

We also modelled (binomial GLM) prey encounter dives (=1) versus non-prey encounter dives (=0) using the habitat suitability model predictors (analysis not presented) but did not find good agreement between the data and the models. This suggests that the data we used (predicted NASC and WMD values) do not explain the variance in prey capture well, which is likely to be better explained by environmental covariates that contribute to successful prey capture, e.g., surface illumination and light attenuation in the water-column (not available in this study).

King Penguin Prey

King penguins from South Georgia likely fed on barracudina close to the seabed and juvenile icefish and painted noties, that typically occupy reefs and kelp beds in the coastal region, on their return trips to the colony (Frolkina et al., 1998). Although relatively infrequent, occurring in just 13% of faecal samples that contained fish DNA from individual chick-rearing penguins, barracudina may make up a relatively large proportion of king penguin diet by mass due to their large size (Olsson and North, 1997). In offshore regions, where the majority of foraging takes place, penguins targeted specific prey types (isolated and larger-scale patches of fish schools, comprised of various species across different life-history stages and depth ranges) within their diving range. Given that king penguins typically require c. 2.3 kg/day of fish (Bost et al., 1997; Putz et al., 1998), our results indicate that penguins must have on average obtained around 10.1 g of prey per dive or 5 g per prey encounter. Since adult K. anderssoni (main prey item of chick rearing penguins) weigh around 4 g at maximum standard length (c. 7 cm) (see Saunders et al., 2020), this would suggest that on average each prey encounter included capture of >1 fish. Penguins likely focused on juvenile (since they do not eat larvae, see Ratcliffe and Trathan, 2011) fish species during the daytime that occurred in the mid-to-upper epipelagic zone (e.g., painted noties/icefish), but also probably fed on adult K. anderssoni both during the day and at dawn/dusk during fish vertical migration. Increases in shallow epipelagic prey density were related to relatively shallow and dense (high echo energy) DSLs, which is in agreement with the findings of Fielding et al. (2012) that acoustic backscattering intensity was positively related to the number of fish schools. This association can be explained by the strong coupling, via DVM and ontogenesis, between highly dynamic shallow water biomass (in the form of ephemeral schools and larvae/juveniles) and mesopelagic biomass (Pinti et al., 2021).

Electrona carlsbergi and K. anderssoni are known to be the dominant schooling myctophid species close to South Georgia, and occur at depth of between 60 and 300 m (Saunders et al., 2013). King penguin chick rearing success has been found to increase when K. anderssoni form a large component of their diet (Olsson and North, 1997). We found that 88% of the chick-rearing birds had K. anderssoni in their diet, whereas no E. carlsbergi were found. Saunders et al. (2013) reported that E. carlsbergi, which are typically located in water associated with the Antarctic Polar Front (Hulley, 1981; Zasel’sliy et al., 1985; Saunders et al., 2014), were absent or in low abundance in some years, whereas K. anderssoni occurred more regularly every year. Since king penguins in guard stage do not feed at the Antarctic Polar Front, the absence of E. carlsbergi in the diets of birds in this study is likely to be a consequence of the constraints on foraging ranges during the guard stage of the breeding season (Olsson and North, 1997), but it is also possible that E. carlsbergi were consumed but not detected due to the relatively low sample size in our diet study.

Layers and associated patches of swimbladdered fish have been observed during the daytime in the Scotia Sea (Saunders et al., 2015; Lourenço et al., 2017) and at Kerguelen (Béhagle et al., 2017). Prey patches at Kerguelen, which were observed shallower than 180 m, were likely comprised of K. anderssoni and Protomyctophum spp. and disappeared at night owing to DVM (as observed here). The depth range of K. anderssoni around Kerguelen is between 100 and 150 m (Bost et al., 2002; Charrassin et al., 2004) and overlaps with the depth range of prey type 2 identified in this study, classified as K. anderssoni through diet analysis. This prey type may also include Protomyctophum spp., since Protomyctophum choriodon have been commonly found in king penguin diets, particularly during February and March (Olsson and North, 1997). Therefore, foraging behaviour of king penguins at South Georgia appears similar to that found at Kerguelen as is the apparent association between targeted prey patches and DSLs.

Predator–Prey Interactions

We found good agreement between the vertical distribution of acoustically detected fish schools and that of king penguin prey encounters at fine vertical scales (10s of m), and were able to map king penguin foraging habitat based on predicted 3D prey distributions at local horizontal scales (10 km) after controlling for regional-scale (10s–100s km) effects related to distance to colony and bathymetry. Within the habitat available to king penguins, defined to be 1.1× maximum foraging range from their colony, prey resources (full water-column NASC) were found to be relatively low when compared to other similar sites (e.g., Kerguelen, see Supplementary Figure 5). King penguins concentrated their diving effort in an area with relatively high levels of predicted NASC toward (but not reaching) the Antarctic Polar Front, likely containing relatively high densities of prey. Vlietstra (2005) found that in a similar resource-limited environment, fine-scale overlap of seabird and prey distributions increased when compared to resource-rich environments. Since prey distributions can shift significantly both seasonally and annually, e.g., presence/absence of E. carlsbergi around South Georgia (Saunders et al., 2013), potential king penguin foraging habitat may change from being a resource-limited to a resource-rich environment. However, previous work has shown that king penguins have foraged in a similar area (Scheffer et al., 2010), suggesting that it may be characteriesd by predictability of resources and low variation in prey availability. It follows that the more constrained predators are by their environment (i.e., the harder it is to find prey compared with random search) the more tightly coupled the relationship between prey and predators will be. Perhaps the most important foraging constraint for king penguins is vertical accessibility, and studies of seabirds in the Humboldt Current have found that this is more important in predicting seabird foraging habitat selection than prey density throughout the water column (Boyd et al., 2015). It is perhaps no coincidence that DSL depth coincides with the maximum dive depth of king penguins: myctophids may go to these depths to reduce predation risk, as opposed to tracking isolumes (Røstad et al., 2016) or being physiologically constrained by swimbladder function (Bone et al., 1995). A relatively small shallowing of DSLs (10s of m), which is predicted by 2100 (Proud et al., 2017), could open up a very large resource to king penguins that is perhaps more reliable than the presently targeted ephemeral fish schools close to the surface. This would of course reduce prey variance at fine and local scales and make it very difficult to predict predator distributions (Logerwell et al., 1998), i.e., the spatio-temporal relationship between predator and prey would weaken as there would be a large number of areas with an abundance of food but no penguins. Under this scenario, availability of prey to king penguins would become analagous to present prey availability to elephant seals now: the very deep diving capability of elephant seals gives them access to the entire mesopelagic zone and the DSLs there, and so it is very difficult to predict elephant seal dive locations from modelled prey distribution patterns because prey is, to an extent, ubiquitous (McMahon et al., 2019). The degree to which DSLs become shallower around South Georgia in the future will depend on whether their depth is driven by environmental conditions (e.g., light) or predation risk – if the latter, they might not change.

Implications for Marine Spatial Planning and Conservation Management

Protection of key Southern Ocean predators and prey is becoming increasingly important as fishing of lower trophic levels (e.g., squid and toothfish; Agnew et al., 2005; Kock et al., 2007) continues and other fisheries expand (e.g., Antarctic krill; Nicol et al., 2012). In this study, we have linked high intensity (and likely high biomass) DSLs with penguin habitat suitability. Penguins are likely drawn to these locations to feed on DSL-associated schools that occur within their dive range. This result supports the approach taken by Hindell et al. (2020) for defining AESs, which assumes that predator aggregations are indicative of elevated productivity and biomass at lower trophic levels. Since DSL intensity is readily predictable at local scales, conservation management could use modelled DSL intensity as a proxy for important predator foraging grounds to designate AESs. However, strong correlation between predators and lower-trophic level biomass is not ubiquitous (Boyd et al., 2015), and is likely to vary by species, season, access to habitat and prey distribution. Similarly, drivers of prey abundance are not always good predictors of predator distribution patterns, particularly in anthropogenically altered environments (Grémillet et al., 2008; Sherley et al., 2017). Therefore, it may not be prudent to define AESs over large areas based on ocean basin scale correlations. Instead, a regional approach, which considers seasonal variations in local populations (e.g., changes in predator distribution due to breeding stages, etc.), might provide a more informed basis for designating MPAs: the habitat suitability model developed in this study could be used to identify AESs within and beyond the SGSSI MPA (see Warwick-Evans et al., 2018; Handley et al., 2020).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://zenodo.org/rec ord/3587678#.YPArB-hKiF4 and https://doi.org/10.17630/d92c1af3-6149-4bb5-aa14-6fb00ef3a85a for the acoustic and penguin tracking data respectively.

Ethics Statement

The animal study was reviewed and approved by the Government of South Georgia and the South Sandwich Islands (regulated activity permit number 2016/048), the British Antarctic Survey and the University of St Andrews.

Author Contributions

RP collected the echosounder data, wrote the manuscript, and performed data analysis. CLG collected the echosounder data and king penguin data, performed preliminary analysis of king penguin and echosounder data, and produced a thesis chapter relating to the study. RBS collected the penguin data and contributed to writing the manuscript. AK helped with the accelerometery data analysis and contributed to writing the manuscript. YR-C provided guidance on data analysis methods and contributed to writing the manuscript. NR conceived the habitat preference analysis and contributed to writing the manuscript. SJ performed the DNA sequencing. AW performed the DNA extraction and contributed to writing the manuscript. JPYA provided the accelerometery tags and contributed to writing the manuscript. RAS provided literature information on diet study species and contributed to writing the manuscript. PGF collected echosounder data, performed the calibration at South Georgia, and contributed to writing the manuscript. LB assigned frontal zones to track positions and contributed to writing the manuscript. ASB conceived the study, planned the field activities and contributed to writing the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The at-sea data collection and 50% of CLG’s Ph.D. studentship was provided by the Swiss Polar Institute as a grant ‘Unlocking the Secrets of the False Bottom’ to ASB. The School of Biology, University of St Andrews, funded the other 50% of CLG’s studentship. Work at South Georgia was supported by the Natural Environment Research Council’s Collaborative Antarctic Science Scheme (CASS-129), a grant from the Trans-Antarctic Association grant to RBS, and a British Antarctic Survey Collaborative Gearing Scheme grant to RBS and ASB. ASB and RP were supported in part by UKRI/NERC under grant NE/R012679/1.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the staff at the British Antarctic Survey base at King Edward Point (South Georgia), Quark Expeditions and the crew and staff of the Ocean Endeavour and the FPV Pharos South Georgia for their help with the fieldwork logistics. We also thank the Swiss Polar Institute and the ACE foundation for funding our ACE project, and all our colleagues who assisted with acoustic data collection at sea: Matteo Bernasconi, Inigo Everson, and Joshua Lawrence. We thank Yves Cherel for fruitful discussion on the role of prey patches for king penguins in the Kerguelen region. We also thank C. Ribout and the Centre for Biological Studies of Chizé for conducting the sexing analyses of the birds.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.745200/full#supplementary-material

Footnotes

1. ^ http://marine.copernicus.eu

References