Distribution and spatio-temporal biomass trends of red mullets across the Mediterranean Mediterranean demersal resources and ecosystems: 25 years of MEDITS trawl surveys

Summary: The present work examines the spatio-temporal biomass trends of Mullus barbatus and Mullus surmuletus in the Mediterranean Sea through the analysis of a time series of data coming from the Mediterranean International Trawl Surveys (MEDITS), accomplished annually from 1994 to 2015. The biomass of both species showed clear declining trends below 150 to 200 m depth, which were steeper in the case of M. barbatus. Increases in temporal biomass trends were observed for M. barbatus from 2008 onward in most geographic sub-areas (GSAs), while stability was mostly observed for M. surmule- tus. For both species , dynamic factor analysis revealed similarities among neighbouring GSAs and the subsequent cluster analysis identified two major GSA groups corresponding to the eastern and western basins of the Mediterranean. Overall, the results suggested that the combined effects of fishing and environmental conditions determine species abundance variations, but the relative importance of each component may vary among areas.


INTRODUCTION
Red mullet (Mullus barbatus L., 1758) and striped red mullet (Mullus surmuletus L., 1758) are demersal fish distributed throughout the Mediterranean and Black seas, as well as in the eastern Atlantic, from the North Sea to Senegal (Fischer et al. 1987). They inhabit depths down to 300 to 400 m, with red mullet preferring muddy bottoms, while striped red mullet is generally found in bottoms with heterogeneous granulometry and even in Posidonia beds.
Both are species of high commercial value and are the main targets of several Mediterranean demersal fisheries (Farrugio et al. 1993, Relini et al. 1999, Papaconstantinou and Farrugio 2000. Red mullet in particular is one of the most important targets of several bottom-trawl fisheries operating along the continental shelf. It is generally considered that its stocks undergo high fishing pressure, but recent studies suggest that they are in better condition than other commercially exploited demersal species (Vasilakopoulos et al. 2014, Tserpes et al. 2016, Cardinale and Scarcella 2017. Various studies have examined the distribution of M. barbatus and M. surmuletus in Mediterranean areas and have commented on their bathymetrical distribution and abundance variations (e.g. Machias et al. 1998, Tserpes et al. 1999, Lombarte et al. 2000. However, the employment of different sampling schemes has not allowed for direct comparisons among areas that could provide a much broader picture of the species distribution. A previous analysis of data coming from the in-ternational bottom-trawl survey (MEDITS) identified Mediterranean-wide abundance trends for both species for a time period of 7 years . The analysis, however, did not focus on the regional level, so no information regarding similarities/dissimilarities among areas was provided.
In the present work, we analyse a longer series of data from the MEDITS surveys with the aim of identifying biomass trends at the geographical sub-area (GSA; Fig. 1) level established by the General Fisheries Commission for the Mediterranean (GFCM), as well as to detect general patterns prevailing at a larger spatial scale. Our findings are discussed in relation to the fisheries exploitation patterns, the environmental features of the area and the biological characteristics of the species.

Sampling
Data were collected in the framework of the MED-ITS bottom-trawl surveys accomplished in the period 1994-2015 at 15 GSAs of the Mediterranean (Fig. 1). Sampling included annual sampling at pre-defined stations following a standardized protocol (see Bertrand et al. 2002, Spedicato et al. 2019. Trawling was performed by means of the experimental net GOC-73 with a 20-mm square mesh cod-end, and selection of stations was based on a depth-stratified sampling scheme that included five depth strata: 10-50, 51-100, 101-200, 201-500 and 501-800 m. Collected data included number, weight, gonad maturation stage and total length measurements for a wide range of fish, cephalopod and crustacean species. In the present study we analysed biomass data for red and striped red mullets by station, expressed in terms of kg per square km of swept area (kg km -2 ).

Data analysis
The analysis utilized generalized additive modelling (GAM) techniques (Hastie and Tibshirani 1990) to investigate the effects of station position and depth on species biomass index. The main advantage of GAM over traditional regression methods is its capability to model non-linear relationships (a common feature of many ecological datasets) between a response variable and multiple explanatory variables using non-parametric smoothers. In the present case, the non-linear predictors included sampling position (entered as the latitude-longitude interaction), depth and year of sampling. A preliminary analysis indicated that the presence of red mullets at depths over 400 m was negligible, so the stations below that depth were not considered in the analysis. The final data set included a total of 16716 hauls.
Given the existence of stations with no presence of red mullets, and based on the diagnostic randomized quantile residual plots (Foster and Branvington 2013), Tweedie GAM models following a compound Poissongamma approach were applied. This avoids multiplestage modelling of zero-inflated data and allows the probability of presence and the non-zero sampled quantity to be modelled jointly (Shono 2008, Lecomte et al. 2013. The smoother function used for the year and depth effects was a penalized cubic regression spline and model fitting was accomplished using the mgcv library (Wood 2006) under the R language environment. The procedure automatically selects the degree of smoothing based on the generalized cross-validation score, which is a proxy for the model's predictive performance. However, to avoid dubious relationships, the model was constrained to be at maximum a quartic relationship. Hence, the maximum degrees of freedom for each smoothing term, measured as number of knots (k), was set to 4 (i.e. k =5 in the GAM formulation). For the longitude-latitude interactions the smoother function used was thin plate regression splines, which are low-rank isotropic smoothers without knots, which provide a sensible way of modelling interaction terms, avoiding problems associated with knot placement (Wood 2003). GAM models assumed a log link function and were fitted for each GSA separately. Deficiencies of the fitted models were diagnosed by means of randomized quantile residual plots (Foster and Branvington 2013).
In the next step, dynamic factor analysis (DFA, Zuur et al. 2003) was applied to the estimated year effects from the GAM models, in order to identify common time series trends among GSAs. The method is based on fitting and comparing models assuming a linear combination of N common trends that are modelled as smoothing curves by means of a Kalman filter/smoother algorithm. In the present case, models containing 1 to 3 common trends were tested and the "best" candidate model was defined based on the Akaike information criterion (AIC), which provides a balance between the model fit and the parameters used. The model with the smallest AIC value was selected. Model fit was accomplished under the R language environment using the MARSS library (Holmes et al. 2012). As the data were previously standardized, an "identity" structure for the error covariance matrix R was assumed. Subsequently, hierarchical cluster analysis was used in order to facilitate identification of common GSA groups based on the factor loadings from the common trends estimated from the DFA. The Euclidean distance was used to represent dissimilarity among GSAs and clustering was done using the average linkage method. All statistical inferences were based on the 95% confidence level.

M. barbatus
The highest values of mean standardized biomass (kg km -2 ) were observed in Corsica (GSA 8) and Crete (GSA 23) and the lowest ones in the most western areas (Table 1).
The GAM models applied to the biomass indexes explained a large proportion of the total variance (>55% in the majority of the GSAs), while depth and latitudelongitude interactions were always significant. The effect of year was significant in all GSAs apart from 8 and 20 (Supplementary material Table S1). Variations were observed among GSAs regarding the estimated temporal biomass pattern, but in most GSAs of the central and eastern part of the Mediterranean, increasing trends were observed from 2008 onward. A decreasing trend was observed in the Balearic Islands (Fig. 2). Generally, biomass was found to decrease with depth, and this was more evident at depths below 200 m (Fig. 3).
DFA model fitting to the 15 time series was accomplished assuming 1 to 3 common trends (CT). The AIC values for the one-, two-and three-CT models were 668.39, 666.75 and 686.48, respectively. Thus, the  Broken lines indicate two standard errors and the relative density of data points is shown by the "rug" on the x-axis. model with two CTs was considered the best solution.
The first CT showed a pattern including a sequence of sharp decreases/increases. Specifically, biomass drops were observed up to 2000, followed by an increasing trend until 2008 and a decrease afterwards. (Fig. 4).
The second CT showed a continuous increasing trend throughout the period examined. The fitted values suggested that most time series were fitted well to the model (Supplementary material Fig. S1). The factor loadings for the DFA model indicated which CT was related to which GSA time series (Fig. 4). The scores indicated roughly that the first CT was important for GSAs 1, 5, 22, 23 and 25, while the second was related to GSAs 7, 8, 9-11 and 18-20. GSAs 10 and 17 were influenced almost equally by both CTs (opposite sign for the first CT) but with low relevance, given the low scores (GSA 10 in particular). Cluster analysis revealed two major GSA groups, each of them including minor subgroups (Fig. 5). All groups were delimited by the geographical position of the GSAs, i.e. the major groups separated eastern to western GSAs, while subgroups included nearby GSAs.

M. surmuletus
The highest values of mean standardized biomass (kg km -2 ) were observed in GSAs including islands, with the Balearics (GSA 5) having by far the highest rate (Table 2). Other GSAs with relatively high values were Crete (GSA 23), the Aegean Sea (GSA 22), Corsica (GSA 8) and Sardinia (GSA 11).
The GAM models applied explained in most cases more than 30% of the biomass variance, and depth and station position (longitude-latitude interactions) were always significant. The effect of year was significant in only five of the GSAs (5, 11 and 16-18), located mostly in the central Mediterranean (Supplementary material  Table S2). In these cases, a declining trend was observed in the latest years (Fig. 6). Generally, biomass was found to decrease progressively with depth, but in several GSAs this was more evident at depths below 250 to 300 m (Fig. 7).
DFA model fitting to the 15 time series was accomplished assuming 1 to 3 CTs. The AIC values for the one-, two-and three-CT models were 723.78, 712.84 and 718.15, respectively. Thus, the model with two CTs was considered the best solution. The first CT showed a continuous increasing pattern throughout the period examined, while the second was rather dome-shaped, with the peak in 2003-2004 (Fig. 8). The fitted values suggest that most time series were fitted well to the model (Supplementary material Fig. S2). The factor loadings scores indicate roughly that the first CT is more relevant for GSAs 8, 10, 19, 22 and 23 (opposite sign), while the second one is mostly related to GSAs 5, 7, 17 and 25 (opposite sign). GSAs 11, 16 and 18 are influenced almost equally by both CTs, while GSA 20 has very low factor loadings (not shown in Fig. 8) suggesting negligible relevance with the identified CTs.
Hierarchical cluster analysis applied on the factor loadings revealed two major groups, each of them in- cluding minor subgroups (Fig. 9). Similarly to M. barbatus, all groups were delimited by the geographical position of the GSAs, i.e. the major groups separated eastern to western GSAs, while subgroups included nearby GSAs.

DISCUSSION
The common sampling protocols followed by the MEDITS surveys give the opportunity of performing Mediterranean-wide studies of the species distribution patterns following harmonized analysis, which facilitates the detection of large-scale events. In the current work, the use of GAM techniques for computing standardized biomass rates removed effects that could bias nominal indexes (Hilborn and Walters 1992) and allowed direct comparisons among areas and identification of CTs. Utilizing a Tweedie distribution model, which has been reported to perform statistically better than a delta-type two-step model (Shono 2008), also allowed simultaneous analysis of zero and positive values and facilitated interpretation of the results.
As the MEDITS surveys are accomplished at the end or just after the reproduction period of red mullets, full recruitment of the newly born individuals to the surveyed fishing grounds cannot be traced. In fact, the MEDITS surveys are accomplished during the May-July period, while recruitment occurs from August to October (Suau andVives 1957, Reñones et al. 1995). For this reason, the identification of spatio-temporal differences through the analysis of a biomass index was preferred, since such an index is less sensitive than one expressed in terms of numbers to spatio-temporal fluctuations owing to the irregular presence of recruits.
Our results indicate that biomass of both species shows clear declining trends after 150-to 200 m depth, and this decline is sharper in the case of M. barbatus. This finding is in line with previous studies suggesting that M. surmuletus has a wider bathymetric range than its congeneric species (Lombarte et al. 2000). In addition, it has been found that M. barbatus was more abundant than M. surmuletus in all examined areas except the Balearic Islands, where the opposite situation occurs. However, the abundance of M. surmuletus may have been underestimated, owing to its preference for rough bottoms not always accessible to the bottomtrawl gear.
For both species, differences in the spatial distribution pattern can be attributed, at least to a certain extent, to the different biotic and abiotic conditions prevailing in each area. Past studies have suggested that certain geo-morphological and seafloor characteristics affect the distribution of red mullets. M. surmuletus prefers rough substrates and narrow continental shelves, while M. barbatus is more abundant on muddy bottoms and areas where the continental shelf becomes wider (Hureau 1986, Fischer et al. 1987, Lombarte et al. 2000. In addition, eco-morphological studies suggest the existence of adaptive morphological and anatomical characteristics that allow M. barbatus to exploit better than its congeneric species resources from muddy and turbid bottoms (Lombarte and Aguirre 1997). M. surmuletus seems to have more developed hearing sensitivity (Lombarte 1992, Aguirre andLombarte 1999), which would give greater adaptation to relatively noisy environments, such as reefs or inshore waters.
Regarding temporal biomass patterns, increases were observed for M. barbatus from 2008 onward in several GSAs. Such increases are in line with recent assessment studies suggesting that several Mediterranean red mullet stocks are in healthy condition, or at least in much better situation than stocks of species inhabiting deeper waters (Vasilakopoulos et al. 2014, Tserpes et al. 2016, STECF 2016, Cardinale and Scarcella 2017. It is likely that the implementation of Council Regulation (EC) No 1967/2006, which introduced additional trawling prohibitions in coastal areas, has contributed to this increase (Tserpes et al. 2011).
Given that recruitment of both species mainly occurs at bottoms very close to the coast, at depths ranging from 10 to 50 m (Suau andVives 1957, Levi et al. 2003), the establishment of satellite-based vessel monitoring systems through Council Regulation (EC) No 1224/2000 has also favoured such increases by discouraging illegal fishing operations in coastal areas. Furthermore, the mesh-size specifications introduced through the aforementioned regulation improved size selectivity of the bottom-trawl fishery by increasing the length at first capture of M. barbatus by 30% to 50% (Sala et al. 2015), thus contributing to overall biomass increases of the species, and most likely to the number of larger spawners.
Several past studies have identified a positive effect of female age/size on larva survival, concluding that the maternal effect on egg quality is very important for recruitment success (Trippel et al. 1997, Vallin and Nissling 2000, Nazari et al. 2009). In line with these findings, a recent study on the reproductive potential of M. barbatus in the southern Adriatic Sea showed that fish length was positively correlated with egg size and plasmatic concentration of vitellogenin, suggesting that survival rates of larvae coming from bigger/older females are higher than those hatched from eggs pro- duced by smaller/younger individuals (Carbonara et al. 2015). Thus, it is expected that abundance increases of larger spawners, due to selectivity improvements, would lead to overall biomass increases of the species. The temporal biomass pattern is different for M. surmuletus, which, showed stability in the majority of the GSAs, similarly to what has been reported in a past study that analysed a shorter time series of MEDITS data in various Mediterranean areas . Given that M. surmuletus also inhabits rough bottoms not exploited by trawlers but accessible to other fishing gears, it is likely that the aforementioned trawling prohibitions have less positive impact on it than on M. barbatus. In the Balearic Islands, the only area where M. surmuletus shows very clear decreasing trends throughout the period examined, the species is widely distributed in rocky and gravel bottoms, which predominate on the shelf of the archipelago, and it is an important target species for the fisheries operating in the area (Quetglas et al. 2012).
Apart from changes in the fisheries exploitation pattern, environmental changes may have favoured the increasing trends observed in certain areas. In the Aegean Sea for instance, increasing abundance trends for several demersal species have been partially attributed to recruitment increases owing to increased productivity (Tserpes and Peristeraki 2002). Certainly, the intensity of fishing activities always plays an important role, and this may explain why areas of low fishing pressure, such as GSA 23 (Crete), show relatively high abundances (Table 1) although they are characterized by low primary production (Peristeraki et al. 2017). Levi et al. (2003), suggested that positive sea surface temperature (SST) anomalies (warmer waters than average) could enhance recruitment of M. barbatus, probably through a reduced upwelling regime, which-although it may involve lower productivitywould result in lower offshore transport of early fish stages, a key factor for the survival of juvenile stages. Maravelias et al. (2007) reported that red mullet abundance in the Aegean Sea (eastern Mediterranean) was consistently higher in areas with warmer waters near the seabed. Fiorentino et al. (2008) suggested that SST increases in the waters of northwestern Sicily during the pre-recruitment phase of red mullet could have a positive effect on the stock, contributing to the benefits derived from trawling bans. Under such scenarios, the warming trends of the Mediterranean waters (e.g. Nykjaer 2009, Vargas-Yáñez et al. 2008) may have favoured the recruitment of M. barbatus in several areas. Although there is a lack of relevant physiological studies, the differences in the spatial distribution pattern of the species (Kaschner et al. 2016) suggest that M. barbatus is more thermophilic than its congeneric. Consequently, sea-water warming may have contributed more to population increases of M. barbatus than to those of M. surmuletus.
Although it is reasonable to assume that the combined effects of fishing and environmental conditions determine species abundance variations, the relative importance of each component may vary among areas. Examination of stock dynamics in areas with different environmental characteristics and various levels of fishing pressure would help to clarify the magnitude of each component, as well as the way they interact (Rouyer et al. 2008).
DFA revealed similarities among neighbouring GSAs in terms of temporal abundance trends, and the subsequent cluster analysis identified two major GSA groups corresponding to the eastern and western part of the Mediterranean basin (Figs 5 and 9). As both species have a pelagic larval stage, larval dispersal may, at least to a certain extent, explain similarities in abundance trends among nearby areas. In general, the spatial connectivity among fish populations has been traditionally attributed to larval dispersal and adult movements (Grüss et al. 2011). Given the low mobility of red mullets and the fact that their eggs and larvae are pelagic, with larvae being able to remain 25-35 days in the water column until they find a suitable habitat to settle (Macpherson and Raventos 2006), it can be assumed that larval dispersal and survival determined by environmental conditions are of critical importance in the context of the spatial biomass differences observed in the current study. In line with this, Gargano et al. (2017) reported that the connectivity between M. barbatus population sub-units inhabiting the European and African shelves of the Strait of Sicily was low, and high quantities of larvae were transported outside the study area. Based on the biological features of the species and the patterns of the sea-water currents, the authors suggested that larvae and pre-recruits lost from the Strait of Sicily can reach suitable coastal areas to recruit even far from their natal area, supporting the existence of a meta-population structure, sensu Kerr and Goethel (2014), of red mullet in the Strait of Sicily and the adjacent areas.
Genetic studies give a more complex pattern regarding spatial connectivity of red mullet populations in the Mediterranean. Maggio et al. (2009), working on six microsatellite loci of M. barbatus from 14 different sites of the central Mediterranean, reported the isolation of the Adriatic samples and a weaker substructuring of the populations in the Gulf of Lions, the Tyrrhenian Sea, the Strait of Sicily and the Ionian Sea. Recently, Matić-Skoko et al. (2018) working on 13 microsatellite loci from 15 Mediterranean sites found differences in the genotype of M. barbatus inhabiting the eastern and western coast of the southern Adriatic Sea, while no differences were reported for the northern and middle Adriatic populations. Surprisingly, no gene barrier was observed between the Balearic Sea and the eastern Mediterranean, while the Tyrrhenian Sea specimens were grouped together with the northern and middle Adriatic populations. In the case of M. surmuletus a lower level of genetic heterogeneity was found between populations. The authors evidenced a clear sign of genetic bottleneck in both species that was attributed to demographic changes caused by a combination of high fishing pressure, habitat fragmentation and naturally occurring fluctuations in population size.
Discrepancies between population dynamics and genetic studies in questions related to the spatial con-nectivity of populations are most likely due to the different temporal scale underlying the two approaches. While results of genetic studies reflect dynamics at the evolutionary scale, classical population dynamics depict variations at the ecological scale, which is shorter and limited to rather recent events.
Our findings regarding spatial association could be used for management purposes to assist the identification of population units that behave as a single stock with the same population dynamics (Cope and Punt 2009). Spatial connectivity, however, depends on numerous parameters and interdisciplinary studies are necessary to provide robust conclusions on the structure of the stocks and support the rational management of fishery resources.

SUPPLEMENTARY MATERIAL
The following supplementary material is available through the online version of this article and at the following link: http://scimar.icm.csic.es/scimar/supplm/sm04888esm.pdf