4.1. Study area
The study area encompasses the Bay of Biscay and adjacent waters in the North Atlantic Ocean, from 43 to 50˚ North and 0 to 10˚ West. The width of the continental shelf increases from South to North (from 30 to 180 km). The oceanic circulation is characterised by a weak anti-cyclonic circulation in the central zone and becomes cyclonic near the continental shelf. The water column of the Bay of Biscay is divided into four major water masses: (1) between 100 and 600 m, the water column has the characteristics of the central waters of the North Atlantic Ocean; (2) between 600 and 1500 m, Mediterranean waters flow from Gibraltar; (3) between 1500 and 3000 m, there are the deep waters of the Northeast Atlantic and (4) beyond 3000 m the deep Antarctic waters flow northward [35]. In this work, only waters from 0-2000 m were considered relevant as deep divers seldom forage deeper than 2000 m.
4.2. Data collection and collation
In this study, we used a part of the dataset assembled in Virgili et al. (31; Fig. 4), we only considered beaked and sperm whale sightings and effort data recorded in the Bay of Biscay. Visual shipboard and aerial surveys performed by seven independent organisations between 1998 and 2016 were assembled (details of the surveys are listed in Appendix A). A single common dataset was created, aggregating all survey datasets standardised for units and formats. Effort data were linearized and divided into 5 km segments using ArcGIS 10.3 [47] and the Marine Geospatial Ecology Tools software [48].
Cetacean sightings were recorded following line-transect methodologies that allow Effective Strip Width (ESW) to be estimated from the measurement of the perpendicular distances to the sightings [49], except for the JNCC-ESAS surveys that used a 300‐m strip‐transect methodology. We built a model to estimate the ESW for beaked whales and sperm whales, depending on observation conditions and survey types (following 31).
To account for the difficulty to identify them at species level, beaked whale species were pooled into one group including Cuvier’s beaked whales (Ziphius cavirostris), mesoplodonts(Mesoplodon spp.)and northern bottlenose whales (Hyperoodon ampullatus).
A total of 113 sightings of beaked whales and 52 sightings of sperm whales were assembled for the present study (Fig. 4b and c). Aggregated data represented about 150,400 km of on-effort transects (i.e. following a transect at specified speed and altitude or height above sea level, with a specified level of visual effort) of which 53% was carried out by boat and the rest by plane (Fig. 4a, Table 2).
Moran’s and Geary’s indexes were calculated to ensure there was no spatial autocorrelation in the data using the ‘spdep’ R-package [50].
Table 2. Effort performed by platform type and Beaufort sea-state per sector in the North Atlantic Ocean. This table presents the total effort conducted in each sector broken down by platform type and Beaufort sea-state. Beaufort sea-state values reported with decimals in the surveys were rounded up. For the analyses, all segments with Beaufort sea-state > 4 were excluded.
|
Total survey effort (km and %)
|
Total aerial effort (km)
|
Total shipboard effort (km)
|
Total effort by Beaufort sea-state class (km)
|
0-1
|
1-2
|
2-3
|
3-4
|
4-7
|
|
STUDY AREA
|
150,400
|
79,100
53 %
|
71,250
47 %
|
35,400
24 %
|
41,400
27 %
|
37,500
25 %
|
23,500
16 %
|
12,600
8 %
|
|
4.3. Data processing
4.3.1. Delimitation of depth classes
To determine whether the presence of deep-divers would be related to surface oceanographic processes or to processes taking place in the water column, four depth classes were delimited according to the water masses reported in the Bay of Biscay (see 2.1) in association with the dive profiles of deep-divers (Fig. 5). Environmental variables were extracted for each of the selected depth classes.
The first class is the “surface zone”; the deep-diver distribution may be related to mechanisms at the surface or these mechanisms may influence processes at depth. The second class is the epipelagic zone (named “0-200 m”), extending vertically from 0 to 200 m. We delineated this class between 0 and 200 m because most of beaked whale recovery dives occur between 0 and 200 m [51] and the upper 200 meters of the open ocean correspond to the euphotic zone where primary production is concentrated. Between 200 and 2000 m, we considered two water masses, the central waters of the North Atlantic Ocean between 200 and 600 m (names “200-600 m”) and the Mediterranean waters between 600 and 2000 m (named “600-2000 m”). From 1500 to 2000 m we considered no changes in the water masses and grouped the masses into one class.
4.3.2. Static and oceanographic variables
To model species distributions, it is necessary to extract environmental variables. We considered static and dynamic variables that can affect the distribution of beaked whales and sperm whales (Table 3). Static variables remain stable over time and are independent from depth classes while dynamic variables were extracted in each depth class and varied over time.
Table 3. Candidate environmental predictors used for the habitat-based density modelling. A: https://www.gebco.net/; 30 arc-second is approximately equal to 0.004°. B: Harris et al. [53). C: Copernicus (https://resources.marine.copernicus.eu/?option=com_csw&task=results&pk_vid=f7b9e59b6db4a7541611582309455951). All dynamic variables were extracted or computed for each depth class (surface, 0-200 m, 200-600 m and 600-2000 m). Abbreviations used in the following article are defined here in bold, d1-d2 refers to the depth classes e.g. 200-600 means between 200 and 600 m.
Variables used in the study with abbreviations and units
|
Original
Resolution
|
Sources
|
Effects on pelagic ecosystems of potential interest to deep-divers
|
Static
|
|
|
|
Depth (m)
|
15 arc sec
|
A
|
Deep-divers feed on squids and fish in the deep-water column
|
Slope (°)
|
15 arc sec
|
A
|
Associated with currents, high slopes induce prey aggregation or enhanced primary production
|
Roughness (m) -
|
15 arc sec
|
A
|
A high roughness indicates an important escarpment and a greater richness in prey
|
Surface of canyons – CanArea (km²)
|
15 arc sec
|
B
|
Deep-divers are often associated with canyons and seamounts structures
|
Dynamic
|
|
|
|
Mean and standard deviation of temperature – mTd1-d2 and sdT d1-d2 (°C)
|
0.083°, daily
|
C
|
Variability over time and horizontal gradients of temperatures reveal front locations, potentially associated with prey aggregation or enhanced primary production
|
Mean and standard deviation of gradients of temperatures – mGrT d1-d2 and sdGrT d1-d2 (°C)
|
0.083°, daily
|
C
|
Mean and standard deviation of eddy kinetic energy – mEKE d1-d2 and sdEKE d1-d2 (m².s-2)
|
0.083°, daily
|
C
|
High EKE relates to the development of eddies and sediment resuspension inducing prey aggregation
|
For static variables, we extracted the depth at a 15 arc second resolution (≈ 500 m; https://www.gebco.net/) and then computed the slope (inclination of the bottom) and roughness (difference between the maximum and minimum depth of the pixels surrounding the central pixel) with the function ‘terrain’ from the R package ‘raster’ [52]. We also extracted the surface area of canyons listed in the study area because deep‐divers are often associated with canyon structures [53; www.bluehabitats.org).
For dynamic variables, we extracted water temperatures and current vectors (U and V) for each depth class (surface, 0-200 m, 200-600 m and 600-2000 m) from the Copernicus website (https://resources.marine.copernicus.eu/?option=com_csw&task=results&pk_vid=f7b9e59b6db4a7541611582309455951). From these variables, we computed spatial gradients of temperatures and EKE (0.5*(U²+V²)). Gradients of temperatures were calculated as the difference between the minimum and maximum temperature values found in the eight pixels surrounding any given pixel of the grid (function ‘detectFronts’ from the R package ‘grec’, 54). Climatological variables were computed by calculating the means and standard deviations over the study period for each variable of each depth class (Appendix B).
4.4. Habitat-based density modelling
To model the distribution of beaked whales and sperm whales, we fitted GAMs [30] with a Tweedie distribution to account for over-dispersion in the cetacean count data [55] with the ‘mgcv’ R-package [56]. GAMs extend generalised linear models to allow for smooth, nonlinear functions of predictor variables to be determined by observed data rather than by strict parametric relationships [30]. The mean number of individuals per segment was linked to the additive predictors with a log-function with four degrees of freedom. An offset equal to segment length multiplied by twice the ESW, or twice the 300 m-strip for JNCC-ESAS surveys, was included (57; refer to Virgili et al. [31] for the ESW estimation; ESWs were calculated for each combination of platform, class of Beaufort sea-state and class of observation height). We removed combinations of variables with Pearson correlation coefficients higher than |0.5|and tested all models with combinations of one to four variables to avoid excessive complexity [58].
To assess the correlation between the variables we created a correlation matrix using the R package ‘corrplot’ (59; Fig. 6). Many variables with depths greater than 200 m were correlated with each other but also with surface variables (depth < 200 m) and with static variables such as slope, roughness and surface of canyons; showing the importance of considering correlations in the model selection.
The Akaike information criterion (AIC, the lower the better; 60) and Akaike weight (w; ‘akaike.weights’ function from ‘qpcR’ R package; 61) were used for model selection. The Akaike weight of a model is the probability that this model is the best among all tested models [60,62]. If the Akaike weight of a model is high (w > 0.9) it can certainly be identified as the best model and inference can be made from this model alone [43,60]. In contrast, if w < 0.9, a model averaging is recommended, which consists in producing parameter estimates from the weighted average of several models and not from a single model [43]. If w < 0.9, some variables are probably very similar and correlated. In this case, it is not possible to choose a best model among all tested models since the models are equivalent and they must all be integrated to produce an average prediction of species distribution. This can be cumbersome in terms of calculations and difficult to interpret, so if possible, it would be easier to obtain a single model.
To be able to identify this single suitable model, we aimed to approximate the average prediction obtained from the models by the prediction obtained from the model that combined the four most important variables. Following Symonds & Moussalli [43], we determined the importance of each variable by summing the Akaike weights of the models in which the variable was selected and ranked all variables. A model using the first four variables was then fitted by ensuring a non-correlation of the variables (if variables were correlated, the next uncorrelated ranked variable was chosen). A prediction was produced from this model and compared to the average prediction obtained from the models that explained 80% of the total Akaike weight. We considered only the models that explained 80% of the Akaike weight because beyond this threshold, the models were very negligible (very low explained deviances and very high AICs). If the coefficient of determination (R²) of the regression line established between the values of the average prediction and the values of the prediction obtained from the four most important variables is close to 1, predictions are similar and the average prediction of all models can be approximated by the prediction of the model fitted to the four most important variables. The four variables can therefore be used to obtain functional relationships and to predict the relative densities of beaked whales and sperm whales in the Bay of Biscay. If not, all models have to be considered to predict the species densities.
There were not enough data to fit a model by month or by season (the number of sightings in winter was too low) so we fitted models to all data of beaked whales and sperm whales and obtained climatological predictions maps for all seasons combined in the 1998-2016 period, although most of sighting data were collected in summer and the prediction maps most likely reflect the summer species distribution. The uncertainties associated with the predictions were also estimated as the standard deviation associated with the predicted relative densities; high values indicate high errors associated with density estimates. However, it should be noted that the uncertainty associated with the model prediction with the four important variables was certainly lower than the uncertainty associated with the mean prediction of all models and it is therefore underestimated.
Finally, the model fit of each model was assessed thanks to the percentage of explained deviance [30,63], the Root Mean Squared Error (RMSE) which measures the prediction errors and the model accuracy (the lower, the better, 64; ‘qpcR’ R-package, 61) and a visual inspection of predicted and observed distributions [65].