Spatial distribution pattern of European hake, Merluccius merluccius (Pisces: Merlucciidae), in the Mediterranean Sea Mediterranean demersal resources and ecosystems: 25 years of MEDITS trawl surveys

repro-Summary: The present study provides updated information on the occurrence, abundance and biomass distribution patterns and length frequencies of Merluccius merluccius in the Mediterranean Sea, by analysing a time series of data from the Mediterranean International Trawl Surveys (MEDITS) from 1994 to 2015. The highest values of abundance and biomass were observed in the Sardinian Seas. The use of a generalized additive model, in which standardized biomass indices (kg km –2 ) were analysed as a function of environmental variables, explained how ecological factors could affect the spatio-temporal distribution of European hake biomass in the basin. High biomass levels predicted by the model were observed especially at 200 m depth and between 14°C and 18°C, highlighting the preference of the species for colder waters. A strong reduction of biomass was observed since the year 2009, probably due to the strengthening of the seasonal thermocline that had greatly reduced the availability of food. The general decrease in biomass of several stocks of anchovy and sardine, preys of European hake, might be indirectly connected to the decreasing biomass detected in the present study. The length analysis shows median values lower than 200 mm total length of most of the investigated areas.


INTRODUCTION
Multispecies fisheries targeting a large number of species represent a management challenge because several biological and ecological parameters must be incorporated in the assessment (Ekerhovd and Steinshamn 2016). This is the case of the Mediterranean fisheries, where the commercial catches are characterized by a high number of target species, among which the most ecologically and economically important one is the European hake, Merluccius merluccius (Linnaeus, 1758) (Orsi Relini et al. 2002).
M. merluccius is widely distributed in the Mediterranean Sea and the northeastern Atlantic Ocean, from the coasts of Norway and Iceland to Mauritania. In the Black Sea, it only occurs along the southern coasts. It is a major component of the demersal fish assemblages and can reach depths ranging from 30 to 1000 m, although the most abundant catches are observed especially between 50 and 400 m depth (Orsi Relini et al. 2002).
The European hake has always been an important food resource for the population of western Europe throughout history. It is caught in mixed fisheries along with blue whiting, megrim, monkfish, shortfin squid, Norway lobster and deep-water rose shrimp by a multirigged fleet. This demersal species is one of the most important commercial resources for trawl and small-scale fisheries (fixed gillnets, bottom longlines) in the Mediterranean basin (STECF 2015). Trawlers mainly catch specimens below 30 cm total length, whereas bottom longlines and gillnets target adults (Sbrana et al. 2007).
European hake represents about 2% of the average total landings in the Mediterranean Sea (FAO 2016). However, according to the data series, annual catches decreased from a maximum of ca. 52000 tonnes in 1994-1995(Maynou et al. 2003 to 24900 tonnes in the period 2000-2013, indicating an overfishing status in several Mediterranean areas (STECF 2015, GFCM 2016. Several scientists have highlighted the risk of a stock collapse (Russo et al. 2017); indeed, European hake stocks in the Mediterranean Sea suffer the highest fishing pressure, with a fishing mortality rate that is on average five times higher than the target fishing mortality level (FAO 2016).
Since 1994, the MEDITS project (Bertrand et al. 2002) has provided the opportunity to collect biological and abundance data on demersal species at a Mediterranean Basin scale, by means of large-scale experimental trawl surveys and standardized sampling protocols. In this context, Orsi Relini et al. (2002Relini et al. ( ) presented data (1994Relini et al. ( -1999 on abundance, population structure and recruitment pattern of European hake from the Alboran Sea to the Aegean Sea, providing the first integrative picture of European hake over the northern Mediterranean Sea.
The present study aims to provide an update of previous studies, presenting more advanced analyses on several aspects of this species, in terms of both methodological approach and time period of investigation, also analysing several drivers shaping its distribution. Considering MEDITS time series from 1994 to 2015, updated data on occurrence, abundance patterns and mean length of M. merluccius in the northern Mediterranean Sea are presented. Finally, the spatial distribution of biomass of this demersal species as a function of environmental variables was analysed and a predictive map of biomass hotspots was obtained.

Sampling
Since 1994, twenty-two annual bottom trawl surveys have been performed covering the European Union Mediterranean Sea in the framework of the MEDITS project. These surveys mainly aim to obtain estimates of abundance indices and length frequency distributions for a set of target species and have usually been carried out from late spring to mid-summer on demersal grounds from 10 to 800 m depth. A stratified random sampling design was adopted considering five bathymetric strata: 10-50 m, 51-100 m, 101-200 m, 201-500 m and 501-800 m. Details of the survey design and sampling methodology can be found in Bertrand et al. (2002) and in the MEDITS Handbook (Anonymous 2017).
In the present study, catch data (number and weight) and total length measurements of M. merluccius collected in the period 1994-2015 from 16 geographical sub-areas (GSAs) (FAO 2009) belonging to nine countries (Table 1, Fig. 1) were analysed. In GSAs 5, 17, 20, 22 and 23, the time series are not continuous because MEDITS surveys were not performed for some years; however, they were included in the analysis in order to cover the largest geographical patterns possible.
A total of 23748 valid hauls were carried out during the MEDITS surveys ; for each haul, number of specimens and total weight (g) per species were collected. For a list of target species, including European hake, biological sampling consisted in recording size (total length, TL in mm), gonad weight (g), sex and gonad maturation stage.

Data analysis
The frequency of occurrence of European hake (F(%)=N positive hauls/total hauls×100) was computed for each GSA considering the overall depth range (10-800 m). The number of individuals (N) and weight (kg) were standardized to a trawled surface unit (km 2 ) as abundance (N km -2 ) and biomass (kg km -2 ) indices. Net horizontal opening and recorded distances were used to calculate the swept area. Their distributions were represented graphically for each GSA by means of box plots showing medians, quartiles and ranges (maximum and minimum).
One-way analysis of variance (ANOVA) (Yandell 1997) was performed on log-transformed abundance (N km -2 ) and biomass (kg km -2 ), followed by pairwise post hoc Tukey honestly significant difference tests to investigate for significant differences in terms of abundance and biomass values of M. merluccius among GSAs. Normality assumption of abundance and biomass distributions were previously evaluated using the Shapiro-Wilk test.
The biomass indices (BI, kg km -2 ) were modelled using generalized additive models (GAMs), which provides a flexible framework for modelling the relationship between the individual predictors and the dependent variable through smooth functions that can be both linear and nonlinear (Wood 2017).
The explanatory variables considered in the model were the geographic coordinates (latitude and longitude) combined in a bidimensional smoother, the mean haul depth (m), sea surface temperature (°C, SST), temperature at the sea bottom (°C, botT) and the surface chlorophyll concentration (mg m -3 , CHL). These variables were selected for the GAM models because they showed considerable variability throughout the study area and were believed to be easier to interpret biologically than some other variables that were available.
The year of the surveys was also included as a factor. Depth grid was obtained from the EMODnet Bathymetry portal (http://www.emodnet.eu/), while data of the other covariates (SST, botT and CHL) were extracted from Mediterranean Monitoring and Forecasting Centre products available on the Copernicus Marine Environment Monitoring Service (CMEMS) web site (http://marine.copernicus.eu). In particular, SST, and botT data were produced using the physical reanalysis component of the Mediterranean Forecasting System (Simoncelli et al. 2014), while CHL data were produced using the MedBFM biogeochemical model (Teruzzi et al. 2014).
Absence of collinearity among covariates was checked by means of variance inflation factor (VIF) values (Zuur et al. 2010). Only variables with a low Pearson correlation coefficient (r) (|r|<0.5) and a low VIF value (VIF<3) were considered not correlated.
BIs were modelled using Gaussian distribution with a logarithmic link function, according to the following specification: log (biomass) = α + s 1 (lon, lat) + s 2 (depth) + s 3 (SST) + + s 4 (botT) + s 5 (SAL) + s 6 (CHL) + factor (year) where α represents the intercept, s i (i=1..., 6) are the smooth functions of the covariates and factor is the step function for the categorical variable year. It is worth noting that, in the model, the smooth function of spatial coordinates represents a residual spatial variation, which is spatially structured. The selection of the best model was realized by minimizing the Akaike information criterion value with the introduction of each variable in the analysis. Residuals of the model were tested with the Anderson-Darling normality test (Gross and Ligges 2015), and the symmetry was assessed by estimating the skewness index of the distribution (Hair et al. 2010). The Spearman rho correlation coefficient was carried out to test the presence of trends in residuals.
For predictions, a 0.03×0.03° resolution grid was used. The bathymetrical extension of the grid was within 10 to 800 m depth. The prediction grid was used to sample the informative layers of the environmental variables (monthly data) used in the model. The grid's points were also coupled with year and coordinates variables in order to estimate the annual predictive maps throughout the GAM model. In particular, the most representative predictive maps, for the years 1995, 2006, 2011, 2013 and 2015, were chosen to represent biomass values (kg km -2 ) of European hake, and to check the model robustness in space, the annual predictions of the model were coupled with the estimation of the mean stratified BIs per GSA obtained following the mean stratified standardization formulas by Souplet (1996).
In order to identify the hotspots of the European hake biomass distribution pattern, prediction maps for the last three years (2013-2015) of the time series were averaged in a single map. This raster map was used to extract the 80 th , 90 th and 95 th percentiles of distribution to draw the actual hotspot contours.
Statistical analyses were performed using R version 3.5.2 (R Development Core Team 2018). In particular, GAMs were implemented by means of the mgcv R package (Wood 2017). Biomass maps were processed and produced with QGIS software (QGIS Development Team 2017).
For each GSA, the length frequencies of M. merluccius were graphically represented by means of box plots of median, quartile, minimum and maximum TLs. An ANOVA test was applied to investigate significant differences among GSAs in terms of size distribution.
The GAM model of biomass explained 82.2% of the variability in the data (deviance explained). All the smoothing terms were significant (p<0.05). The degrees of freedom estimated for all smoothers are reported in Figure 3.
Although residuals of the model tested with the Anderson-Darling normality test (p<0.05) were not normally distributed, the model was considered sufficiently symmetric because its skewness index (0.023) was within the ±1 range. Moreover, the Spearman rho test carried out on the model's residuals found no significant trend (p>0.05).
Plots of the smoothers and diagnostic plots of the residuals resulting from the model are presented in Figures 3 and 4. The model indicated that high BI levels of European hake were mainly detected in the northern Spain area, in the Gulf of Lions, on the western side of Sardinia, in the northern Adriatic Sea and in the northern area of the Aegean Sea. The highest biomass of European hake was identified by the model at 200 m depth, although the optimal depth range can extend up to 400 m depth. In relation to sea surface temperature, fairly stable biomass values were observed between 20°C and 22°C. Considering the bottom temperature's spline shape, the partial effect of the variable shows an increase in the response variable corresponding to temperatures of between 14°C and 18°C, while a progressive reduction of the biomass was expected at higher bottom temperatures. Also, chlorophyll concentration was linked with European hake presence: the partial effect of the CHL variable, as is shown by the spline shape, describes an increase in biomass at chlorophyll concentrations higher than 0.6 mg m -3 , and-albeit for few values-another peak was also evident around higher values than 1.7 mg m -3 . The biomass index showed a fairly stable level over time up to 2002 and large fluctuations thereafter. In particular, an increase of the index recording the highest predicted values  GSAs. After this period, a reduction of the BI index was observed, with a relative peak in 2013 and the lowest values in 2011. In relation to the year's effect (Fig. 3), the most representative predictive maps for the years 1995 (start, medium level), 2006 (high level), 2011 (low level), 2013 (high and recent) and 2015 (low and recent) are shown in Figure 5. For the same years, the observed biomass values from MEDITS data were compared with the predictions of the model; comparable tendencies of biomass values were observed for the different GSAs (Fig. 5).
The areas with the greatest abundance of European hake were located in the Gulf of Lions, the eastern Adriatic Sea, the western side of Sardinia, the south of Sicily and the northern Aegean Sea. This distribution pattern was confirmed and better evidenced by the map showing the biomass hotspots of the European hake for the years 2013-2015 (Fig. 6) calculated at the 80 th , 90 th and 95 th percentiles of the BI distribution.
The length frequencies of M. merluccius are summarized in the box plots shown in Figure 7. In the western Mediterranean, the largest median value was found in Corsica (230 mm TL), while the smallest median values were detected in northern Spain (125 mm TL) and the Gulf of Lions (130 mm TL). Median values lower than 200 mm TL were also observed in the northern Alboran Sea and in the Balearic Islands. In the central Mediterranean, the median values were below 200 mm TL in all areas (the Ligurian and Tyrrhenian seas, Sardinian waters, south of Sicily, the Adriatic Sea and the northwestern Ionian Sea). A prevalence of median values larger than 200 mm TL was observed in the eastern Ionian Sea except the Cretan Sea (160 mm TL). ANOVA results also suggested a significant difference among GSAs in terms of median size (F=1684, p<0.001).

DISCUSSION
The present study provides updated information on the occurrence, abundance and biomass distribution patterns and length frequencies of M. merluccius in the Mediterranean Sea on a large spatial scale during spring-summer over the last two decades. In particular, the use of a GAM model, in which standardized indices of biomass were analysed in relation to the effect of environmental variables, explained how ecological factors could affect the spatio-temporal distribution of European hake biomass in the area.
The results confirm the wide spatial distribution of the species, in terms of both bathymetric range and occurrence, in accordance with previous studies in several Mediterranean areas (Orsi Relini et al. 2002, Katsanevakis et al. 2009). The significant differences among GSAs in terms of both abundance and biomass highlight a high variability in the catches of this demersal species in the whole area, in agreement with the findings by Orsi Relini et al. (2002).
All environmental factors were significant drivers affecting European hake distribution. The application of the GAM model, which considers the effect of some environmental variables, has significantly improved knowledge of the spatio-temporal variations of M. merluccius in the Mediterranean Sea. The high biomass levels predicted by the model in northern Spain (the Catalan Sea), the Gulf of Lions, the western side of Sardinia, south of Sicily, the eastern Adriatic Sea (mainly Croatian waters) and the northern Aegean Sea are in accordance with previous studies carried out in the Mediterranean Sea (Orsi Relini et al. 2002, Hidalgo et al. 2008, Druon et al. 2015. Furthermore, the predictive model showed a biomass peak at 200 m depth, most likely correlated to the presence of nursery grounds in many Mediterranean areas (Orsi Relini et al. 2002, Druon et al. 2015. The optimal depth range for BI was observed between 200 and 400 m, beyond which a consistent decrease in biomass was evident throughout the time period. These results are in agreement with previous studies using the MEDITS protocol, which indicate that abundance of European hake increases especially at depths ranging from 100 to 400 m, but decreases at depths greater than 500 m (Orsi Relini et al. 2002, Katsanevakis et al. 2009, Garofalo et al. 2018). This demersal species shows an optimal bottom temperature range of between 14°C and 18°C, with noticeably lower biomass for higher bottom temperatures. The GAM model confirms the results obtained by De Pontual et al. (2013) concerning the preference of the species for cold waters, as it has a limited tolerance to warm temperature.
The CHL covariate also has a significant effect: an increase in biomass occurs with increasing CHL concentration, which might be considered a proxy of primary production, thus indicating areas where favourable feeding conditions would attract the species (Brosset et al. 2017). M. merluccius is an opportunistic feeder (Hidalgo et al. 2008) that is able to feed on a wide range of trophic resources. In the central Mediterranean, Carpentieri et al. (2005) reported that European hake size classes are differentiated along food niche dimensions according to different prey sizes or different prey taxa. During its early demersal life, the species feeds on crustaceans (mostly euphausiids and mysids); subsequently the juvenile hakes migrate from the nursery areas to the parental stock, and when they reach between 18 and 32 cm TL, they shift their diet towards pelagic and necto-benthic fish (Carrozzi et al. 2019). The predicted biomass hotspots of the European hake could be linked to areas characterized by high availability of crustaceans and small pelagic fish consumed mostly by individuals of age 0 and large specimens of European hake respectively.
This study highlights that the effect of environmental factors such as the range of preferred bottom temperature (14°C-18°C) and the increase in chlorophyll concentration could offer favourable conditions for a productive habitat, enhancing the prey availability and consequently increasing the biomass levels of European hake in Mediterranean areas. The results confirm the potential distribution of juvenile European hake habitats reported by Druon et al. 2015, with the highest biomass levels of M. merluccius occurring in nursery areas located near the shelf break in the northern Catalan Sea, the eastern Gulf of Lions, southwestern Sardinia, the northern Tyrrhenian Sea, southeastern Sicily, the semi-enclosed part of the eastern Ionian Sea and fragmented areas of the Aegean Sea. The prediction of both important nursery habitats (Druon et al. 2015) and biomass hotspots in the northern Aegean Sea is not revealed by MEDITS survey data (Tserpes et al. 2008). The low presence of juveniles in the northern Aegean is rather due to a time lag between the MEDITS sampling period and the seasonal recruitment peak in autumn-winter (Hellenic Centre for Marine Research, unpublished data). The high biomass observed in the present study may confirm the presence of a nursery area in the southern Adriatic Sea (Pomo Pit), as hypothesized by Arneri and Morales-Nin (2000). The high biomass of small pelagic fish in the Mediterranean basin (Palomera et al. 2007, Giannoulaki et al. 2011) offers a great lure of food for the species. In the southern-central Mediterranean Sea, Garofalo et al. (2018) reported that high small pelagic fish biomass related to the high primary production had led to the presence of large adult European hakes in the Gulf of Gabes. The strong relationship between productive fronts and demersal resources especially linked to fish preys of low trophic level has been observed in other marine areas (Alemany et al. 2014).
The predictive maps showing the European hake biomass index along a temporal scale highlighted a strong reduction of biomass since the year 2009 according to FAO (2016), which reports a decrease in European hake landings in Spain, Italy and France since 2009. Also, in the ecological niche modelling developed by Druon et al. (2015) on a Mediterranean scale, a substantial decrease in biomass of recruits and adults in 2010 and 2011 was detected, probably due to the strengthening of the seasonal thermocline that had greatly reduced the availability of food. Also, a recent report by the GFCM (2019) highlights that the increase in sea bottom temperature could likely be responsible for the decrease in hake nurseries in recent years.
Environmental factors such as temperature and food availability can also affect the spatial distribution of the biomass of European hake through indirect effects. The general decrease in biomass of several stocks of anchovy and sardine observed at a Mediterranean scale (Vasilakopoulos et al. 2014) might be connected to the decreasing M. merluccius biomass detected in the present study along the time series. The effects of climate change (e.g. increase in mean SST, sea surface salinity and water stratification) are favouring smaller size classes in the zooplankton community (Herrmann et al. 2013). Since smaller plankton has lower energy values (Zarubin et al. 2014, Brosset et al. 2017, the condition of small pelagic fish may remain poor or even continue to decline, in turn having a negative effect on the European hake biomass, which would tend to decrease over time due to the reduced quality and/or quantity of food. Other investigations document that the decrease in biomass of small pelagic fish could arise from the combined effect of exploitation and environmental changes, especially in areas where anchovy and sardine were clearly overexploited, such as the Adriatic Sea (FAO 2014). Indeed, such decrease in the biomass of anchovy and sardine stocks could have negative effects on many Mediterranean fisheries, hindering productivity of both small pelagic stocks and species preying on small pelagic fish, such as European hake.
The length analysis confirms, as other studies have reported (Orsi Relini et al. 2002, Ligas et al. 2015, that trawl fisheries in Mediterranean waters mainly target first age classes (0-1) corresponding to sizes of between 100 and 180 mm TL. The smaller median sizes of M. merluccius observed in the western Mediterranean (northern Spain and the Gulf of Lions) and in the central Mediterranean, mostly in the Ligurian and Cretan seas, probably confirm concentrations of juveniles in these areas (Hidalgo et al. 2009), whereas the largest median values, appearing in Corsican waters, the Aegean Sea and Cyprus waters, reflect the lack/ low presence of recruits in these areas (Orsi Relini et al. 2002, Druon et al. 2015 due to the recruitment seasonal peak (autumn-winter), which escapes from MEDITS sampling.
The present study has shown good predictive results highlighting how environmental factors can affect the spatio-temporal distribution pattern of the European hake in the Mediterranean basin. The results emphasize the importance of integrating any existing knowledge about environmental drivers in the assessments of the ecosystems and fisheries, in line with the ecosystem approach to fisheries management and the requirements of the European Commission (EC 2008) to the Scientific, Technical and Economic Committee for Fisheries (STECF) and the International Council for the Exploration of the Sea (ICES). Furthering knowledge of the spatio-temporal distribution pattern of European hake and identifying persistent areas of high abundance of juveniles and spawners in relation to environmental factors are key elements for the management of this important commercial species that has long been overfished and is now subject to multi-annual management plans in several GSAs of the Mediterranean Sea (FAO 2016).