Spatio-temporal variability in the distribution pattern of anglerfish species in the Mediterranean Sea Mediterranean demersal resources and ecosystems: 25 years of MEDITS trawl surveys

Summary: The Mediterranean distributions of two species of anglerfish, the blackbellied anglerfish ( Lophius budegassa ) and the white anglerfish ( Lophius piscatorius ), were analysed from trawl survey data (MEDITS project – Spain, France, Italy and Greece) from 2006 to 2015 implementing a Delta model approach with residuals autocovariate boosted regression trees. Sea bottom temperature (SBT), sea bottom salinity (SBS), bathymetry, slope of the seabed and distance to the coast were considered possible predictors. The results show that the locations with a higher presence, abundance and biomass of L. budegassa are those with a depth range between 150 to 300 m, with an SBT range between 17.5 and 18.5 ° C, and SBS of 37-38 PSU. Similarly, L. piscatorius shows a higher probability of presence, abundance and biomass in location with a bathym- etry range of 200-400 m, an SBT of 17.5 ° C to 18.5 ° C and an SBS of 36.5 to 37.5. Our results identify preference habitats for the anglerfishes in the Mediterranean Sea such as the Aegean Sea, the Gulf of Lions, south and southeast Spain and the northwestern Ionian Sea. In general terms, these findings enhance our understanding of the differences in the spatio-temporal distribution of these two species, providing useful information that can help their fisheries management and conservation.


INTRODUCTION
The anglerfish Lophius budegassa, Spinola, 1807 and Lophius piscatorius Linaeus, 1758 are widely distributed species in the Mediterranean basin and the northeastern Atlantic Ocean (from the British Isles to Senegal) in different bathymetric strata from 50 m up to around 1000 m depth (FAO 2007). They are usually caught by multi-species trawlers, and the commercial importance of their catches varies from the Atlantic Ocean to the Mediterranean Sea, with high variability between areas.
Although both anglerfish are commercial species in both seas, they are more important target resources in the Atlantic than in the Mediterranean Sea, where they are often caught as a by-catch.
However, the European Commission's Scientific, Technical and Economic Committee for Fisheries (STECF) already conducts standard stock assessments of both species in several areas of the Mediterranean Sea-for example, for the L. budegassa stocks located in the western and central parts of the Mediterranean basin, namely in geographical subareas (GSAs) 5, 6, 7, 15-16 and 22 (Vasilakopoulos and Maravelias 2016 for 22;STECF 2012STECF , 2013GFCM-SAC 2012, GFCM-SAC-SCSA 2014.
In the Atlantic, studies of biology, age and growth have been conducted mainly for L. piscatorius (Duarte et al. 1997, Quincoces et al. 1998, 2008, Velasco et al 2008. In addition, its historical catch trends and statistics were evaluated and the exploitation status of its stocks were established (ICES 2007, GFCM-SAC-SCSA 2014, STECF 2012, 2013. In the Mediterranean, Ungaro et al. (2002), Maravelias and Papaconstantinou (2003), García-Rodriguez et al. (2005) and Massaro et al. (2016), studied different aspects of their biology and life history, and Maravelias and Papaconstantinou (2003) and Carlucci et al. (2009) studied their geographical distribution. However, most of these studies were local and used short time series, Ungaro et al. (2002) being the only one with a broad Mediterranean spatial and temporal scale. Moreover, it is important to note that, although the FAO catch statistics in the Mediterranean Sea have remained stable over the last 25 years, the last stock assessments carried out noted that overall the blackbellied anglerfish has an overfished status in the Mediterranean.
Within this context, the aim of this study was to explore the occurrence, biomass and abundance patterns of the two anglerfish species in the Mediterranean Sea on a broad spatio-temporal scale. Indeed, despite the commercial importance of both anglerfish species, many aspects of their distribution and abundance patterns still remain unknown for the Mediterranean. In this paper the Mediterranean International Trawl Survey (MEDITS) datasets from 2006 to 2015 of Spain, France, Italy and Greece were analysed in combination with environmental data using a Delta model approach with residuals autocovariate boosted regression trees (RAC-BRT).
The results of this study will provide new information on the spatial and temporal distribution of these species that will be valuable to fishery managers for a marine spatial planning approach. The implementation of an ecosystem approach to fisheries management and marine spatial planning will include the protection of priority habitats for these species and proper sustainability of their fisheries. The essential prerequisites for achieving this objective include a solid knowledge of species-environment relationships and identification of priority areas using robust analysis of existing information and databases (Pennino et al. 2013).

Species data
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 37-38 PSU. De igual manera el rape blanco (Lophius piscatorius) muestra la probabilidad más alta de presencia, abundancia y biomasa en localidades con profundidades comprendidas entre 200-400 m, una SBT entre 17.5-18.5°C y 36.5-37.5 PSU de SBS. Nuetros resultados identifican como hábitats preferenciales de los rapes el mar Egeo, el Golfo de Lion, el sur y sureste de España y el noroeste del mar Jónico. En términos generales estos hallazgos proporcionan una información que mejora nuestra comprensión sobre las diferencias espacio-temporales de la distribución de las dos especies de rape, proporcionando una información muy útil que puede ayudar en la gestión de su pesquería y conservación.
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.

Environmental data
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)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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.

Model approach
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 calibrationvalidation 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 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 (χ 2 =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 (χ 2 =11.474; p=0.01). The western Ionian Sea showed the highest biomass scores (χ 2 =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 (χ 2 =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 (χ 2 =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 (χ 2 =82.95; p=0.0001). In France, abundance was higher in the  Gulf of Lions than in Corsica (χ 2 =105.6; p=0.0001). In Italy, the abundance scores decreased from Sardinia to the Ionian Sea (χ 2 =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 (χ 2 =6.04; p=0.049).
Individuals in Italy were the largest and in Greece they were the smallest. Spain and France showed intermediate values (χ 2 =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.

White anglerfish (Lophius piscatorius)
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 (χ 2 =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 (χ 2 =543.65; p=0.0001). In Greece the biomass showed significant differences between areas (χ 2 =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 (χ 2 =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 (χ 2 =32.54; p=0.0001) ( Table 1). In Italy the highest abundance scores were found around Sardinia and in the northern Tyrrhenian Sea (χ 2 =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 (χ 2 =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.

DISCUSSION
This study investigates the possible relationships between the occurrence, biomass and abundance of two anglerfish species and environmental variables in the Mediterranean Sea. The use of quantitative models including environmental variables is complementary to previous studies that roughly identified the spatiotemporal variations of these species (Ungaro et al. 2002). In addition, the incorporation of a spatial correlation contiguity matrix as an explanatory variable in the BRT models, using the residuals autocovariate approach, accounted for the autocorrelated nature of the data. This methodology recently developed by Crase et al. (2012) and successively implemented by Escalle et al. (2016) showed good predictive results for both studied species, highlighting its efficiency to deal simultaneously with the spatial correlation, non-linearity and the large proportions of zeros.
The main predictor of population distribution for L. budegassa was depth, while the main predictors for L. piscatorius were temperature and salinity. It should be noted that, for both species, the spatial component was one of the most important variables selected by the models. This component indicates the intrinsic spatial variability of the data after the exclusion of the environmental variables, highlighting that other important spatial processes not taken into account in the model are affecting the distribution, abundance and biomass of these species. Indeed, in demersal species such as anglerfishes, the role of other abiotic factors (such as the type of bottom) could affect the spatial distribution  (Smith et al. 2008). This feature is also directly related to other factors, such as the type of prey, and biotic processes, such as competition, predation and recruitment, which are also spatially structured (Katsanevakis and Maravelias 2009, Stagioni et al. 2013. The occurrence distribution of both species evidenced some differences. In particular, L. budegassa showed a higher degree of occurrence than L. piscatorius, in accordance with past findings (Ungaro et al. 2002, Maravelias and Papaconstantinou 2003, Carlucci et al. 2009). L. budegassa showed a very high occurrence in the eastern and western Mediterranean at depths around the upper slope between 150 and 350 m (continental shelf), while L. piscatorius was mainly present in the Gulf of Lions and in the northern Aegean Sea, showing a deeper distribution between 200 and 600 m (upper part of the continental slope and slope). These differences could be due to their different environmental preferences. Indeed, preferential habitats of L. budegassa were those with a depth range of 150 to 300 m and an SBT range of 16.5°C to 17.5°C and an SBS of 36.5-37.5, while preferential habitats of L. piscatorius were the those with a depth range of 200 to 600 m, an SBT of 15.5°C to 16.5°C and an SBS of 37.5 to 38.
L. budegassa showed persistency in its occurrence areas, suggesting that these locations are highly suitable, as observed for other species (Paradinas et al. 2015). Conversely, habitat suitability for L. piscatorius was more heterogeneous.
However, it is worth noting that preferential habitats of both species were observed in areas close to marine canyons (the Alboran Sea and the Levant Sea in Spain; the Gulf of Lions and west of Corsica in France; Sardinia and the western Ionian Sea in Italy; and the southwestern Aegean in Greece) and delta rivers (e.g. the Ebro Delta in Spain, and the deltas of Strymon, Axios and Evros in the northern Aegean), or in hydrographically dynamic areas (Alboran Island and the Aegean Sea). All these locations are characterized by strong flows, ridges and high run-off areas that boost primary production availability and hence create a preferential habitat for many species (Maravelias and Papaconstantinou 2003). In addition, the canyons can act as a kind of ecological refugee or natural reservoirs. Taking into accounts all these considerations, larger and older specimens of both species could survive being excluded from fishing exploitation, helping the renewal of stock in neighbouring grounds where trawling occurs more often (Carlucci et al. 2009). These results are also in line with previous findings of other studies in the same area for other pelagic and demersal species (Bellido et al. 2008, Paradinas et al. 2015. L. budegassa showed higher abundance than L. piscatorius, whereas the opposite was observed for biomass, suggesting the presence of more small-sized individuals of L. budegassa and fewer individuals of L. piscatorius but of larger size. These results could be attributed to lower exploitation rates of L. piscatorius because of its deeper bathymetric preference. This trend is evident in waters around Sardinia, perhaps suggesting that the fishing pressure is still tolerable for the populations in these areas. The low exploitation level of both anglerfish species could be due to the particular composition of the Sardinian fishing fleet. Indeed, in this area the fleet mostly comprises boats belonging to artisanal fishery that typically employ highly selective gears, in terms of both fish species and sizes, while the trawl sector represents only 11% of the overall Sardinian fishing vessels (Anonymous 2011). Another possible explication could be the adoption by the Autonomous Region of Sardinia of specific minimum landing sizes for both species (20 cm total length for L. budegassa and 40 cm total length for L. piscatorius) (D.A.D.A.R.S. NR. 412 DEL 10.05.1995; http://www.sardegnaambiente. it/documenti/19_425_20140214113336.pdf). Our results suggest that these types of species-specific conservation measures are being effective in this area and should probably be extended to other Mediterranean regions.
In the central Tyrrhenian Sea, although the continental shelf is quite extensive, the abundance of L. budegassa is limited due to the high fishing pressure (Massaro et al. 2016). These results are consistent with the size of the exemplars in this area. On the other hand, in the northern Tyrrhenian the individuals are bigger, as in Corsica waters and the Ionian Sea.
Our findings on the abundance, biomass and estimated size trends also suggest that Greece, the Gulf of Lions and the southern Tyrrhenian are possible recruitment areas for L. budegassa, while the Levant of Spain and the Gulf of Lions are possible recruitment areas only for L. piscatorius. The effects of reduction of fishing pressure on the fish population structure have been widely described. One of these effects occurring in several areas is the increase in the mean fish size (Botsford et al. 1997, Bianchi et al. 2010, Smith et al. 2008. Our results are consistent with to this thesis and also with results previously found in the same area (García-Rodriguez et al. 2005). In Alboran Island, because of its distant location and the narrow strip of available fishing grounds, fishing pressure is considerably lower than in other areas. Those conditions could explain the presence of larger individuals of both anglerfish species than in other Spanish waters. On the other hand, waters surrounding the Columbretes Islands in Spain have been already highlighted as a possible nursery area for many demersal species. These islands have been a marine protected area since 1989, and according to Paradinas et al. (2015) they provide a stable, high-quality ecosystem for stocks that could be exporting adults and recruits to adjacent areas. In addition, López et al. (2016) showed that juvenile specimens of both species were mainly located in coastal and continental shelf areas in northern Spain.
Our results of biomass, abundance and mean estimated size show that Italian waters (the northern Tyrrhenian and Ionian seas) concentrate the largest individuals of both species, probably because of the significant feeding grounds in the area, while L. pisca-torius is also found mainly in French waters and in the northern Aegean Sea.
In agreement with our results, Maravelias and Papaconstantinou (2003) also highlighted that the Aegean Sea is an apparent recruitment area for L. piscatorius. Oceanographic studies , Zervakis et al. 2000 indicated that intensive vertical mixing takes place in the Aegean Sea, promoting the post-larval transportation of this species from local recruitment areas to deeper waters.
MEDITS survey data used in this study covered only a short season from late spring to early summer. Consequently, the fitted models can only reflect a snapshot view of the expected relationships and distribution of these species.
Our results could contribute to the characterization of marine protected areas and any other measures of spatial fishing management. Indeed, effective marine spatial planning should include seasonal spatio-temporal dynamics of these species, of course also considering the particular signals from fishery-dependent data to complete this view. These results, ideally also expanded to incorporate multiple species, segregation between species and life stages, could configure a solid approach to implementing an effective marine spatial planning that embraces the requirements for an ecosystem approach to fisheries management.