Species data were collected during the EU-funded MEDITS project (Bertrand et al. 2002) carried out from spring to summer (April to September depending on the country) from 2006 to 2015. The MEDITS project uses a stratified sampling design based on depth (five bathymetric strata: 10-50 metres, 51-100 m, 101-200 m, 201-500 m and 501-700 m) and GSAs. Sampling stations were placed randomly within each stratum at the beginning of the project. In all subsequent years sampling was performed in similar locations. Specifically, this study concerns 13 GSAs including four different Mediterranean countries: namely, Spain (GSAs 01 - Alboran Sea; 02 - Alboran Island; 05- Balearic Islands; and 06 - Levant), France (GSAs 07 - Gulf of Lions; and 08 - Corsica), Italy (GSAs 09 - Ligurian Sea, northern and central Tyrrhenian Sea; 10 - central and southern Tyrrhenian Sea; 11 - Sardinia and 19 - western Ionian Sea) and Greece (GSAs 20 - eastern Ionian Sea; 22 - Aegean Sea; and 23 - Crete). The studied time series covers the period 2006 to 2015 as the sampling methodology is more homogeneous between those years. However, it is worth noting that the Greek dataset shows some important gaps, being available only for 2006, 2008 and 2014 because of problems in Data Collection Framework (DCF) implementation. Nevertheless, they were included in the study because of the quantitative and qualitative importance of that area. The model was balanced to deal with such a lack of some data in the eastern area.

MEDITS data of occurrence (presence/absence), biomass (kg km^–2) and abundance (ind. km^–2) of the two anglerfish species were used as response variables in the statistical models. Biomass data were log-transformed for the statistical models to ensure a normal distribution and to reduce the weight of extreme values.

Anglerfish distributions were estimated using five possible predictors: SBT, SBS, bathymetry, seabed slope and distance to coast. Monthly averages (corresponding to the sampling) of SBS (expressed in practical salinity units) and SBT (expressed in °C) covering the entire study period (2006-2015) were extracted from the EMIS European Seas project (http://mcc.jrc.ec.europa.eu/emis/) with a spatial resolution of 4 km. These predictors were included in the analysis because SBT and SBS are strongly related to marine system productivity affecting nutrient availability, metabolic rates and water stratification (Pennino et al. 2013).

Bathymetric features were derived from the MARSPEC database (http://www.marspec.org) (Sbrocco et al. 2013) with a spatial resolution of 1 km. Bathymetry and distance to coast were considered for analysis as they are among the main drivers of species distribution and have been previously reported as key factors for identifying spatial patterns of many invertebrate and fish demersal communities (Pennino et al. 2013, Lauria et al. 2015, Paradinas et al. 2015). In addition, the slope of the seabed describes its morphology. In particular, low slope values correspond to a flat ocean bottom (usually areas of sediment deposition), while higher values indicate potential rocky ledges (Lauria et al. 2015).

All variables were aggregated at 4×4 km cells in order to have the same spatial resolution and were explored for collinearity, outliers and missing data before being used in the modelling (Zuur et al. 2010). In addition, variables correlation was tested with the Pearson correlation index. Variables were not highly correlated (r<0.5), so all were considered in the further analysis.

In order to identify and predict the preferential habitats of the anglerfish species in relation to environmental variables in the Mediterranean Sea, a spatio-temporal model was applied. In particular, exploratory analysis highlighted that anglerfish species biomass and abundance data have three main issues: (i) a large proportion of zeros, (ii) a non-linear relationship between the response variable and explanatory variables and (iii) spatial autocorrelation. In this study, we implemented a Delta model approach with RAC-BRT to handle with all these issues at the same time (Elith et al. 2008, Crase et al. 2012, Escalle et al. 2016). The Delta approach deals with high numbers of zeros using two modelling stages. In the first one, presence/absence data were modelled using a binomial distribution in order to obtain a prediction of presence probability of the studied species; in the second, biomass and abundance data were modelled only in areas where species were expected to be present. Here, the biomass data were modelled using a Gaussian distribution and the abundance data using a Poisson distribution. For both stages, all environmental variables and possible interaction terms were included, in addition to a yearly temporal factor to account for annual variability. Model selection was performed using both forward and backward approaches and looking at the total explained deviance of each variable. In addition, in order to take into account that some country has more data than another, a weighting factor was used to indicate the number of years for which data were available.

The modelling technique implemented in the Delta approach was that of BRT because they deal with non-linear and non-monotonic relationships between response and explanatory variables. Indeed, BRT combines regression trees and boosting methods to fit complex non-linear relationships between predictors and the response variable (Elith et al. 2008). BRT parameters were selected to optimize the model using the caret (Kuhn 2008) and dismo packages (Hijmans et al. 2014) of the R software (R Core Team 2017). In particular, the optimal number of boosting trees was assessed with the gbm.step function, while the tree complexity of the model was fixed at 2 and the learning rate at 0.01.

Finally, the RAC technique was implemented to deal with the spatial correlation issue during both stages of the BRT Delta approach. Specifically, the spatial autocorrelation was accounted for by adding the autocovariate term to the final selected BRT (after model selection), which represents the influence of neighbour observations on the response variable at a specific location (Escalle et al. 2016). In addition, the Moran index and the Moran statistical test were used to verify the spatial correlation of the data using the spdep R package (Bivand 2013).

Predictive ability of binomial BRT models was assessed using a cross-validation procedure with the PresenceAbsence R package (Freeman 2012). In particular, the dataset was randomly split into two main subsets: a training dataset including 75% of the data and a validation dataset containing the remaining 25%. The training dataset was used to assess the relationship between occurrence data and the environmental variables, while the quality of predictions was evaluated using the validation dataset. This calibration-validation procedure was repeated ten times for each species and the results were averaged over the random subsets. As is common in this procedure, a confusion matrix that records the number of true/false positive/negative cases predicted by the model was generated. From this matrix, the area under the receiver operator curve (AUC) (Pearce and Ferrier 2000) and the true skill statistics (TSS) (Allouche et al. 2006) measures were computed to evaluate the prediction of the models. In particular, the AUC measures the ability of a model to discriminate between sites where a species is present and those where it is absent. AUC ranges from 0 to 1, with values below 0.6 indicating a performance no better than random, values between 0.7 and 0.9 considered useful, and values greater than 0.9 considered excellent (Pennino et al. 2013). Similarly, TSS measures the ratio between specificity (i.e. the proportion of true negative correctly predicted that reflects a model’s ability to predict an absence given that a species in fact does not occur at a location) and sensitivity (i.e. the proportion of true positive correctly predicted that reflects a model’s ability to predict a presence given that a species in fact occurs at a location).

In order to quantify the predictive accuracy of the selected biomass and abundance model, predicted and observed values were compared using the Pearson correlation index.

The mean size of the exemplars was estimated by dividing biomass by number of individuals by country and area.

Finally, in order to detect significant differences among countries and areas in abundance, biomass and size, the non-parametric Kruskal-Wallis test and the Dunn-Bonferroni post-hoc method were applied to the mean values of the time period considered.

Maps were plotted using a depth threshold that respected the depth-stratification of the MEDITS surveys in order to avoid extrapolation in prediction.

The blackbellied anglerfish was caught in 4100 hauls (41.4%) of the 9906 performed during the MEDITS surveys (2006-2015) (Table 1).

The best binomial BRT RAC models retained all the variables, except for distance to coast. The percentage of total deviance explained was about 75% and bathymetry accounted for the vast majority of the relative influence of each variable (21.4%) (Table 2). Year explained 9.1%, slope 8.3%, SBT 8.2% and SBS 6.2%. The rest of the variability was accounted for by the autocovariate spatial term (21.8%) (Table 2). The blackbellied anglerfish showed its optimum occurrence between 150 and 350 m depth on flat sea bottoms, an SBS of 36.5 to 38 and an SBT of 16.5°C to 17.5°C.

Absence of spatial autocorrelation was found (Moran’s index=0.05; p-value<2.2e^–16), highlighting that the RAC method accurately handled this issue (Table 2). Model prediction performance statistics for the binomial BRT model achieved an AUC value of 0.82 (Table 2), which indicates an excellent degree of discrimination between locations in which the species is present and absent. Similarly, the TSS value was of 0.74 which indicates a good agreement between the predicted and observed values (Table 2).

The results showed several hot-spots of occurrence located in south and southeast Spain, the Balearic Islands, the Gulf of Lions, southern Italy (western Ionian Sea) and the Aegean Sea, followed by the waters around Elba island in the Tyrrhenian Sea and on the western Sardinian coast (Fig. 1 and Table 1). Although the general pattern is quite stable, some sites of these hot-spots change slightly among years.

For the biomass the best BRT RAC model selected as significative predictors SBT (17.2%), SBS (10.3%), bathymetry (8.5%) and autocovariate spatial term (12.5%) (Table 2). Slope, distance to coast and year were not retained in the final model. Similarly to occurrence, the bathymetry range of suitability habitats was 150 to 300 m depth with an SBT of 17.5°C to 18.5°C and a SBS of 37 to 38 (Table 2).

The Pearson correlation index showed a value of 0.88, highlighting a good ability of the model to predict biomass values. No spatial autocorrelation was detected (Moran’s index=0.04; p-value<2.2e^–16) for the final BRT biomass RAC model (Table 2).

Figure 2 shows the hot-spots identified with the highest biomass values of the blackbellied anglerfish. The most suitable areas for this species seem to be the northern Aegean and southeast Spain (GSA02 and GSA06) (Fig. 2). The biomass differed among countries (χ²=657.5; p<0.0001), showing no clear tendency from outer to inner parts of the Mediterranean Sea. The highest maximum biomass values were found in Greece (the Aegean Sea), Italy (the central and southern Tyrrhenian and western Ionian seas) and the lowest ones in France (Gulf of Lions) (Table 1). Differences among areas within each country were also detected (Fig. 2 and Table 1).

In Spain, Alboran Island showed higher biomass values than the other areas (p<0.0001). In France, the Gulf of Lions showed higher biomass than Corsica (χ²=11.474; p=0.01). The western Ionian Sea showed the highest biomass scores (χ²=154.09; p=0.0001) and the Tyrrhenian Sea the lowest. In Greece, the biomass showed a decreasing trend from the eastern Ionian to the Crete Sea (χ²=26.71; p=0.0001).

For abundance, the best BRT RAC model explains 65% of the variability with SBT (24.6%) as the main predictor, followed by SBS (15.6%), bathymetry (11.5 %), and autocovariate spatial term (12.3%) (Table 2). Slope distance to coast and year were not retained in the final model. Similarly to biomass, the bathymetry range of suitable habitats was 150 to 300 m, with an SBT of 17°C to 18°C and an SBS of 37-38 (Table 2).

The Pearson correlation index showed a value of 0.83, highlighting a good ability of the model to predict abundance values. No spatial autocorrelation was detected (Moran’s index=0.02; p-value<2.2e^–16) for the final BRT abundance RAC model (Table 2).

Figure 3 shows the identified hot-spots of abundance of the blackbellied anglerfish. The pattern is similar to that of biomass except on the Spanish coast, where abundance is higher.

Statistically significant differences in the abundance scores were found among countries (χ²=788.9; p=0.0001). In particular, Greece and Spain showed higher abundance values than France and Italy (Table 1). Abundance scores decreased from western areas towards the central part of the Mediterranean Sea and increased again sharply in Greece. Differences among areas in each country were also detected (Fig. 3 and Table 1).

In Spain, the Levant area and the Alboran Sea showed higher maximum abundance values than the Balearic Islands and Alboran Island (χ²=82.95; p=0.0001). In France, abundance was higher in the Gulf of Lions than in Corsica (χ²=105.6; p=0.0001). In Italy, the abundance scores decreased from Sardinia to the Ionian Sea (χ²=490.79; p=0.0001). In the Ligurian and northern and central Tyrrhenian seas, L. budegassa was more abundant in the central part of the area (i.e. the eastern Ligurian and northern Tyrrhenian seas) due to the wide extension of the continental shelf. Finally, in Greece the maximum abundance scores decreased slightly from the Ionian Sea to the Cretan Sea (χ²=6.04; p=0.049).

Individuals in Italy were the largest and in Greece they were the smallest. Spain and France showed intermediate values (χ²=813.44; p=0.0001) (Table 1).

Finally, the temporal trends of abundance and biomass of black anglerfish by geographical areas, shown in Figure 4, were consistent with those of the maps. In general, abundance and biomass increased from 2010 on.

The white anglerfish was caught in 2953 hauls (29.81%) of the 9906 performed during the MEDITS surveys (2006-2015) (Table 1). The main occurrence of the species during the whole time series was recorded in the Gulf of Lions (GSA07). Temporal variation of the occurrence was apparent in Spain, where a sharp increase in occurrence was found in 2010 and 2011, followed by a decrease. In Italy, similarly to the blackbellied anglerfish, the occurrence values were steadier during the study period.

The best binomial BRT RAC models retained all the variables except for distance to the coast and slope of the seabed. The total deviance explained was about 72.2%. In particular, SBT was the most important factor (22.2%), followed by year (11.3%), bathymetry (10.5%) and SBS (7.7%) (Table 2). The rest of the variability was accounted for by the autocovariate spatial term (20.6%). The higher probability of occurrence of this species was found between 200 and 600 m depth, with an SBS of 37.5-38 and a SBT of 15.5°C to 16.5°C (Table 2).

No spatial autocorrelation was detected (Moran’s index=0.02; p-value<0.10e^–16), and model prediction performance statistics achieved values higher than 0.70 (AUC=0.78; TSS=0.72), which indicate a good agreement between the predicted and observed values (Table 2).

Figure 5 showed the Gulf the Lions as the area with the highest probability of occurrence for this species, followed by the Aegean Sea and the western coast of Sardinia.

For the white anglerfish biomass model, the final predictors were SBT, bathymetry and the autocovariate spatial term. The total deviance explained was about 74%, with SBT as the main predictor (32%) for explaining the variability of the biomass of this species, followed by bathymetry (20.5%) and the autocovariate spatial term (20.5%). The biomass was higher at locations with a lower SBT, between 18°C and 19.5°C, and a depth of 200 to 400 m (Table 2).

The Pearson correlation index showed a value of 0.76, highlighting a good ability of the model to predict biomass values. No spatial autocorrelation was detected (Moran’s index=0.02; p-value<1.2e^–16) for the final BRT biomass RAC model (Table 2).

Figure 6 shows a differentiated geographical pattern of the white anglerfish biomass, highlighting that Sardinia, the Aegean Sea, the Spanish western Mediterranean coast and the Gulf of Lions are the areas with the highest biomass for this species (Fig. 6 and Table 1).

In France biomass values were higher in the Gulf of Lions than in Corsica (χ²=21.05; p=0.0001) (Table 1). In Italy the biomass importance decreased from Sardinia to the Ionian Sea and from the southern to northern Tyrrhenian Sea (χ²=543.65; p=0.0001). In Greece the biomass showed significant differences between areas (χ²=16.75; p=0.0001). The highest biomass was recorded in the northern Aegean. In Spain the highest biomass values were found in the Levant and Balearic Islands (northeast of Menorca Island). The rest of the areas had similar biomass scores.

Figure 7 shows a differentiated geographical pattern of white anglerfish abundance, highlighting Spain, Italy and Greece as the areas with the highest values (χ²=448.8; p=0.0001). For the white anglerfish abundance model, the final predictors were bathymetry, SBT and the autocovariate spatial term. The total deviance explained was about 67%, with bathymetry as the main predictor (25%) explaining the variability of the abundance of this species, followed by SBT (18.5%) and the autocovariate spatial term (23.5%). Abundance was higher in locations with a bathymetry range of 150 to 300 m and an SBT of 18°C to 19°C (Table 2).

The Pearson correlation index showed a value of 0.77, highlighting a good ability of the model to predict biomass values. No spatial autocorrelation was detected (Moran’s index=0.03; p-value<1.2e^–16) for the final BRT biomass RAC model (Table 2).

In Spain the highest abundance values were registered in the Balearic Islands and the Levant area (χ²=32.54; p=0.0001) (Table 1). In Italy the highest abundance scores were found around Sardinia and in the northern Tyrrhenian Sea (χ²=235.11; p=0.0001). In Greece the biomass showed no significant differences between areas, though in the northern Aegean Sea the scores were slightly higher. Finally, in France, although the occurrence scores were the highest of all the studied Mediterranean areas, the abundance values were very low.

Italy showed the largest sizes and Spain and France the smallest; Greece showed intermediate values (χ²=218.108; p=0.0001) (Table 1).

Finally, the temporal trends of abundance and biomass of white anglerfish by geographical areas shown in Figure 4 were consistent with those of the maps. In general, abundance and biomass showed a slight decrease towards 2015.