Highlights

  • Mapping of distribution of marine predators in one of the largest marine reserves in the world, using static, oceanographic and acoustically derived prey covariates.

  • Preyscapes derived from hydro-acoustic data improve predator niche modelling.

  • Spatial distributions of nine predator species are mapped.

Introduction

Modelling the habitat (for example, thermal preference) of top predators is critical to understand their spatiotemporal dynamics and to develop appropriate conservation strategies. However, the way top predators are monitored largely varies among groups and can make the modelling of their habitat a challenging task. The abundance of commercially exploited species is often assessed using catch data (for example, Lan and others 2018), and that of non-target species, using by-catch data (for example, Escalle and others 2019) and on-board observer data (for example, Lopez and others 2003). Non-targeted species such as cetaceans or seabirds are usually counted from boats or aircraft. Movements of exploited and non-exploited marine predators are also monitored using a variety of tags (for example, Block and others 2011; Abecassis and others 2015; Leroy and others 2016; Ravache and others 2020a). With some tagging methods, behaviour (foraging, migrating, resting) can be identified at various spatial–temporal scales: from a few seconds to several years and from a few centimetres to thousands of kilometres.

Moreover, although vertical and horizontal movements of marine predators are generally assumed to be mainly motivated by foraging behaviour (Lehodey and Maury 2010), proxies of prey availability and abundance are often absent from predator niche modelling (Robinson and others 2011). Such studies are generally restricted to oceanographic (for example, temperature, oxygen, chlorophyll-a concentration) and static (for example, depth, bottom slope) variables. The absence of prey in predator niche modelling results from difficulties in quantifying their abundance and characterizing their dynamics (for example, Escobar-Flores and others 2018; Proud and others 2018), as well as from variability in predator diet (for example, Potier and others 2007; Olson and others 2014; Duffy and others 2015). At the local scale, predator–prey relationships have sometimes been studied using acoustically derived prey distribution (for example, Bertrand and others 2002a; Torres and others 2008; Lezama-Ochoa and others 2010). At the regional scale, assessments of the influence of prey distribution on predator distribution generally do not use in situ prey information but, rather, ecosystem model output from which preyscapes are derived (a spatial layer of prey biomass or abundance; for example, Briand and others 2011; Lambert and others 2014; Miller and others 2018; Green and others 2020; Pérez-Jorge and others 2020). In the absence of in situ prey distribution observations to validate ecosystem model outputs, the reliability of such predicted preyscapes is questionable (Receveur and others 2020a).

In addition, diel vertical migration (DVM) of prey in the water column and the vertical behaviour of hunting predators influence predator–prey interactions (Schaefer and Fuller 2007, 2010; Benoit-Bird and McManus 2012; Choy and others 2017; Houssard and others 2017). Therefore, the vertical preferences and dynamics of organisms have to be taken into account in habitat predictive modelling (Opdal and others 2008; Louzao and others 2019a).

Horizontal and vertical preyscape structures around New Caledonia have been assessed using in situ acoustic data from at-sea surveys (Receveur and others 2019, 2020a). This region is characterized by a high richness of predators (for example, Ceccarelli and others 2013; Laran and others 2016; Borsa and Vidal 2018), which is supported by a diversity of prey (Allain and others 2012; Williams and others 2015; Receveur and others 2020b). This high diversity motivated the creation of the Natural Park of the Coral Sea covering the entire New Caledonia exclusive economic zone (EEZ, 1.3 million km2; Martin and Lecren 2014). The creation of this large Natural Park has emphasized the need for robust scientific information on the pelagic ecosystem status and potential anthropogenic threats (Ceccarelli and others 2013; Gardes and others 2014). In the present study, a statistical model combining environmental covariates and vertically structured prey abundance with top predator abundance and occurrence data was developed for New Caledonian waters. The benefit of including a proxy for prey abundance (the acoustic echo intensity integrated over different depth layers) in predator niche modelling, in addition to that of other marine physical variables, was quantified. Statistical relationships were then used to predict and map the spatial distributions of the nine selected marine predators.

Methods

Statistical models were fitted to indices of predator abundance and occurrence as a function of covariates that included prey distribution and other biotic and abiotic oceanographic conditions. Our domain of investigation covered the New Caledonia EEZ (14°S–27°S, 156°E–175°E, Figure 1).

Figure 1
figure 1

Spatial distribution of environmental covariates across the New Caledonia EEZ averaged across the 2010–2019 period: A bathymetry (in km), B sea surface temperature (SST, in °C), C depth of the 20 °C isotherm (in m) and D chlorophyll-a concentration (in mg m−3). Data sources are given in Table S1. Land is represented in grey, and white areas are waters with a bathymetry shallower than 300 m.

Full size image

Covariates

Five biophysical oceanographic variables and a proxy for biotic prey comprised the six environmental variables used to explore drivers of predator abundance/occurrence (Table S1). All covariates were extracted at the date and position of each predator observation. Bathymetry was extracted from the ZoNéCo database at a 500-m spatial resolution (ZoNéCo, 2013). Sea surface temperature (SST) was extracted from the NOAA OI SST high-resolution dataset (Reynolds and others 2007). The depth of the 20 °C isotherm was extracted from the Armor3D dataset (Guinehut and others 2012). Surface chlorophyll-a concentration was extracted from GLOBCOLOUR (Saulquin and others 2011). SST, depth of the 20 °C isotherm and chlorophyll-a concentration were extracted on a ¼-degree spatial grid, with a weekly temporal resolution (Figure 1).

Prey Data

The preyscapes were derived from the 3-D statistical model developed by Receveur and others (2019), which relates micronekton vertical distribution obtained from acoustic data to environmental conditions in the Coral Sea region. Acoustic data collection, statistical methodology and validation of micronekton biomass estimates in the New Caledonia EEZ are detailed in Receveur and others (2019) and in Supplementary 2. Acoustic data, covering more than 17,500 km of line transects, were collected during six cruises from 2011 to 2017 (Figure S3 and Table S3; Allain and Menkes 2011), using an EK60 echo sounder (SIMRAD Kongsberg Maritime AS, Horten, Norway). After processing, 38 kHz acoustic data were echo-integrated to provide the nautical area scattering coefficient (NASC, in m2 nmi−2), a proxy for micronekton biomass (Maclennan and others 2002; Irigoien and others 2014; Proud and others 2017). It is important, however, to keep in mind that gas-filled organisms (for example, fishes with a swim bladder) are largely dominant in the NASC38 kHz signal, whereas other organisms without gas are poorly represented (Davison and others 2015).

The final dataset was then composed of 16,715 acoustic vertical profiles deployed between 10 and 600 m, excluding transition periods at dawn and dusk. A machine learning model (XGBoost, a machine learning regression tree) was fitted to link the acoustic vertical profiles to environmental explanatory variables: norm of the wind, 0–600 m mean temperature, 0–600 m mean oxygen, 0–600 m mean salinity, surface chlorophyll, bathymetry and sun inclination (Figure S4, left column). This model enabled us to robustly predict (with an success rate of 87% of in cross-validation and Figure S5) a proxy for the vertical micronekton biomass (NASC profiles) for the entire New Caledonia EEZ.

NASC profiles were predicted on a regular ¼° spatial grid and on a weekly temporal resolution (Figure S4, right column) for day and night periods. Predicted NASC profiles were then integrated to obtain the day and night preyscapes over four vertical layers that were defined according to predators’ behaviour (see “Statistical analysis for niche modelling” section): 10–30 m (sub-surface; ‘sub_surf’), 10–200 m (that is, epipelagic layer; ‘epi’), 200–400 m (that is, upper mesopelagic layer; ‘up_meso’) and 400–600 m (that is, lower mesopelagic layer; ‘low_meso’) (Figure 2).

Figure 2
figure 2

Spatial distribution of micronekton biomass index NASC (nautical area scattering coefficient) values (in m2/nmi2) during the day (top row) and during the night (bottom row) integrated over four vertical layers and on average for the 2010–2019 period: 10–30 m (sub-surface, first column), 10–200 m (epipelagic, second column), 200–400 m (upper mesopelagic, third column) and 400–600 m (lower mesopelagic, fourth column) across the New Caledonian EEZ. Colour scale is different among panels. Land is represented in grey, and white areas are waters with bathymetry shallower than 300 m.

Full size image

Predator Data Collection

Tuna and Dolphinfish

Predator fish abundance was inferred from fishery daily catch rates (Table 1, Figure 3). Daily catch rates (number of fish per hundred hooks) of albacore tuna (Thunnus alalunga, ALB), yellowfin tuna (Thunnus albacares, YFT) and dolphinfish (Coryphaena hippurus, DOL) were obtained from the logbooks of all longliners operating in the New Caledonian EEZ from 2010 to 2018 (Anonymous 2018). The logbooks provided information on a daily basis: name of vessel, date and location (latitude and longitude) of the beginning of the fishing set, number of hooks and number of specimens captured per species. The logbook data from the New Caledonia fishing fleet are considered highly reliable; 100% of logbooks are received by the Direction des Affaires Maritimes, fishing locations recorded in logbooks are verified against Vessel Monitoring System data for each set, while catch species and quantities are verified by on-board fisheries observers and port samplers.

Table 1 Predator Variable Summary Indicating Taxon, Species Code, Logo used in Figures (Credit: see caption of Figure 3), Numbers of Observations (N) with Number of Occurrences (that is, When Individuals Were Detected for Cetaceans and When the Behaviour at the Position was ‘Foraging’ for Seabirds), Observation Unit, Data Source, Time of the Day, Sampling Years (Y) and Sampling Months (M, with 1 to 12 indicating January to December)
Full size table
Figure 3
figure 3

Raw data for the nine predators included in the study. Top row: catch per unit of effort (CPUE: number of fish caught per 100 hooks) for albacore (ALB), yellowfin tuna (YFT) and dolphinfish (DOL). Middle row: counts (grey points show absence of sighting) of Delphininae (DELPH), Globicephalinae (GLOB) and Ziphiidae (ZIPH, see section 2.1 for details on species). Bottom row: foraging (0: do not feed and 1: feed) behaviour positions obtained from GPS loggers in wedge-tailed shearwater (WTSH), Tahiti petrel (TAPET) and red-footed booby (RFBO). Seabird breeding sites are given in Figure S2. Credit: ©Les Hata, SPC (fish), ©Youngmi Choi, SPC (cetacean), ©Tubenoses Project, Hadoram Shirihai (WTSH and TAPET) and ©Andreas Ravache, IRD (RFBO). Land is represented in grey.

Full size image

Cetaceans

Cetacean counts were obtained from a large-scale aerial survey of the New Caledonia EEZ made between October and December 2014 (REMMOA survey: Laran and others 2016; Table 1). More than 30,000 km line transects were surveyed over a total surface area of 542,300 km2 (39% of the EEZ, Figure 3) at a ground speed of 90 knots and an altitude of 600 feet. A line transect methodology was used to count cetaceans (Buckland and others 2001) from high-wing double-engine aircraft equipped with bubble windows enabling vertical observation of the sea surface. Cetacean group size was recorded at the lowest possible taxonomic level, alongside observation conditions (Beaufort Sea state, glare severity, turbidity, subjective condition and cloud coverage). Cetacean observations were summed over 10-km segments of effort with homogenous weather conditions. Only segments with good weather conditions (that is, Beaufort Sea state ≤ 3 and subjective condition medium, good or excellent) were kept. Perpendicular distances from the track line were obtained from clinometer measurements by observers.

As the number of observations per species was relatively small, we grouped species into energetic guilds based on morphological features and diving capabilities (Lambert and others 2014; Mannocci and others 2014a, 2014b). These authors considered diving ability as a proxy for energetic requirements to classify tropical odontocete species into three guilds: Delphininae (DELPH), Globicephalinae (GLOB) and Ziphiidae (ZIPH). The DELPH guild was composed of small Delphininae [Stenella longirostris (9% of sightings), S. attenuata (9%) and unidentified Stenella spp.] and large Delphininae including Tursiops truncatus (21%), T. aduncus (19%), Lagenodelphis hosei (7%) and unidentified large Delphininae. The GLOB guild was composed of Globicephala macrorhynchus (52%), Grampus griseus (16%), Pseudorca crassidens (8%), Peponocephala electra (3%) and unidentified small and large Globicephalinae (3% and 16%, respectively). The ZIPH guild included Ziphius cavirostris (41%), Indopacetus pacificus (18%), Mesoplodon spp. (20%) and possibly other, unidentified beaked whale species.

Seabirds

Seabird positions were obtained from GPS-tracking surveys in 2012 and 2015 (in May, June and November) for the red-footed booby (Sula sula, RFBO; 21 individuals equipped) and in 2017, 2018 and 2019 (in January, June and from August to December) for the Tahiti petrel (Pseudobulweria rostrata, TAPET; 21 individuals equipped), both present year-round in New Caledonian waters (Table 1; Borsa and others 2010). The wedge-tailed shearwater (Ardenna pacifica, WTSH; 68 individuals equipped), present in the area during the breeding season only, that is, between October and May (Weimerskirch and others 2020), was tracked in 2017, 2018 and 2019 (in March and April).

GPS receivers were attached to adult WTSH during the chick-rearing period at five study colonies (Figure S2). The Pindaï study colony is located on New Caledonia’s main island [colony size (N) = 11,000 breeding pairs; see Table S2 for detailed references]. Three other study colonies (Tiam’bouene, N = 11,520; Canards, N = 340; and Mato, N = 2000) were lagoon islets, located less than 20 km from the main island. The farthest study colony was Loop islet in the Chesterfield atoll (N = 11,000), located 600 km NW of New Caledonia’s main island. Details concerning the deployment of the GPS receivers have been presented previously (Ravache and others 2020a; Weimerskirch and others 2020).

GPS-tracked TAPET were from the Mato (N > 20 breeding pairs) and Nemou (diffuse colony of a few hundred breeding pairs) study colonies, the latter being located in the eastern lagoon of New Caledonia (Figure S2; Ravache and others 2020b).

GPS-tracking experiments on adult RFBO were conducted on colonies from the Chesterfield atoll (N = 3540) in 2012 and Surprise islet (N = 2883) (Entrecasteaux atoll, 200 km north to New Caledonia’s main island) in 2015 (Mendez and others 2017, 2020).

The geographic distribution of the WTSH and TAPET study colonies was deemed representative of the populations of these two species across the New Caledonian EEZ (see Supplementary Table S2). The two RFBO colonies sampled were among the largest within the New Caledonian EEZ and were located in its northern half; two other major colonies located in the south-eastern part of the New Caledonian EEZ (Walpole and Hunter; Table S2) were not sampled.

The expectation–maximization binary clustering (EMbC) algorithm (Garriga and others 2016) was used to classify seabird movements in four behavioural categories. EMbC is a robust multivariate clustering algorithm based on trajectory sinuosity and speed, allowing identification of specific behaviours. It assumes that wide turning angles between consecutive locations indicate foraging behaviour, while fast and straight bouts are associated with commuting movement. One of the four following behaviours was assigned to each GPS position using the EMbC method: resting (low speed, small turning angle), commuting (high speed, small turning angle), extensive search (high speed, wide turning angle) or intensive foraging (low speed, wide turning angle). Movement bouts including at least three consecutive ‘intensive foraging’ locations (and ‘intensive foraging’ or ‘extensive search’ for TAPET, Ravache and others 2020b) were then conserved as ‘foraging presence’ locations, while the other locations were considered as ‘foraging absence’ locations for further analysis.

Statistical Analysis for Niche Modelling

Generalized additive models (GAMs; Hastie and Tibshirani 1995) and generalized additive mixed models (GAMMs) were used to examine relationships between predator abundance/occurrence (the response variable) and a suite of covariates (that is, bathymetry, SST, chlorophyll-a concentration, depth of the 20 °C isotherm d20 and NASC, Table S1). Absence of collinearity between pairs of covariates was tested using Spearman’s rank correlation test (below 0.5, Louzao and others 2011). All variables were smoothed with cubic splines with a maximum knot number set to eight for tunas and four for the seabirds and cetaceans due to the limited number of data, to prevent overfitting. Residual checks were conducted for the nine models using quantiles (Dunn and Smyth 1996): residual variance was constant, residuals were independent, and neither predicted values and residuals nor explanatory variables and residuals were statistically linked.

Models were fitted successively, excluding each of the remaining terms one at a time, to assess the percentage of variation explained by a predictor and to classify predictors according to the magnitude of their influence in the model.

Models for Tunas and Dolphinfish

For fishes, the response variable of the three GAMs was the number of specimens caught. The logarithm of the number of hooks was added as an offset to consider the variability in fishing effort among trips. A random effect by vessel was also fitted to deal with the variability in catchability among vessels. Year was added as a continuous explanatory variable to assess possible temporal change, and spatial smoothing in latitude and longitude was added to account for spatial autocorrelation (Table 2).

Table 2 Details of Habitat Distribution Models Tested: Predator Taxon, Response Variable, Link Function, Distribution, Offset to Correct for Sampling Effort, Explanatory Variables Other than the Nautical Area Scattering Coefficient (NASC), NASC Variable and Details of the Nested Model
Full size table

As the exact time of the catch by commercial longliners was unknown, NASC values integrated by depth class were averaged over day and night. Only the epipelagic layer was considered for DOL (Furukawa and others 2014; Lin and others 2019); the epipelagic and the upper mesopelagic layers were chosen for YFT (Dagorn and others 2006; Schaefer and others 2009), while the epipelagic and the upper and lower mesopelagic layers were used for ALB (Williams and others 2015).

Models for Cetaceans

For cetaceans, the response variable of the three GAMs was the number of individuals counted along the flight track per approximately 10-km-long segment. To estimate the sample area, effective strip width (ESW) was used, estimated as the point at which the number of schools detected beyond is equal to the number of schools that are missed within (Marques and Buckland 2004). The area sampled for each single segment was the length multiplied by twice the corresponding associated ESW. The log-transformed sampled area was added as an offset in the model to account for the variability in sampling effort. ESW for these datasets was estimated by Laran and others (2016).

As aerial surveys were conducted during the day, NASC values were used for the day period only and included the epipelagic, upper and lower mesopelagic layers (Lambert and others 2014).

Models for Seabirds

For seabirds, GAMMs were used to account for autocorrelation between consecutive behaviours in an individual. An autocorrelation structure of order two with a random effect fitted by individual was nested to deal with the correlation structure (Wood 2006; Dormann and others 2007). The response variable of the three GAMMs was the occurrence of foraging behaviour.

Using the GPS time stamp, we calculated sun inclination (Michalsky 1988; Blanc and Wald 2012) to estimate the period of the day (negative sun inclination values for night and positive values for day) of each spatial location. NASC values were then used for the corresponding moment of the day (day or night), and the sub-surface layer (10–30 m) was used, that is, the vertical layer used by these species to forage (Table 2). The first 10 m was missing as it is not recorded by the echosounder and therefore not present in the acoustic modelling (Receveur and others 2019).

Prediction Protocol

To avoid predictive extrapolation, predictions were limited to sampled environmental variable ranges and hence the presence of blank areas in the prediction maps, particularly visible for seabirds. Moreover, to avoid temporal extrapolation, predictions were made only for months when data were collected (given in Table 1).

The predictions were made on a regular 0.25-degree spatial grid based on environmental and NASC values, on a weekly temporal resolution, for all weeks of surveyed months between 2010 and 2019. For visualization, the predictions were then averaged over three-month periods (DJF: December, January and February; MAM: March, April and May; JJA: June, July and August; SON: September, October and November) over the 2010–2019 period. For the NASC values, the same vertical layers and night/day periods (that is, night-day mean for fish, day for cetaceans and night and day separated for seabirds) were used as for the fitted models (Table 2). As the fine-scale robust bathymetry used was available only for the New Caledonian EEZ (ZoNéCo 2013), predictions were limited to this EEZ.

The three fish models were used to predict catch values with an offset of number of hooks set to the mean value over the studied years (1960 hooks/longline set). To obtain catch per unit of effort (CPUE), a proxy for fish abundance, predicted catches in each grid cell were divided by the offset.

For cetaceans, the number of individuals was predicted with an offset of sampled area set to a mean value (4 km2). To obtain density of individuals (no./km2), predicted numbers of individuals in each grid cell were then divided by the offset.

Finally, for the three seabird species, the model directly predicted the probability of occurrence of foraging behaviour. The breeding parents are spatially limited to the region they can reach by their flying ability and by the time that they can spend away from their chick. Therefore, the maximal distance travelled from the colony was calculated based on the datasets for each species (774 km for WTSH, 1079 km for TAPET and 320 km for RFBO), and the predictions were spatially limited to the circle with a radius equal to these maximal distances. For WTSH and TAPET, the predictions were limited by this maximal distance around the sampled colonies. For the RFBO, three well-known colonies were added in the south of the EEZ (Figure S2), and the predictions were limited around these five (two sampled and three well-known) colonies.

All analyses were conducted using R version 3.5.0 (R Core Team 2016), and GAMs were done using the mgcv package (Wood 2017). The R script used for modelling and predictions is given in Supplementary 3.

Results

Environmental and Prey Layers

Bathymetry is complex in the New Caledonia EEZ (Figure 1A), with deep waters in the north and the south-east (~ 5 km), shallower waters in the south-western corner and south of the main island (~ 3 km) and the deep New Hebrides trench in the east (> 7 km deep; Stewart and Jamieson 2018).

On average, for the 2010–2019 period, relatively warmer SSTs occurred in the north (~ 24 °C) and colder SSTs in the south (~ 18 °C) (Figure 1B). Mean depth of the 20 °C isotherm was about 160 m in the south and about 210 m in the north with the deepest values (~ 230 m) located west of the Chesterfield reefs (Figure 1C). Chlorophyll-a concentration was high close to coasts and reefs and generally higher in the south of the main island compared to the north.

As for SST and the 20 °C-isotherm depth, NASC spatial distribution showed a strong latitudinal gradient, especially during the night, with values in the south almost twice as high as values in the north on average for the 2010–2019 period (Figure 2). During the day, between 0 and 30 m, the largest NASC values were close to the coast of the main island and in the south. The 200–400-m integrated NASC showed similar spatial patterns. The 0–200-m and 400–600-m integrated NASC was higher in the south than in north.

Predator Data

ALB, YFT and DOL datasets were composed of 14,941, 17,420 and 17,106 longline sets, respectively, collected between January 2010 and December 2018. The number of fish caught per 100 hooks, calculated by set, ranged between 1 and 26.5 (ALB), 1 and 9.1 (YFT) and 1 and 3.3 (DOL). The EEZ was fished across its entirety, except for small regions in the south-eastern corner and in the extreme north (Figure 3).

Cetacean data were distributed more sparsely, with spatial and temporal observation coverage lower than for fish. There were around 2500 segments, but with a low number of sighting occurrences (between 31 and 40 according to the cetacean guild). The encounter rates (number of individuals per sampled km2) ranged from 0 to 26.7 for DELPH, from 0 to 163.8 for GLOB (with only two rates larger than 10) and from 0 to 1.5 for ZIPH (Figure 3).

For seabirds, temporal sampling effort was limited by the breeding season, when animals were present in New Caledonia (Table 1). The number of recorded GPS positions varied substantially among species, with more than 11,500 observations for WTSH (including 3934 eating occurrences), more than 6500 for RFBO (including 1547 eating occurrences) and more than 4500 for TAPET (including 1586 eating occurrences).

Environmental Drivers of Predator Distribution

Importance of Prey Distribution

Total deviance explained by abundance/occurrence models ranged from 7.3% (DELPH) to 31.7% (GLOB) (Table 3). The NASC variable had a relatively small influence in all models, although it ranked second in WTSH and RFBO, third in YFT, fourth in DOL and fifth in ALB (Table 3). Relationships to NASC were weak and non-significant for DELPH, GLOB, ZIPH and TAPET (Figure 4) with very low deviance explained by this variable (Table 3).

Table 3 Total Deviance Explained (Second Column) by Predator Abundance/Occurrence Models and for the Five Explanatory Variables Ordered by the Importance of Their Influence
Full size table
Figure 4
figure 4

Modelled responses to the prey biomass index NASC (nautical area scattering coefficient, in m2/nmi2) variations: CPUE for tunas and dolphinfish, counts of individuals for cetaceans and foraging probability for seabirds. Colours indicate different vertical layers of prey for tuna and cetaceans and different times of day for seabirds. Solid grey ribbons correspond to the confidence limits (± 2SE) of the model. Stars coloured by the factor (vertical layer or moment) indicate significant relationships (p value < 0.01). The ticks on the x-axis indicate the position of the observed values. Species codes on the y-axis are detailed in Table 1. ©Les Hata, SPC (fish), ©Youngmi Choi, SPC (cetacean), ©Hadoram Shirihai/The Tubenoses Project (WTSH and TAPET) and ©Andreas Ravache, IRD (RFBO).

Full size image

ALB abundance was positively linked to NASC values in the epipelagic and lower mesopelagic layers: the higher the NASC values, the larger the predicted CPUE (Figure 4), and negatively linked in the upper mesopelagic. DOL abundance was positively linked to epipelagic NASC values. YFT relationships with epipelagic and upper mesopelagic NASC were negative, with a steeper slope for the epipelagic than for the upper mesopelagic NASC. For ALB and DOL, NASC influence was relatively small compared to the total deviance explained, while it was five times higher for YFT (Table 3).

WTSH foraging probability was linearly positively related to 10–30 m NASC values during the day, while the relationship exhibited a bell shape during the night (Figure 4). The NASC variable was the second most important driver of this species’ foraging probability as it explained 6.7% of the deviance (Table 3). RFBO foraging probability was negatively linked to NASC values during the day (Figure 4), and NASC relative importance in RFBO model ranked second with 5.7% explained deviance (Table 3). The flat relationship between NASC and TAPET foraging probability demonstrated the non-significance of the relationship (Figure 4).

Other Oceanographic Drivers

For the predatory fishes, the four other environmental covariates were almost all statistically significant but their influence ranked differently across species (Table 3). YFT and DOL CPUEs were most influenced by chlorophyll-a concentration and ALB CPUE by SST. Bathymetry was the second influencing factor for ALB and YFT, and SST ranked second in DOL.

Optimal environmental values varied among fish species (Figure 5): higher ALB catch was predicted for intermediate bathymetry (around 4 km), whereas YFT catch was highest for shallower bathymetry (shallower than 2 km), and DOL catch was highest for bathymetry deeper than 4 km depth. Maximum YFT catch was predicted for intermediate chlorophyll values (0.2–0.4 mg m−3), whereas DOL catch was highest for low chlorophyll-a concentration values (below 0.1 mg m−3). The relationship with chlorophyll for ALB was almost flat. Cold SST (below 22.5 °C) was favourable to ALB catch, whereas warm SST increased YFT and DOL catch (about above 27 °C). ALB and DOL catch increased with 20 °C-isotherm depths greater than 200 m depth; conversely, YFT catch was higher when the 20 °C isotherm was shallower than 200 m depth (Figure 5).

Figure 5
figure 5

Modelled responses of predators (by row) to four environmental variables (by column). The solid grey ribbons correspond to the confidence limits of the model (± 2SE). Stars indicate the significance level: ** indicates highly significant (p value < 0.01) and * indicates slightly significant (0.01 < p value < 0.1). The ticks on the x-axis indicate the observed values’ position. Species codes are detailed in Table 1.

Full size image

DELPH abundance was mostly influenced by chlorophyll-a concentration; GLOB abundance by the 20 °C-isotherm depth; and ZIPH abundance, by chlorophyll-a concentration (Table 3). Higher DELPH abundance was predicted in waters with low chlorophyll-a concentration (below 0.1 mg m−3) and cold SST (< 25 °C). GLOB abundance was predicted to be higher when the 20 °C isotherm was deep, when the depth was either shallower than 2 km or deeper than 4 km, and in waters with low chlorophyll-a concentration. ZIPH abundance was predicted to be higher in waters with intermediate chlorophyll-a concentrations (about 0.8 mg m−3), around 3 km deep, with warm SST and with a 20 °C-isotherm depth at around 150 m (Figure 5).

Foraging behaviour of the three seabird species was primarily influenced by bathymetry, but with contrasting optimal values: WTSH mostly fed where waters were deeper than 2 km, while TAPET and RFBO fed in waters shallower than 2 km (Figure 5). Secondarily, WTSH foraging behaviour was positively influenced by high NASC values (above 5 m2 nmi−2) integrated between 10 and 30 m deep, and when 20 °C-isotherm depth was shallower than 150 m. TAPET was predicted to more likely feed in waters with chlorophyll-a concentration > 0.2 mg m−3. RFBO foraging behaviour was also influenced by SST, with more feeding predicted to occur in warm waters in general.

Spatial Distribution Predictions

Predator Spatial Distribution

The preferred habitats for ALB, YFT and DOL, indicated by predictions of high CPUE, did not spatially overlap (Figure 6). Higher ALB CPUE was generally predicted west of New Caledonia’s main island. The largest CPUE was predicted west of the Chesterfield reef, except for a small area of low CPUE directly surrounding the Chesterfield and Bellona atolls. (Reef and island names are indicated in Figure 1A.) Conversely, high predicted YFT CPUE was predicted in this small area. Higher YFT CPUE was predicted close to the lagoon, around the main island, especially in the north, around the d’Entrecasteaux atolls, and to a lesser extent around the Loyalty Islands. Finally, the channel between the main island and the Loyalty Islands was predicted to be the most favourable for high DOL CPUE. Predicted standard errors were larger in the south-east of the New Caledonia EEZ and lower than 10% for the rest of the EEZ (Figure S6).

Figure 6
figure 6

Spatial predictions of fish abundance index by quarter in the New Caledonia EEZ: catch per unit of effort (CPUE; number of fish caught per 100 hooks) of albacore (ALB), yellowfin tuna (YFT) and dolphinfish (DOL), with different scales according to species, and on average for the 2010–2019 period. DJF: December, January and February; MAM: March, April and May; JJA: June, July and August; SON: September, October and November. For each quarter, the three months were used for making predictions, as all months were sampled. Land is represented in grey. Reef and island names are indicated in Figure 1A. White areas (excepting waters with bathymetry shallower than 300 m identified in Figure 1) have no predictions because no extrapolation was performed.

Full size image

Higher DELPH abundance in the EEZ was predicted to the south-east and south of the main island between September and November, while a hotspot was predicted for this cetacean guild in the northern part of the EEZ between December and February (Figure 7). GLOB were predicted to occur most frequently in the extreme north of the EEZ and around the Lansdowne bank. ZIPH were predicted to occur around Chesterfield reefs and between the main island and the Loyalty Islands between September and November and especially between the main island and the Loyalty Islands between December and February. For the three cetacean guilds, predicted standard errors were high, on average larger than 50% for DELPH and GLOB, certainly linked to the (low) performances of the models (Figure S7).

Figure 7
figure 7

Spatial predictions of cetacean abundance by quarter in the New Caledonian EEZ: individual counts of Delphininae (DELPH), Globicephalinae (GLOB) and Ziphiidae (ZIPH), with different scales according to species, and on average for the 2010–2019 period. DJF: December, January and February; MAM: March, April and May; JJA: June, July and August; SON: September, October and November. Italic months in the right top corners indicate the sampled month and therefore the month used to do predictions. Reef and island names are indicated in Figure 1A. White areas (excepting waters with bathymetry shallower than 300 m identified in Figure 1) have no predictions because no extrapolation was performed.

Full size image

WTSH were predicted to feed mostly north of 19°S when the species is present in New Caledonian waters. TAPET mostly fed in waters south of the main island in shallow waters (close to the lagoon), south of Bellona and above the Lansdowne bank, between Chesterfield and the main island (Figure 8). Predicted standard errors for WTSH and TAPET were below 20% (Figure S8). No clear foraging hotspot was identified for RFBO: predictions showed scattered areas around the EEZ mid-latitudes, always close to the two studied colonies; but predicted standard errors were higher than 60% during the period extending from June to November.

Figure 8
figure 8

Spatial predictions of seabirds foraging probability by quarter in the New Caledonian EEZ: foraging occurrence of wedge-tailed shearwaters (WTSH), Tahiti petrels (TAPET) and red-footed boobies (RFBO), with different scales according to species, and on average for the 2010–2019 period. DJF: December, January and February; MAM: March, April and May; JJA: June, July and August; SON: September, October and November. Italic months in the right top corners indicate the sampled month and therefore the month used to do predictions. Reef and island names are indicated in Figure 1A. White areas (excepting waters with bathymetry shallower than 300 m identified in Figure 1) have no predictions because no extrapolation was performed.

Full size image

Discussion

Previous studies on marine predators in the tropical south-western Pacific have focused on specific predator groups including seabirds (Pandolfi-Benoit and Bretagnolle 2002; Borsa and others 2010; McDuie and others 2015; Mendez and others 2017, 2020; Ravache and others 2020a, 2020b Weimerskirch and others 2020), cetaceans (Garrigue and others 2015; Laran and others 2016; Derville and others 2019, 2020), sharks (Tirard and others 2010; Boussarie and others 2018) and albacore tuna (Briand and others 2011; Williams and others 2015) in separate publications. This study represents the first attempt to reach an integrated, quantitative, multispecies view of marine predator trophic ecology across the Coral Sea pelagic ecosystem. Prey abundance as estimated through acoustic surveys was tested together with oceanographic covariates to improve predator distribution models. Preyscape significantly affected five predator distribution patterns out of nine, but their influence remained small relative to other predictors and they did not rank as the most important variable for any species.

Sources of Uncertainty Linked to Predator Data

Fishing vessels in New Caledonia are longliners with similar fishing capabilities (Anonymous 2018). By including a random ‘vessel’ effect, we considered the potential variability of fishing power between vessels. By including the effort (that is, number of hooks) as an offset, we transformed catches into CPUE, which was then used as a proxy for fish abundance. Although large differences in CPUE can be erroneously measured where only small differences in abundance exist (Kleiber and Maunder 2008), our modelling approach incorporating a large dataset of fishing operations is likely a reliable proxy for fish abundance.

Compared to CPUE data, the three-month aerial survey REMMOA, which only occurred once in New Caledonian waters, only gives a snapshot of animal presence and partial distribution over a short period. Cetaceans are highly mobile predators; they may have been travelling between favourable habitats when spotted and therefore the oceanographic conditions attributed to that group were not necessarily optimal foraging habitats. Furthermore, aerial surveys cannot detect the presence of animals when they are diving, which then represent false absence values in the model. Improving our niche modelling of cetacean habitats may require additional aerial survey effort to better characterize their distribution.

For seabirds, uncertainty was linked to the clustering algorithm that identified behaviour based on trajectory sinuosity and speed. Here, by selecting only the intensive foraging observations occurring for more than three consecutive positions, we ensured that the animals were robustly classified as foraging.

In the present study, we mostly selected adult individuals for fish (through catch) and seabirds (through GPS tracking). For future studies, it would be important to include juvenile, immature and non-breeding animals which represent a significant part of a population. Immatures and non-breeders may use different habitats from adults (Mendez and others 2016, 2020; Houssard and others 2017; Weimerskirch and others 2020). However, this objective would require the use of different sampling techniques and the development of modelling strategies mixing different sampling methods (for example, aerial count for seabirds together with tracking methods) for one species.

Acoustic Micronekton Layers for Predator Niche Modelling

Our approach relied on a large dataset to robustly predict preyscape at a large spatial–temporal scale. For a given observed acoustic value, the predicted acoustic value varied by ± 10 dB (Receveur and others 2019; Figure S5), which could result in large biomass changes (Benoit-Bird 2009; Proud and others 2019). However, predicted acoustic spatial and temporal variations matched observed variations (see Figure 11 of Receveur and others 2019).

The use of the 38 kHz frequency favoured the detection of organisms with gas-filled organs (for example, fishes with a swimbladder) while excluding organisms without gas (Davison and others 2015), which may bias our perception of the spatial variation of prey abundance (Proud and others 2018, 2019). Nevertheless, given the large diversity in predator diet (Olson and others 2014; Duffy and others 2015; Williams and others 2015) and the difficulty in transforming acoustic signal into species-specific biomass (Davison and others 2015; Proud and others 2019), we used NASC values as a proxy for the relative abundance of prey available to all predators. Separating the acoustic signal into species groups (such as ‘fish’, ‘gelatinous’ or ‘crustacean’) would allow linking each predator to its favourite prey. We assumed that a large part of the weakness of the trophic relationships was linked to this acoustic value, which does not fully represent the specific diet of each predator.

Another acoustic limitation was linked to the fine horizontal scale (for example, prey aggregations in shoal) that was removed from the acoustic signal in the statistical framework, although it is an important component of foraging behaviour (Benoit-Bird and Au 2003; Bertrand and others 2014; Grados and others 2016).

More positively, one strength of our acoustic modelling framework was its ability to average relative prey abundance over different vertical layers. Prey layers were defined according to the predators’ vertical feeding habitats as reported from the literature. The differences in spatial patterns during the day of integrated NASC values over depth layers 10–30 m and 10–200 m (Figure 3) confirmed the relevance of extracting specific vertical layers to analyse prey–predator relationships. Integrating the vertical dimension is a necessity to define the 3-D oceanographic habitats used by marine species (Lambert and others 2014; Brodie and others 2018; Louzao and others 2019a).

In previous studies, two approaches were used to implement prey abundance in predator niche modelling: (1) by simultaneously collecting prey data and observing megafauna (for example, Doray and others 2018; Louzao and others 2019b; Putra and Mustika 2020) and (2) by using SEAPODYM, an end-to-end ecosystem model (Lehodey and others 2008, 2010) prey output (Briand and others 2011; Lambert and others 2014; Miller and others 2018; Green and others 2020). Our study suggests a third alternative where prey abundance is estimated during previous acoustic surveys focusing on low trophic levels, with limitations detailed above, and then compared to predator observations from other surveys. Such an approach could be useful for predator sampling techniques with which it is impossible to collect prey data at the same time (for example, aerial survey and telemetry). Our approach enabled us to study predator–prey relationships at the regional scale based on in situ prey data, but required acoustic observations at this vast spatial scale. The use of a machine learning algorithm to extrapolate the acoustic signal to the whole EEZ helped make robust predictions, with an 87% success rate in cross-validation. Receveur and others (2021) validated for the Coral Sea region the extrapolation of acoustic signal in the epipelagic layer by comparing it to the outputs of an ecosystemic model (SEAPODYM). Conducting more acoustic surveys through, for example, routine acquisition by commercial ships, and improving the relationships between acoustic signal and species composition and biomass, would help refine prey distribution models which, in turn, would improve prey implementation into predator niche models.

Predators Habitat

Tunas and dolphinfish were strongly influenced by SST: YFT and DOL preferred warm SST, whereas the optimal SSTs for ALB were colder (Bertrand and others 2002b; Kleisner and others 2010; present results). ALB’s relationships to NASC values were in agreement with previous results found in the New Caledonia EEZ using the SEAPODYM outputs as prey data source (Briand and others 2011). Briand and others (2011) showed a positive relationship between ALB CPUE and epipelagic layer prey biomass and found a negative relationship to upper mesopelagic prey biomass, similar to our results (Figure 4).

To explain the negative relationship to upper mesopelagic prey, these authors assumed that the increase in prey abundance induced a competition for baited hooks located in the mesopelagic layer. This hypothesis was confirmed by a study on albacore tuna in French Polynesia (Bertrand and others 2002b). This hypothesis would also explain the negative relationships of YFT catch values to NASC values. ALB, YFT and DOL exploit comparable foraging resources and have similar niches (Teffer and others 2015), especially YFT and DOL in the Pacific (Young and others 2010). Here, we instead observed a relatively weak distributional overlap between the three species (except for DOL and YFT east of the New Caledonian main island, Figure 6A) and different relationships to environmental predictors. We hypothesize that the data obtained at the scale of the Pacific Ocean for ALB, YFT and DOL show roughly similar distributions but that at a smaller scale, the spatial distributions of the three species can be distinguished based on their habitats. For instance, Bertrand and others (2002a; b) showed that tunas were more abundant in waters rich in prey and with favourable hydrological conditions at a regional scale (for example, Polynesian EEZ scale) and were attracted by small prey patches at smaller scale in the south-east Pacific. We confirmed this result here for ALB, but found that YFT was more abundant in waters with relatively low prey abundance. The 20 °C-isotherm depth was relatively important to explain the presence of the two tuna species. It is a proxy for vertical habitat change, but it is strongly linked to catchability changes. Therefore, part of the abundance changes linked to the 20 °C-isotherm depth may also reflect catchability changes rather than a true ecological relationship.

Chlorophyll-a concentration and the 20 °C-isotherm depth (GLOB and ZIPH), together with SST (DELPH), influenced the distribution of the cetacean predators in New Caledonia, as previously shown elsewhere in a tropical region (Mannocci and others 2014a). The three bell-shaped relationships with the chlorophyll-a concentration identified here, were in agreement with previous published work in the circumtropical belt (Mannocci and others 2015). This study also emphasized the non-negligible influence of the depth of the minimum dissolved oxygen concentration, a variable that we did not include. The absence of this parameter may explain the relatively low model performances of the three cetacean guilds of the present study. The positive linear relationship between ZIPH abundance and SST was also already demonstrated in the South Pacific gyre (Mannocci and others 2014a). NASC values were non-significant for the three guilds in our study. This must be linked to the acoustic signal. Indeed, as discussed in “Acoustic-Micronekton Layers for Predator Niche Modelling” section, the acoustic signal encompasses a large variety of micronekton species dominated by gas-filled organisms. However, cetacean diet species are likely poorly represented in the acoustic signal, explaining the absence of a significant relationship. The SEAPODYM micronekton model includes these cetacean diet species, as Lambert and others (2014) found significant trophic relationships for cetacean; thus, that may be a more suitable approach to study cetacean trophic relationships. Lambert and others (2014) also noted that delphinids mostly forage at night. The use of day-only NASC values in our modelling may also limit the strength of the fitted relationships.

Finally, cetacean species were grouped based on morphological and taxonomic characteristics, to delimit energetic guilds, while other taxonomic groups were studied at the species level. This grouping leads to a mix of different habitats while adjusting a single habitat model, which might have contributed to the low explained deviances. However, these guilds also correspond to energetic considerations and, therefore, the grouped species have a similar habitat. In addition, this strategy is needed to increase the robustness of the model considering the low number of sightings by species. Finally, this approach allowed the full exploitation of aerial survey data in which identification to the species level was not always possible.

Seabird relationships to oceanographic conditions were in agreement with previous results in the western Coral Sea (McDuie and others 2018). Among the variables tested, bathymetry had the strongest influence on all three species, as already demonstrated for WTSH in the Indian Ocean (Catry and others 2009). Relationships to the 20 °C-isotherm depth were mostly negative for WTSH and RFBO. Such relationships may reflect the fact that at a large scale, shallower 20 °C isotherms characterize regions with usually higher primary production. Moreover, the 20 °C-isotherm depth is linked to the position of meso-scale eddies, which have been shown to influence foraging behaviour of seabirds (Tew Kai and others 2009; Tew Kai and Marsac 2010; Jaquemet and others 2014; McDuie and others 2018). Indeed, meso-scale eddies potentially structure vertical and horizontal prey abundance, with large variability of demonstrated effects (Sabarros and others 2009; Godo and others 2012; Della Penna and Gaube 2020). More detailed spatial data and localization of eddies would be required to understand these relationships. The difference in NASC relationships between day and night in WTSH suggested that WTSH can locate prey concentration near the surface during the day, whereas other feeding techniques might be used during the night. This confirmed the results of Ravache and others (2020a), who found that WTSH feeding behaviour was higher during the full-moon nights due to higher luminosity. TAPET did not target areas with high prey abundance but foraged over shallow waters with high chlorophyll-a concentrations. The absence of a relationship to NASC values confirms results showing that TAPET have a scavenging foraging behaviour and mostly target isolated prey (Ravache and others 2020b). For the RFBO, the radius of the prospective habitat reached is smaller than for the two other seabirds species (Mendez and others 2017, 2020; Ravache and others 2020b; Weimerskirch and others 2020; present results). This can explain the large confidence interval predicted around the NASC relationship: the RFBO individuals found their food at proximity and are not able to prospect on larger scale to target high prey concentration.

Conclusion

We developed a niche modelling approach to evaluate the role of environmental predictors, including preyscape and other oceanographic conditions, to characterize the spatial distribution of nine marine predators in the south-west Pacific region. The inclusion of the different preyscapes based on in situ data was especially innovative and interesting. Our approach illuminated different trophic interactions and behaviours (for example, different across the species considered, different across the vertical layers, different between day and night). To further develop this approach, a variety of additional taxa (for example, sea turtles, sharks) could be added and would provide a more exhaustive overview of predator–prey relationships in the pelagic ecosystem.