Assessing cetacean surveys throughout the Mediterranean Sea: a gap analysis in environmental space

Heterogeneous data collection in the marine environment has led to large gaps in our knowledge of marine species distributions. To fill these gaps, models calibrated on existing data may be used to predict species distributions in unsampled areas, given that available data are sufficiently representative. Our objective was to evaluate the feasibility of mapping cetacean densities across the entire Mediterranean Sea using models calibrated on available survey data and various environmental covariates. We aggregated 302,481 km of line transect survey effort conducted in the Mediterranean Sea within the past 20 years by many organisations. Survey coverage was highly heterogeneous geographically and seasonally: large data gaps were present in the eastern and southern Mediterranean and in non-summer months. We mapped the extent of interpolation versus extrapolation and the proportion of data nearby in environmental space when models calibrated on existing survey data were used for prediction across the entire Mediterranean Sea. Using model predictions to map cetacean densities in the eastern and southern Mediterranean, characterised by warmer, less productive waters, and more intense eddy activity, would lead to potentially unreliable extrapolations. We stress the need for systematic surveys of cetaceans in these environmentally unique Mediterranean waters, particularly in non-summer months.


Results
Spatiotemporal coverage of surveys. We aggregated line transect surveys conducted by 12 organisations, including universities, consultancy companies and non-governmental organisations, resulting in a total of 302,481 km of effort (Table 1, Fig. 2). Aerial surveys represented 149,225 km of effort and shipboard surveys represented 153,256 km of effort.
Survey effort was concentrated in the northwestern Mediterranean and comparatively patchy in the eastern and southern Mediterranean (Fig. 3a). The Algero-Provençal basin was the Mediterranean subregion that received the largest amount of survey effort, followed by the Tyrrhenian Sea/eastern Ligurian Sea and the Alborán Sea/Strait of Gibraltar (Table 2). In comparison, the Levantine Sea and Aegean Sea subregions received only 2.1% of the effort each.
Aerial surveys were conducted off Spain west of the Balearic Islands, in the Pelagos Sanctuary (a marine protected area for marine mammals in the Ligurian Sea established by international treaty 34 ), in the Tyrrhenian Sea, in the Strait of Sicily (especially around the Maltese Islands), in the northern Ionian Sea, and in the Adriatic Sea (Fig. 3b). Shipboard surveys were conducted primarily in the northern Alborán Sea, in the western Ligurian Sea, in the eastern Ionian Sea, and around the Maltase islands. Shipboard survey effort was very patchy in the rest of the Mediterranean Sea (Fig. 3c).
In terms of inter-annual survey coverage, two peaks were apparent: a small peak in the early 2000s and a larger peak between 2009 and 2013 ( Fig. 4a)  was exceptionally low in 2015. The northern Alborán Sea, the northern Gulf of Lions and Ligurian Sea, the area around Maltese Islands, and the Hellenic Trench were covered by surveys in most years, reflecting long-term monitoring efforts ( Supplementary Fig. S1).  In terms of intra-annual survey coverage, the amount of survey effort was largest in July and August and lowest in February and April (Fig. 4b). Most surveys were conducted from May to September. Survey effort was spread across the Mediterranean Sea from May to September, but restricted to the northwestern and central Mediterranean the rest of the year ( Supplementary Fig. S2).

Surveying entities Platform Surveyed years Surveyed subregion 1 Total effort (km) References
Environmental coverage of surveys. Univariate environmental coverage. The surveys spanned broad ranges of values for all considered static covariates (depth, slope, distance to seamounts, and distance to canyons), resulting in limited univariate extrapolation throughout the Mediterranean Sea (Table 3) (Supplementary Figs S3  and S4).
The surveys spanned narrower ranges of values for the considered dynamic covariates (sea surface temperature, primary productivity, and eddy kinetic energy), resulting in a larger univariate extrapolation (especially in non-summer months) ( Table 3).
The extent of univariate extrapolation would be largest with the model including sea surface temperature only. In winter and spring, extrapolation would be needed to predict cetacean densities in warmer waters of the southeastern Mediterranean and in colder waters of the northern Mediterranean (Gulf of Lions, Adriatic Sea, and northern Aegean Sea) (Fig. 5a, Supplementary Figs S5 and S6). In summer and autumn, extrapolation would occur in warmer waters of the Levantine Sea and of the Tunisian Plateau (but with a reduced extrapolation compared to winter and spring).
Extrapolation in more productive coastal waters would occur year-round on the Tunisian Plateau and in non-summer months in the northern Adriatic Sea and off Egypt. In spring and autumn, extrapolation would occur in the less productive offshore waters of the southeastern Mediterranean (Fig. 5b, Supplementary Figs S7 and S8).
The extent of univariate extrapolation would be lowest with the model including eddy kinetic energy only. Extrapolation would occur in waters characterised by more intense eddy activity in the Alborán Sea and Levantine Sea, mostly in winter and spring (Fig. 5c, Supplementary Figs S9 and S10).
Multivariate environmental coverage. As expected, the extent of extrapolation was larger with combinations of covariates than with individual covariates (Table 3). For example, the extrapolation extent with the combination of sea surface temperature, primary productivity, and eddy kinetic energy was larger than the extrapolation extents with these three covariates individually.
For the model including all four static covariates, extrapolation would be limited to 3.7% of the Mediterranean Sea (Table 3). Extrapolation would mostly occur in offshore waters of the central Mediterranean and off Egypt (Fig. 6a). The proportions of data nearby in the multivariate space defined by static covariates were lowest on the Tunisian Plateau, in the northern Adriatic Sea, and in offshore waters throughout the Mediterranean Sea (Fig. 6b). Differences between extrapolation metrics (the extent of extrapolation and the proportion of data nearby) are presented in the Methods and discussed further below. For the model including all three dynamic covariates, extrapolation would occur in 49.9% of the Mediterranean Sea on average (minimum: 8.1% in September; maximum: 55.1% in February) ( Table 3). Extrapolation would be widespread in the eastern Mediterranean and in parts of the western Mediterranean (e.g., the Alborán Sea and the Ligurian Sea), especially in winter and spring (Fig. 7a, Supplementary Fig. S11).  In summer, extrapolation would be mostly limited to the Levantine Sea. The proportions of data nearby in the multivariate space defined by the three dynamic covariates were low in most of the eastern Mediterranean and in the northern and southern parts of the western Mediterranean (Fig. 7b, Supplementary Fig. S12). In autumn, higher proportions of data nearby appeared in a mid-latitude band. Some areas such as the Alborán Sea and the northern Adriatic Sea had low proportions of data nearby year-round.
For the model including all static and dynamic covariates, extrapolation would increase to 80.1% of the Mediterranean Sea on average (minimum: 55.5% in July; maximum: 96.5% in February) ( Table 3). Extrapolation would be widespread in all seasons except in the western Mediterranean in summer (Fig. 8a, Supplementary  Fig. S13). Overall, the proportions of data nearby in the multivariate space defined by all static and dynamic covariates were lower compared to the proportions of data nearby in the multivariate space defined by dynamic covariates only (as seen by the comparison of colour scales between Figs 8b and 7b). In winter and spring, the proportions of data nearby were low throughout the southeastern Mediterranean and in offshore parts of the western basin (Fig. 8b, Supplementary Fig. S14). In summer, the proportions of data nearby were overall higher in coastal areas, particularly to the north. In autumn, the proportions of data nearby were low in most of the eastern Mediterranean and in the northwestern Mediterranean.

Discussion
We aggregated over 300,000 km of cetacean line transect surveys suitable for estimating abundance that were conducted in the Mediterranean Sea within the past 20 years by numerous organisations. Survey coverage was heterogeneous geographically and seasonally: large data gaps were present in the eastern and southern Mediterranean in all seasons and elsewhere in non-summer months. Survey coverage was also heterogeneous in environmental space. Surveys covered a number of static bathymetric features (e.g., coastal areas, continental slopes, submarine canyons and seamounts), but not the full range of dynamic oceanographic conditions found within the   Table 3. Spatial extent of extrapolation (i.e., the percentage of cells of the study area where extrapolation occurred) with single covariates and combinations of covariates. For dynamic covariates, the mean extent of extrapolation averaged over the 12 month period is provided, followed by the minimum and maximum monthly extents in parentheses. SST: sea surface temperature; PP: primary productivity; EKE: eddy kinetic energy. Mediterranean Sea. In particular, Mediterranean waters characterised by comparatively warmer temperatures, lower productivity and higher eddy activity were poorly surveyed for cetaceans. This raised the prospect that cetacean density models fitted to environmental covariates would have to be extrapolated in order to provide predictions for the entire Mediterranean Sea.
To assess extrapolation in environmental space, we mapped two metrics: the extent of extrapolation versus interpolation-a binary metric indicating whether a point is within the convex hull enclosing the sampled environmental space-and the proportion of data nearby in multivariate environmental space-a continuous metric derived from Gower's distance. The first metric distinguishes predictions in unsampled geographic space that are within sampled environmental space (interpolation) from predictions in unsampled geographic and  environmental space (extrapolation). It can be used to restrict predictions to the sampled environmental space by overlaying a mask on the prediction map or by flagging extrapolations in unsampled environmental space 13,35 ; the standard errors of such extrapolations can then potentially be inflated until new data are acquired. The second metric provides a quantitative measure of the reliability of extrapolation in multivariate environmental space. It can be used as a measure of prediction uncertainty when maps of species distributions are needed for an entire region (e.g., to support management decisions). This metric is useful for differentiating predictions that are just outside the sampled environmental space from those that are far outside (and thus less reliable). It also identifies predictions that occur in sparsely sampled regions of environmental space despite being classified as interpolations by falling within sampled environmental space. For example, this situation was observed in the northern Adriatic Sea where low proportions of data nearby would occur although interpolation is predominant (see Fig. 6). Using the binary metric that distinguishes interpolation from extrapolation as the only measure of the reliability of predictions would lead to an overly high degree of confidence in predictions in the northern Adriatic Sea.
Mesgaran et al. 15 proposed a modification of the Mahalanobis distance to assess deviations between the prediction and the calibration dataset using the mean and linear correlations between covariates 15 . We did not follow their approach because the covariates we selected a priori were not highly correlated (see Methods). Conn et al. 36 assessed extrapolation using a model-based convex hull derived from a generalization of Cook's independent variable hull. Because our study objective was to conduct a gap analysis of survey effort in environmental space (rather than to fit models), we adopted a model-independent approach. We found that the extent of extrapolation was larger for combinations of covariates than for individual covariates. This is a common issue when predicting SDMs beyond the scope of the data: as more covariates are included in a model, new individual covariate values and new combinations of covariates appear, rapidly increasing the extent of extrapolation 13 . For example, in a study of small pelagic fish in the Gulf of Lions, models that included single covariates yielded less than 5% extrapolation while models that included five covariates exceeded 50% extrapolation 37 . For the Mediterranean cetacean surveys, a model that included all seven covariates resulted in 80% extrapolation; this model would be overly complex (in terms of having too many parameters) for reliably predicting cetacean densities at the regional scale. To predict species distributions or densities beyond the scope of the calibration data using SDMs, we recommend choosing models with a small number of covariates to limit the extent of extrapolation in environmental space.
Our analysis showed that the extent of extrapolation would be lower with static covariates than with dynamic covariates. Advantages of static covariates for modelling cetacean densities include generally low measurement errors and easy acquisition and processing leading to prompt model fitting. However, static covariates are rather indirectly related to cetacean densities. For example, steep seafloor slope may aggregate cetaceans via oceanographic processes that enhance prey availability and accessibility, assuming cetaceans are involved in foraging activities. Indirect ecological relationships are known to transfer poorly to new geographic regions because ecological processes underlying cetacean distributions may differ between surveyed and unsurveyed regions 38,39 .
By design, a gap analysis does not require the observation data (from which the response variable is derived) because it aims at assessing how representative sampling effort is of the environmental space where predictions are to be made. However, for the gap analysis to be useful, the environmental covariates to be tested should be the ones that are most likely to be selected when fitting the models. This can be done a priori (as in this study) by using sets of covariates that are ecologically-relevant and likely to be selected in future cetacean density models. This could also be done a posteriori for a given case study by using the set of covariates that has been actually selected in the model selection process. In this latter case, the response variable would be needed, not for the gap analysis itself, but to help identify the set of most relevant environmental covariates to be tested. When it comes to model development, it is important to keep in mind that (1) sample size of species observations affects model accuracy and precision 40 and (2) the selected model strongly influences predictions, particularly when predicting outside the sampled range of the covariate data 41,42 .
Our study showed that if cetacean density models were calibrated on existing line transect surveys and predicted across the Mediterranean Sea, widespread extrapolations in environmental space would occur, leading to potentially unreliable predictions. However, it should be noted that the surveys available for this exercise were limited to those compatible with distance sampling, which is used to estimate absolute abundance. For certain management applications, relative abundance predictions will suffice. There are a number of additional surveys in the Mediterranean Sea that were not aimed at estimating absolute abundance with distance sampling but that followed similar methodology e.g., 43,44 and that could be included in models of relative abundance. Including such surveys would likely reduce some of the gaps in environmental space.

Conclusion and Perspectives
Systematic survey efforts for cetaceans have previously been recognized as heterogeneous across the Mediterranean Sea 18,45,46 , but this heterogeneity was never quantified. By aggregating surveys conducted throughout the Mediterranean Sea, we were able to identify gaps in the spatiotemporal and environmental coverage of cetacean surveys for the first time.
We strongly recommend that future surveys be conducted in the eastern and southern Mediterranean where waters are characterised by warmer sea surface temperature, lower productivity, and higher eddy kinetic energy. These environmental conditions are found nowhere else in surveyed regions; they are unique to the eastern and southern Mediterranean and were under-represented in our dataset. Conducting surveys there would be particularly valuable in non-summer months because oceanographic conditions differ the most from those of other areas surveyed in those months. We also recommend additional survey effort in the northern Mediterranean in non-summer months, for example in the Gulf of Lions in April, which is characterised by colder waters than other areas surveyed in that month.
A Mediterranean-wide survey is scheduled to take place in the summer 2018 thanks to the long-standing initiative of ACCOBAMS (Agreement on the Conservation of Cetaceans in the Black Sea, Mediterranean Sea and Contiguous Atlantic Area). This regional survey will be extremely helpful to fill existing data gaps by collecting data synoptically with a consistent methodology. Unfortunately, this survey will not be sufficient to fill data gaps in non-summer months in which environmental conditions (and by extension cetacean distributions) are widely different. Achieving a representative sample of environmental space is critical for surveys intended for habitat-based density modelling of cetaceans. An analysis of environmental space early in the survey planning process-for example assessing the environmental representativeness of a set of potential survey designs-can help ensure model-based predictions of cetacean densities are derived with limited extrapolation for the study region.
The natural next step of this study is to incorporate cetacean observations obtained from the aggregated surveys to build models and predict cetacean densities across the Mediterranean Sea. The resulting density maps could be related to maps of human activities such as ship traffic and underwater noise, which are particularly harmful to cetaceans in the Mediterranean Sea 30,47 , in order to identify unsurveyed areas where populations may be more at risk. Predicted density maps could also be used for delineating areas of potentially high cetacean densities in poorly surveyed regions of the Mediterranean Sea, and thus be utilised as part of the Important Marine Mammal Area (IMMA) process of International Union for Conservation of Nature.

Study area. The Mediterranean Sea is a semi-enclosed water body connected to the Atlantic Ocean by the
Strait of Gibraltar, to the Black Sea by the Bosphorus, and since 1869 to the Red Sea by the Suez Canal. It is divided into a western and an eastern basin by a central ridge between Sicily and the Tunisian-Libyan coast. The Mediterranean Sea is mainly characterised by narrow continental shelves, steep slopes and extensive abyssal plains. It includes a variety of submarine canyons, mostly located along the continental slopes in the north. It also includes approximately one hundred seamounts ( Supplementary Fig. S15), known to affect the distribution of pelagic species, including cetaceans [48][49][50] . The Mediterranean is an oligotrophic sea characterised by salty and nutrient-poor waters 51 .
Circulation in the Mediterranean Sea is mainly driven by water flow through the Strait of Gibraltar, freshwater inputs from the main rivers (Nile, and to a lesser extent, Po, Rhone and Ebro), wind stress, and thermohaline and topographic features 52 . Atlantic Surface Water flows into the Mediterranean Sea through the Strait of Gibraltar and circulates in a cyclonic (counterclockwise) direction (Fig. 1). Water flow along the southern coasts generates short-lived mesoscale anticyclonic eddies (e.g., the eddy field off Algeria). To the north, water flow creates persistent cyclonic gyres (e.g., the Lions gyre) associated with upwelling of nutrient-rich waters that result in enhanced primary productivity 52 . As it moves eastward, surface water evaporates and becomes saltier, warmer and poorer in nutrients, resulting in a gradual decline in phytoplankton biomass and productivity from west to east 16,53  Phytoplankton biomass and primary production have marked seasonal cycles in the Mediterranean Sea 54,55 . Phytoplankton blooms are primarily initiated in winter and spring by wind stress, causing mixing and nutrient uplift to surface layers. Uplift of nutrients also occurs at cyclonic eddies. While blooms are markedly seasonal and intense in the northwestern basin (e.g., in the Gulf of Lions), they are often sporadic and subject to significant inter-annual variability in the eastern basin. Stratification occurs in summer, resulting in a lower and more homogeneous phytoplankton biomass across the Mediterranean Sea.
The spatial extent for this study corresponded to the entire Mediterranean Sea, excluding enclosed lagoons and estuaries (surface area: 2.5 million km 2 ). The study area was projected to a custom Lambert Azimuthal Equal Area projection in ArcGIS 10.2.2 and gridded into 10 × 10 km cells for analysis, representing a compromise between the various resolutions of environmental covariates (see Supplementary Tables S1 and S2). The Mediterranean study area was divided into 8 distinct ecological subregions following previous studies 18,56 (Fig. 2).
Line transect surveys. We aggregated data from visual shipboard and aerial line transect surveys that followed distance sampling methodology, i.e., allowed the calculation of perpendicular distances to observed groups of cetaceans and thus the estimation of absolute abundances 57 . We considered both occasional large-scale surveys and recurrent small-scale surveys aimed at abundance estimation (homogeneous coverage probability was not a requirement for this analysis). We included only line transect surveys conducted after October 1997 to ensure the full set of environmental covariates was available (primary productivity data was available from October 1997; see Supplementary Table S2), and until April 2016.
Spatiotemporal coverage of surveys. For each survey, we created in ArcGIS a shapefile of on-effort tracklines with associated latitude, longitude and date. We merged shapefiles of all survey tracklines and added them to a file geodatabase. We then intersected and spatially joined the tracklines to the 10 × 10 km grid of the study area. To examine the spatiotemporal coverage of surveys, we summed survey effort per grid cell in the entire study area and created subsets of surveys per month, year, and platform (ship or aircraft).

Environmental covariates.
We examined the environmental coverage of surveys with respect to four static covariates (depth, slope, distance to canyons, and distance to seamounts) and three dynamic covariates (sea surface temperature, primary productivity, and eddy kinetic energy) known to be important for cetaceans and likely to be included in future habitat-based density models (Supplementary Tables S1 and S2).
We created rasters of static and dynamic covariates (data sources detailed in Supplementary Tables S1 and S2) and resampled them to the 10 × 10 km grid of the study area (using nearest neighbour interpolation implemented in the Marine Geospatial Ecology Tools software 58 ) in ArcGIS.
We used monthly contemporaneous and climatological resolutions for dynamic covariates 59 . For surveyed grid cells (i.e., representing the "available" dataset; see below), we used the month and year when the survey took place (contemporaneous resolution). For all grid cells of the study area (i.e., representing the "prediction" dataset; see below), we used a monthly climatological average describing long-term oceanographic conditions in the Mediterranean Sea. Doing so allowed us to compare surveyed oceanographic conditions with long-term oceanographic averages (rather than oceanographic conditions of a particular year). To obtain monthly climatological rasters, we averaged monthly rasters per grid cell from 1997-2016. These climatological means smooth out the inter-annual variability of oceanographic processes while maintaining the strong intra-annual variability which characterizes the Mediterranean Sea 52,54,55 .
Environmental coverage of surveys. To analyse the environmental coverage of surveys, we examined two multivariate environmental metrics: the extent of extrapolation versus interpolation and the proportion of prediction points near available data points.
To evaluate the extent of extrapolation versus interpolation if models calibrated on the available survey data were used to predict cetacean densities in the entire Mediterranean Sea, we used the convex hull approach 37,60,61 . The convex hull of a set of points is defined as the smallest convex set that contains these points. Predictions inside the convex hull are interpolations while predictions outside the convex hull are extrapolations. The convex hull of a single covariate is the interval between the minimum and maximum data points (i.e., a univariate environmental envelope). The convex hull of two covariates is a polygon with vertices at the extreme points of the data points (i.e., a bivariate environmental envelope). A convex hull can be defined for any number of covariates, although visualization becomes difficult for more than two covariates and computational power starts to limit calculations for more than ten covariates 61 . As stressed by Authier et al. (2017), the assessment of the convex hull does not require model fitting (only effort and environmental covariate data are required) 37 .
Here, environmental characteristics of the surveyed grid cells during the month and year of each survey represented the "available" dataset and monthly climatological environmental conditions of all grid cells of the study area represented the "prediction" dataset. We considered the seven individual covariates enumerated above and the three following combinations of covariates: • all static covariates: depth, slope, distance to seamounts, and distance to canyons; • all dynamic covariates: sea surface temperature, primary productivity, and eddy kinetic energy; • all static and dynamic covariates: depth, slope, distance to seamounts, distance to canyons, sea surface temperature, primary productivity, and eddy kinetic energy.
An examination of Spearman's rank correlation coefficients 62 showed that these covariates were not highly correlated (Spearman's rho <|0.6|) in both available and prediction datasets. To further quantify extrapolation in environmental space, we calculated Gower's distances (G 2 ) for the three combinations of covariates listed above. As for the convex hull, the calculation of Gower's distances does not require model fitting. The Gower's non-parametric distance between points i and j belonging to the available and prediction datasets, respectively, is defined as the average absolute distance between the values of these two points in each dimension, divided by the range of the data [61][62][63] . With K environmental dimensions, the Gower's distance formula is as follows: where r k is the difference between the largest and the smallest values of the available dataset for the k th covariate. Because of range standardization, covariates have equal contribution to the Gower's distance. We considered a prediction point to be in the neighbourhood of available data points if it was situated within a radius of one geometric mean Gower's distance of all pairs of available data points 61 . We stress that this definition of neighbourhood relates to the multivariate environmental space defined by the considered environmental covariates, not geographic space. The larger the proportion of prediction points near available data points in this multivariate environmental space ("proportion of data nearby" for brevity), the more reliable the extrapolation. We used the WhatIf package (version 1.5-6) 64 in R (version R-3.1.1) 65 to calculate convex hulls and Gower's distances. We projected the extent of extrapolation versus interpolation and the proportion of data nearby to geographic space for mapping.
Data availability statement. The aggregated survey dataset and R code to reproduce the gap analysis are made available via the Dryad repository (https://doi.org/10.5061/dryad.4pd33).