The known sensitivity of beaked whales to anthropogenic sounds emissions [1-5] has stimulated worldwide the interest about the ecology of this species. A key element in the success of managing species of conservation concern is the knowledge of the species’ habitat use and preferences. Knowledge on the habitat use and preferences of beaked whales has been acquired worldwide in the last twenty years [6-13]. Several studies have reported a strong relationship with the topographic features of the sea bottom as animals were regularly observed over the continental slope in waters up to 2000 m of depth [11,14-16] and near submarine canyons [17]. Very few of these studies [10,17] have considered habitat predictors other than topographic or remote sensed parameters.

Cuvier’s beaked whale (Ziphius cavirostris G. Cuvier, 1823), is the only beaked whale regularly inhabiting the Mediterranean Sea area, where this species has been found associated with continental slope and with submarine canyons and seamounts areas [18-28]. In particular, the Genoa canyon area, located in the Ligurian sea (north western Mediterranean Sea), has been identified as a key habitat for Cuvier’s beaked whales [13].

Submarine canyons represent heterogeneous environments not only for their complex topography, but also for the fact that they are characterized by peculiar patterns of hydrography, flow sediment transport and accumulation [29]. Physical conditions in canyons are in fact often distinct from the surrounding shelf and slope, and can affect the structure of canyon communities [30]. Physical oceanographic conditions inside canyons, such as accelerated currents, upwelling and dense-water cascades, are often caused by topographic and climate forcing. These phenomena can be responsible for increasing suspended particulate matter concentrations and transport of organic matter from coastal zones to the deep ocean, enhancing both pelagic and benthic productivity inside canyon habitats [31]. Consequently, ocean dynamics may affect the distribution and behaviour of Cuvier’s beaked whales and the organisms forming the food web upon which the whales feed, through mechanisms such as concentration of nutrients and increase of primary production or aggregation of biomass due to convergence of water masses [9,31-33]. The opportunity to investigate how ocean dynamics influence the habitat selection of deep diving species is often limited by the availability of systematic oceanographic data collected in concurrence with marine mammal surveys. However, by collecting detailed oceanographic data, it is possible to identify additional key habitat characteristics. The aim of this study is therefore to correlate Cuvier’s beaked whale presence to the in situoceanographic measurements collected to a depth of about 1200 m during two dedicated research cruises conducted in 2002 and 2011 in the Genoa canyon area. This is the first time that measured oceanographic data are employed to model beaked whale habitat in the Mediterranean Sea. The specific objectives of this study are: 1) To evaluate whether dynamic parameters (i.e., temperature, salinity, sound velocity, dissolved oxygen, fluorescence and turbidity) may be used as potential predictors of the Cuvier’s beaked whale habitat; 2) To validate the prediction accuracy of the oceanographic-based models on an independent dataset collected in the same area; and 3) To improve the understanding of beaked whale key habitats characteristics in the Genoa canyon area.

Study area

The Gulf of Genoa (Figure 1) is located in north-western portion of the Ligurian Sea and is contained within the “International Sanctuary for the Protection of Mediterranean Marine Mammals’’ also known as “Pelagos Sanctuary”. Several canyons characterize the Gulf of Genoa with very steep slope gradients extending from the shelf break to a depth of about 2000 m. The Genoa canyon is the largest and the most northern canyon of the western Mediterranean Sea. Genoa canyon has its axis oriented northeast–southwest, with its two main canyons, the Polcevera and Bisagno, found in the head. These two canyons exhibit a linear along-axis topographic V shaped profile, are more than 700 m deep, 20 km wide and about are 60 km in length. Their steep walls suggest they are strongly affected by land sliding processes [34]. Directly east of this region is a wider canyon with a wide shelf to its south. The western part of the valley has a steep slope and several small canyons cut it. To the southwest, the canyon system descends into a deep abyssal plain. This large submarine valley (hereafter called “Genoa canyon area”) forms a boundary for the predominant circulation [35]. The circulation in the Ligurian basin consists of a basin-wide cyclonic gyre [36,37], which extends over the upper 500 m and can spread out to the west in to the Catalan Sea.

Figure 1: a) The Gulf of Genoa study area. b) Western Mediterannean Sea: Pelagos Sanctuary 638 borders (thick black lines) are shown.

Two major northward-flowing currents influence the general circulation in the basin: The Eastern Corsican Current that brings Tyrrhenian water in to the Ligurian Sea through the Corsica Channel and the Western Corsican Current which is part of the large cyclonic circulation. When these water masses join together north of Corsica, they form the Northern Current, also called the Liguro-Provencal- Catalan Current, which flows along the slope of the Italian, French and Spanish coasts [38]. A frontal region is commonly found at the limit of the cold core of the cyclonic Ligurian Sea circulation that separates the warmer, less saline peripheral zone near the coast from the denser and colder water in the central zone of the Ligurian Sea [39-41]. The occurrence of this permanent frontal area leads to an upwelling of deeper waters in the central zone and downwelling along the peripheral zone [39]. The upwelling provides the surface layer with inorganic nutrients that support a spring primary production which is higher than the average for the western Mediterranean waters and allows the fresh organic material produced near the surface to sink to the depths, enriching the entire water column and enhancing biological processes [40,35]. The Gulf of Genoa is north of the frontal area, however, the distribution of autotrophic biomass suggests an enriched zone in the southern part of the study area [35].

Data collection

From 1999 to 2011, more than 10 large-scale interdisciplinary sea trials were conducted in the Mediterranean Sea under the framework of the NATO Undersea Research Center (NURC)¹ Marine Mammal Risk Mitigation Project (hereinafter called MMRM). This long-term study was designed to better understand the distribution of the cetacean population in the Mediterranean Sea, to develop analytical models to predict species presence and to test and improve passive acoustic monitoring methods. The majority of the at sea campaigns, entitled SIRENA, were conducted in the Ligurian Sea area. In this study data collected from the NATO Research Vessel (NRV) Alliance in the Genoa canyon area in summer 2002 (SIRENA02/phase B) and 2011 (SIRENA 11/phase A) are analysed. Among the different sea trials conducted in the Ligurian Sea, these two campaigns were selected for the analysis because the sampling station protocol, oceanographic stations locations and time period were considered comparable.

Survey details of dataset used in the analysis are summarized in Table 1. The 2002 set has been used to construct the Cuvier’s beaked whale presence/absence model while the 2011 data set is used as independent evaluation set.

Table 1: Summary of the data collected during research cruises SIRENA 02 (calibration set) and SIRENA 11 (evaluation set). *Km conducted with 4 visual observers and Beaufort sea state less than 4.

SIRENA ‘02/ phase2 was conducted in July 2002 covering an area of about 11,000 km² (see Figure 2a). Oceanographic measurements were taken from the surface to a maximum depth of 1500 m, using a Conductivity, Temperature, Depth (CTD, Sea-Bird 911) and auxiliary sensors installed on a Rosette frame deployed from the N R/V Alliance at 21 pre-defined stations (Figure 1a). At each station, the following oceanographic data were collected: conductivity (mS/cm) to determine salinity, temperature (°C), density (kg/m³) derived from temperature and salinity, pressure to determine depth at which measurements are taken, fluorescence (ug/l) as a proxy for chlorophyll-a, dissolved oxygen (ml/l), turbidity an optical parameter collected with a turbidimeter as indicator of the scattered light due to suspended particles (FTU, formazine turbidity unit) and finally sound velocity (m/s) computed from salinity, temperature and depth measurements.

Figure 2a: SIRENA 02 dataset. Ziphius cavirostris sightings CTD stations and the corresponding 641 stations number are shown.

Cetacean visual sightings were collected using two big-eye verticalscale binoculars (Fujinon, 25*150, MTSX, field of view 2.75) mounted on the flying bridge of the NRV Alliance (at a height of 16 m above sea level) on the forward port and starboard sides while the ship was transiting at an average speed of 5-6 knots (i.e., about 10-12 km/h) in sea state conditions corresponding to a Beaufort scale lower than 4. Cetacean sightings were collected during daylight hours by a team of experienced observers that rotated through three observation positions (port, center, and starboard) and a recording position. Observers scanned the sea surface from 90° port to 90° starboard, where 0° is on the track line, with big-eye or regular binoculars (Fujinon 7*50 FMT/MT Field of view 7°30’). For every species encountered, time, bearing, radial distance, species, group size, behaviour, sighting cue, and swimming orientation (aspect) data were recorded. Environmental data, including sea state and effort status, were recorded every 30 minutes or more frequently if changes in conditions occurred. All the visual sighting data were recorded through dedicated data logging software. In the Genoa Canyon area a total of about 323 kilometres (174 nm) of track were surveyed under positive conditions (i.e., visual team on effort and Beaufort sea state lower than 4) and 24 sightings of Cuvier’s beaked whale were made.

The SIRENA11/phase A trial was conducted in 2011 during the time period and in the same area (Figure 2b and Table 1) investigated in 2002 and by used the same methods and protocols employed during SIRENA02 trial. Line transect surveys were conducted using passive acoustic and visual methods to determine the presence and absence of cetaceans. However, in the study presented here, only the visual data collected were employed in the validation model process to allow a comparison with the sea trial conducted in 2002. A total of about 870 km (470 nm) were covered in positive conditions resulting in 14 beaked whale sightings (i.e., 13 Cuvier’s beaked whales, and 1 undetermined beaked whale).

Figure 2b: SIRENA 11 dataset. Ziphius cavirostris sightings CTD stations and the corresponding

During the same time period, oceanographic data were collected at 18 stations for the same variables that were collected in 2002.

Oceanographic measurements collected during the two research campaigns were calibrated and processed by NURC. The variables considered in this study were conductivity, temperature, salinity, density, dissolved oxygen, fluorescence, turbidity and sound velocity. Supplementary Figures S1 and S2 in Annex A shows the vertical profile of each of these variables as a function of depth for both SIRENA cruises. To better understand the spatial pattern of the studied variables a series of plots, cross-sectional to the canyon area, were developed following the orientation northeast–southwest of canyons main axis. Additional east–west section plots were produced (see Supplementary Figures S3-S7 for oxygen section plots only). For the SIRENA 02 calibration dataset 39 water column statistics were computed, station by station, based on the 8 oceanographic variables available: (conductivity (mS/cm), salinity (PSU), density (kg/m3), temperature (C°), fluorescence (μg/l), dissolved oxygen (ml/l), turbidity (FTU) and sound velocity (m/s)). Table 3 summarizes the oceanographic variables derived from the CTD data used in the analysis. Minimum and maximum values (Min, Max) and the corresponding depths (Depth min, Depth max) were computed for every variable. Concerning the turbidity profiles the average (mean), median and standard deviation (SD) values were calculated. Gradients were computed as the ratio between the difference of the maximum and minimum values and the difference between the corresponding depths [i.e., (Max-Min)/(Depth Max-Depth Min)]. In addition, for the dissolved oxygen, the upper gradient (Ox Up gradient) was computed as the ratio between the difference of maximum and surface dissolved oxygen values and the depth of the maximum dissolved oxygen [i.e., (Max Ox-Ox at the Surface/Depth max)]. Fluorescence at 200 m was selected as the minimum measured fluorescence value. The depth of the 13.8°C isotherm was selected as the base of the thermocline. A Principal Component Analysis (PCA) was applied to the 39 parameters dataset to reduce the dimensionality of the CTD data matrix, and to select the most representative features of the water column profiles. For a better interpretation of the reduced components, the principal components previously extracted, were rotated applying a Varimax rotation criterion (i.e., the PCA solution was turned into a Factor Analysis solution), maintaining their orthogonality (i.e., the factors being still uncorrelated). By looking at the factor loadings matrix (i.e., the correlations of the original variables with the extracted factors) the most meaningful parameters within each factors were identified. Parameters lying on the same factor share the same information. The number of factors to retain was chosen on the basis of the “eigenvalue higher than 1” criterion (i.e., only the factors explaining more than the variance of one of the original variables were retained).

The Genoa canyon area was divided into 308 cells of about 5 × 4 km by means of GIS Spatial Analyst tool. Visual effort was evaluated in terms of kilometres of track line per cell unit. Only the sightings reported on-effort and in favourable conditions were considered. Encounter rates were then calculated for each cell unit as the number of sightings per kilometre surveyed under favourable conditions.

To test whether cells were spatially auto-correlated, the Moran’s Index of the cell encounter rates was computed. According to Legendre et al. [42] the presence of spatial correlation may bias tests of significance for correlation or regression, when both the animal response variable and the environmental predictor are spatially autocorrelated. Moran’s Index provided an assessment that ensured the chosen grid size was large enough to have the cell units encounter rates not spatially correlated. The factors extracted through the PCA/FA analysis of the CTD data were used to describe the hydrological characteristics of the cells. Factor values, available as point observations at the stations, were interpolated by using an Inverse Distance Weighted algorithm (IDW²) [43] allowing attribution of the CTD factor values to every cell unit. The cell mean values of every CTD factor, were used as covariates in a subsequent regression analysis.

Regression analysis

Binary logistic regression analysis [44,45] was employed to correlate the beaked whales presence/absence data to the oceanographic predictors. Initially the cell means, calculated for every CTD factor were used as covariates. Afterwards, a new regression analysis was carried out, substituting the factor values with the most representative (i.e., the variables with the highest factor loadings) CTD parameters. To balance the number of absence observations, presence data were weighted on the basis of encounter rates, according to Azzellino et al. [24,25]. In order to select the best set of predictors, a forward stepwise method based on the significance of the Wald statistic [46] was used. Each predictor was tested for entry into the model one by one, based on the significance level of the Wald statistic. After each entry, variables that were already in the model were also tested for possible removal. The variable with the largest probability greater than the specified threshold value was removed, and the model re-estimated; the procedure stopped when no more variables met the entry or removal criteria or when the current model was the same as the previous. The classification performance of the logistic models was evaluated in terms of confusion matrices [47]. The threshold value (i.e., cut-off) for the presence classification was set at a 0.5 probability.

Model validation: Sirena 11 dataset

Predictions for the Genoa canyon area were made using the habitat model developed on the 2002 data set (evaluation set). Model validation was then conducted using as independent evaluation set the data collected in the 2011 campaign and by evaluating the agreement between predictions and actual observations.

A 2 × 2 cross tabulation was created between observed and predicted presence/absence data (i.e., the confusion matrix according to Kohavi and Provost [47]). The confusion matrix can be used to describe predictive performance of models [47,48]. The sensitivity (i.e., the true positive fraction, the proportion of the positive observations in agreement with the presence predictions over the total positive observations) was used to evaluate the accuracy and the reliability of the models. For the validation set a more conservative threshold value, set at a 0.75 probability, was assumed for the presence prediction.

PCA/FA of CTD data.

As shown in Table 2, the 91% of the original variance was summarized by 8 factors, retained on the basis of the “eigenvalue higher than 1” criterion.

Table 2: Initial eigenvalues, variance explained, and cumulative variance explained for the principal component analysis (PCA) solution and the factor analysis (FA) rotated solution.

These 8 factors were interpreted on the basis of the factor loadings (Table 3) as follows:

Table 3:Factor analysis matrix with the variables used in the analysis. Each Factor loading represents the correlation between the original variables and the factor. Higher correlation coefficients are presented in bold.

Factor 1 (F1) is positively correlated with the maximum values of density, salinity, and with depth of the maximum salinity and minimum dissolved oxygen. F1 is also inversely correlated with the minimum measured fluorescence value (at a depth of about 200 m), and with the gradients of density and salinity. It accounts for 25.78% of the total variance.

Factor 2 (F2) is directly correlated with the minimum values of conductivity, sound velocity and temperature (13.8 °C) and with the depth of the maximum fluorescence. F2 is also inversely correlated with maximum values of dissolved oxygen and fluorescence. F2 accounts for 20.9% of the total variance.

Factor 3 (F3) is correlated with the maximum values of temperature, conductivity and sound velocity, close to the surface. F3 accounts for 15.8% of the total variance.

Factor 4 (F4) is the directly correlated with the mean, median and SD values of the turbidity. F4 accounts for 7.8% of the total variance.

Factor 5 (F5) is positively correlated with minimum salinity and the depth of the minimum salinity. It accounts for 6% of the total variance.

Factor 6 (F6) is positively correlated with dissolved oxygen surface values and inversely correlated with the oxygen upper gradient. F6 accounts for 5.8% of the total variance.

Factor 7 (F7) is positively correlated with the depth of maximum oxygen and the depth of the minimum sound speed. F7 accounts for 4.7% of the variance.

Factor 8 (F8) positively correlated with the minimum turbidity values, F8 accounts for 4% of the variance.

As described, the data available from the 21 CTD stations were interpolated through a IDW algorithm to a grid of 486 cells of about 3 nautical miles (0.05°) of cell size. Moran’s Index tested on the Cuvier’s beaked whale encounter rate showed that cells were not spatially auto correlated (Moran coefficient: I: 0.133962; SD: 0.032712).

The 8 factors were evaluated as predictors through the logistic regression analysis. The cell mean values, calculated for every factor, were used as covariates in the models. The stepwise procedure selected “Factor 7” (i.e., water column oxygen feature) as the best predictor (see Table 4). As shown in Table 5 the model showed a good accuracy in predicting both presence and absence cells (72.2%).

Table 4: Results of the binary logistic regression analysis model in the calibration set (SIRENA 02). Presence/absence of Cuvier’s beaked whale were correlated with cell means, calculated for every factor.

The following statistics are shown: B: Unstandardized regression coefficient; S.E.: Standard error of B; Wald statistic for the included parameter; df: Degrees of freedom; P: Level of significance.

Table 5: Confusion matrix of the model using factors cell means as predictors (Table 4).*Zc01 indicate Cuvier’s Beaked whale presence (1) and absence (0) cells.

Afterwards, a new regression analysis was carried out by substituting the factor with the original variables having the highest factor loading on the selected factor. In this second regression, the depth (m) of the maximum oxygen, (factor loading: 0.879) and the depth (m) of the maximum sound speed (factor loading: 0.689) were used as covariates. The strongest predictor of Cuvier’s beaked whale presence was the depth of the maximum oxygen value (hereinafter called: “Depth Ox Max”) being directly correlated to the whale presence (Table 6). In this case the overall model accuracy was 61.1%, showing a higher accuracy for the presence cells (66.7%) than for the absence cells (55.6%) see Table 7. Basically, as shown in Figure 3, Cuvier’s beaked whale presence in the Genoa canyon area is associated with cells showing deeper maximum oxygen layer with respect to the absence cells where the maximum oxygen layer was shallower.

Table 6:Results of the binary logistic regression analysis model performed in calibration set (SIRENA 02) using depth of the maximum oxygen value (Depth Ox max) as covariate.

Table 7: Confusion matrix of the models performed using depth of the maximum oxygen (Depth Ox Max) as predictor (Table 6).*Zc01 indicate Cuvier’s Beaked whale presence (1) and absence (1) cells.

Figure 3: Spatial interpolation of depth (m) of maximum oxygen layer (Depth Ox Max) collected during SIRENA 02 sea trial (calibration set). Ziphius cavirostris sightings, CTD stations and 1000 (light blue), 1400 (dark blue), 2000 (black) m depth contours are shown.

For model validation, the maximum oxygen depths were calculated station by station from the CTD data collected during the SIRENA 11/ phase A campaign and interpolated by using the IDW algorithm to the same grid that was used for the 2002 analysis.

Presence/absence predictions for the Genoa canyon area were made using the model based on the “Depth Ox Max” predictor developed on the 2002 data set (calibration set). Such predictions were overlaid to the Cuvier’s beaked whale observations collected during the Sirena 11/ phase A campaign (i.e., independent evaluation set) (Figure 4). The accuracy of the model was assessed through a cross-tabulation analysis (Tables 8). A Chi-square test showed a significant association between model prediction and the observations (χ2 = 4.267; df = 1; P-level < 0.05). The model showed good performance of classification, being slightly more accurate in classifying absence cells (87.5%) rather than presence cells (62.5%). This might be the effect of the low sample size (8 presence cells) used for the evaluation set.

Table 8: Cross-tabulation analysis of the model predictions versus the Cuvier’s beaked whale observations of the evaluation set (SIRENA 11).

Figure 4: Spatial interpolation of depth (m) of maximum oxygen layer (Depth Ox Max) collected during SIRENA 11 sea trial (evaluation set). Ziphius cavirostris sightings, CTD stations and 1000 (light blue), 1400 (dark blue), 2000 (black) m depth contours are shown are shown.

In this study, in situ oceanographic measurements collected in Genoa canyon area in concurrence with Cuvier’s beaked whale sightings, were analysed with the aim to identify the mechanisms that may influence the habitat selection of beaked whales at the canyon scale. This is the first time that measured oceanographic data were employed to model beaked whale presence in the Mediterranean Sea. The only other available example in the literature concerns the northern Adriatic Sea where CTD data were used to predict the habitat of common bottlenose dolphin, Tursiops truncatus [32]. In situ oceanographic measurements such as temperature and salinity were successfully employed to model beaked whale spatial density in the Pacific Ocean [10,49]. In that study the extent of the investigated spatial scale pointed out the effect of sea surface temperature and thermocline strength differences as boundary conditions. In the Ligurian Sea, the strong relationship of the species presence with the Genoa Canyon area [21,24-27] supports the hypothesis that this particular area has a special ecological meaning for Cuvier’s beaked whale, probably related to the prey availability within the canyon. In situ oceanographic measurements provided new elements for identifying the ecological mechanisms that can influence beaked whale habitat selection at canyon scale. The results of this study suggest that dynamic variables derived from oceanographic measurements (i.e., CTD data) may be employed as predictors of Cuvier’s beaked whale presence in Genoa canyon area. Particularly, the depth (m) of the maximum concentration of dissolved oxygen (Depth Ox Max) was found to be a significant predictor of beaked whale presence, providing good performance of presence/absence classification in the calibration dataset (SIRENA 02) and showing a good model accuracy for the evaluation set (SIRENA 11). The fact that model accuracies were comparable for both the calibration and evaluation data set proves the reliability of the predictor “Depth Ox Max” for the data collected 9 years apart. Higher probabilities of beaked whale presence are in fact directly associated with the areas showing the deepest maximum oxygen layer. Dissolved oxygen profiles for this region show that the maximum dissolved oxygen concentration occurs within the first 50 m of the water column, than it decreases with depth, the minimum concentration layer occurring around 400 m, and then showing a low and gradual increase in concentration up to the maximum recording depth (about 1200 m). It is widely known that beaked whales are deep diving species, performing deep foraging dives to feed on deep water food resources. Cuvier’s beaked whale diving profiles in the Ligurian Sea [50] showed foraging activity in mesopelagic to bathypelagic water depths (613–1297 m). The stomach contents of 3 Cuvier’s beaked whales stranded along the Ligurian coast consisted of digested mesopelagic cephalopod beaks principally of the Histioteuthidae family, specifically Histioteuthis reversa and H. bonnellii and other cephalopods species such as Octopoteuthis spp. (Octopoteuthidae), Galiteuthis armata (Cranchiidae), Chiroteuthis. veranii (Chiroteuthidae) e Ancistroteuhis lichtensteinii (Onychoteuthidae) [51]. From this perspective, the predictor “Depth Ox Max” selected in this study does not reflect the foraging depths generally used by the species, but may be reasonably considered a proxy for vertical exchanges in the water masses. In fact, dissolved oxygen in the ocean is known to be a useful tracer for tracking the movement of water-types and water-masses [52]. It is known that submarine canyons may modify flow, shelf-slope exchanges of water and material [53,54] and this combination may influence the transport of particulate organic matter that influence productivity. The degree of complexity of the sea bottom, relative position, size and morphology of the canyons as well as wind stress strength, affect the local circulation, altering current flows favouring the rise or sinking of the water masses [55,56]. In addition, canyons environments are interested by phenomena such as dense shelf water cascade that transport large amounts of water, sediment and organic matter, downslope in the canyon to the deep ocean, reshaping submarine canyon floors and affecting the deep-sea environment [57]. In the Genoa canyon area, a deeper maximum oxygen layer could suggests a sinking of productive surface waters that brings oxygenated waters deeper due to the motion of the water mass in the canyon (coastward downwelling) stimulated by the frontal area. This downwelling motion allows the fresh organic material produced near the surface to sink to depth, thus enriching the water column and enhancing the biological processes. It is widely accepted that vertical sinking of organic material is one of the main source of nutrients to the deep-sea floor [56]. The results obtained on seawater samples collected during SIRENA 02 [35] in the Genoa canyon, suggest that the transfer of organic matter to depth, driven by downwelling, may sustain intense microbial activity in the mesopelagic layer. Moreover, the same authors found increases in ectoenzymatic activity in the mesopelagic layer that was linked to the presence of active microbial assemblages which was thought to enhance the recycling potential of the food web. According to Misic and Fabiano [35] the productive events, that are not dependent on surface productivity processes, can be associated by increase of turbidity at depth. To investigate this aspect, the mean values of turbidity between the depths of 400 m and 1200 m (range of possible feeding habitat of beaked whale), were calculated in both calibration (SIRENA 02) and validation set (SIRENA 11). It was shown that the areas with higher presence probability (i.e., probability higher than 75% based on the Depth Ox max model predictions), are also characterised by a higher mean (Mann Whitney U = 2634; P-level < 0.1) and maximum (Mann Whitney U = 2401; P-level < 0.05) turbidity with respect to the lower probability areas (i.e., probability lower than 75% based on the Depth Ox max model predictions) as shown in Figure 5.

Figure 5: SIRENA 02, calibration set. Mean turbidity values (FTU) for beaked whale presence and absence cells (MW test U = 2643; Plevel < 0.1).

The same pattern was also observed in the evaluation set (Mann Whitney U = 9317; p-level < 0.05) (Figure 6). Thus, the depth of the maximum oxygen layer acts as a proxy for areas where vertical motions support the enhancing of biological processes at the different trophic levels. The trophic enrichment suggested by the increase of the turbidity in the canyon’s deep waters may be an attractor for beaked whale prey, such as the mesopelagic cephalopods, and be a key factor in the beaked whale habitat selection within the canyon.

Figure 6: SIRENA 11, evaluation set. Mean turbidity values (FTU) for beaked whale presence and absence cells (MW test U = 9317; Plevel < 0.05).

This study provided the unique opportunity to investigate how ocean dynamics might influence the habitat selection of Cuvier’s beaked whale in one Mediterranean key habitat, the Genoa Canyon area. In situ oceanographic measurements were employed as predictors to model beaked whale presence in the Canyon area. The depth of the maximum oxygen was found directly correlated with beaked whale presence. The relevance of such predictor was also validated based on an independent dataset which was generated 9 years later in the same area. We assume that depth of the maximum oxygen might be a proxy of downwelling regions where energy from the richer surface waters is transferred to the deep waters, enhancing the productivity in the meso- and bathy-pelagic environment where beaked whales feed on their prey. These results suggest that dynamic predictors may act as proxy of the macro-scale features that characterise beaked whale habitat. Further research should be dedicated to identify more accessible (e.g. remote sensed parameters) proxies of such dynamics, enabling the confirmation of these results and improving the identification of beaked whales habitat in different areas of the Mediterranean Sea.

The authors acknowledge and appreciate the support of the NRV “Alliance’s” Captains and crew for their passionate support during the Sirena/MED trials. Additionally, the authors would like to acknowledge the MMRM Project, all NURC (CMRE) personnel that worked at the project during the years and all the colleagues from different research organization who participated in the data-collection effort. The authors want to thank Dr. James Eckman, for the support and interest demonstrated about the topic of this study, Jerry Olen of SPAWAR Systems Center Pacific for his continued support and Philippa Dell that reviewed the paper for English language. This work was principally supported by the Office of Naval Research (Grant N00014-10-1-0533) and also partially supported by the EC FP7 Marie Curie Actions People, Contract PIRSES-GA-2011-295162–ENVICOP project (Environmentally Friendly Coastal Protection in a Changing Climate).

¹NURC, has recently changed name. Now is CMRE, Centre for Maritime Research and Experimentation.

²The IDW interpolator assumes that each input point has a local influence that diminishes with distance. So the points closer to the processing cell have greater weight than more distant points.