Stability of the relationships among demersal fish assemblages and environmental-trawling drivers at large spatio-temporal scales in the northern Mediterranean Sea and 25 years of surveys

Summary: Trawling pressure and environmental changes may affect the composition of fish assemblages. Our knowledge on large spatio-temporal patterns of demersal fish composition remains incomplete for the Mediterranean Sea. We investigated (1) the spatio-temporal stability of demersal assemblages, (2) the relationships between these assemblages and potential structuring factors (trawling pressure and environmental conditions) in order to assess the dynamic of the assemblage structure at the scale of the northern Mediterranean Sea. We analysed a dataset of 18062 hauls from 10 to 800 m depth performed annually during the last two decades across 17 Geographical Sub-Areas (GSAs) (MEDITS program). A multi-table analysis (STATICO-CoA) evidenced a strong inter-GSAs stability in the organization of assemblages, with specificities for some GSAs. The most stable structuring factors were linked to combined gradients of chlorophyll a , phytoplancton carbon biomass and temperature, inversely correlated with depth, salinity and nutrient gradients (axis 1 of the STATICO-CoA compromise, 93.74% of the total variability). A common pattern linking the distribution of species to these environmental gradients was evidenced for most of the 17 GSAs. Estimate of trawling pressure showed a minor role in the organization of the assemblages for the spatial scale and years investigated (axis 2, 4.67%).


INTRODUCTION
Due to its geological evolution and environmental conditions, the Mediterranean Sea is characterized by relatively high biodiversity (Coll et al. 2010). However, it is also one of the world's most threatened marine regions due to environmental change and fishing activities (Ben Rais Lasram et al. 2010, Coll et al. 2010. These drivers might affect the structure and composition of species assemblages (Menge andOlson 1990, Greenstreet andHall 1996), in particular those of demersal fish targeted by intensive commercial bottom trawling, such as in the northern Mediterranean Sea (Rochet et al. 2010). However, our knowledge on the impact of trawling and environmental forcing on demersal assemblages in the Mediterranean is still largely incomplete. On the one hand, previous studies analyzing fish composition at large scale in this sea have provided insights based on gridded species presence-absence maps drawn on the basis of expertbased knowledge (e.g. Mouillot et al. 2011, Coll et al. 2012). However, these works were not necessarily focused on demersal fish assemblages caught by trawl-ing. In addition, they were mainly restricted to indirect presence-absence data estimated at coarse spatial grain and they could not investigate temporal trends. On the other hand, some previous studies investigated large spatio-temporal trends in the Mediterranean Sea on the basis of fish abundance data collected by standardized scientific surveys. However, these studies were mainly focused on diversity components, such as species number (Gaertner et al. 2007) evenness, taxonomic, phylogenetic and/or functional diversity , Granger et al. 2015, Brind'Amour et al. 2016), but they did not deal with the structure and composition of demersal fish assemblages per se. Indeed they mainly focused on diversity indices which thus summarize in several single and synthetic values some features of assemblages, while some multi-variate analyses allows to assess the identity and contribution of species to patterns.
Thus, the aim of this study is to complete our knowledge in this domain by assessing large-scale spatial and temporal patterns of the structure and composition of demersal fish assemblages sampled by trawling over the northern Mediterranean Sea. More specifically, the Summary: Trawling pressure and environmental changes may affect the composition of fish assemblages. Our knowledge on large spatio-temporal patterns of demersal fish composition remains incomplete for the Mediterranean Sea. We investigated (1) the spatio-temporal stability of demersal assemblages, (2) the relationships between these assemblages and potential structuring factors (trawling pressure and environmental conditions) in order to assess the dynamic of the assemblage structure at the scale of the northern Mediterranean Sea. We analysed a dataset of 18062 hauls from 10 to 800 m depth performed annually during the last two decades across 17 Geographical Sub-Areas (GSAs) (MEDITS program). A multi-table analysis (STATICO-CoA) evidenced a strong inter-GSAs stability in the organization of assemblages, with specificities for some GSAs. The most stable structuring factors were linked to combined gradients of chlorophyll a, phytoplancton carbon biomass and temperature, inversely correlated with depth, salinity and nutrient gradients (axis 1 of the STATICO-CoA compromise, 93.74% of the total variability). A common pattern linking the distribution of species to these environmental gradients was evidenced for most of the 17 GSAs. Estimate of trawling pressure showed a minor role in the organization of the assemblages for the spatial scale and years investigated (axis 2, 4.67%). main objectives are to assess both the stable and variable parts of the spatio-temporal relationships between the species composition of demersal fish assemblages and the environmental-trawling variables at large scale in the Mediterranean Sea. The analyses were based on species abundance data collected from 1999 to 2015 within a standardized scientific protocol consisting in 18062 sampled hauls performed along more than 20000 km of coastline over a wide bathymetric zone (10-800 m), seagrass excluded (Posidonia oceanica meadows, etc.).

Survey methods and faunal data collection
Data was collected from annual bottom trawl surveys performed in May-July from 1999 to 2015 over the continental shelf (10 to 200 m depth) and the continental slope (200 to 800 m) of the Mediterranean Sea within the MEDITS scientific program (Bertrand et al. 2002a, MEDITS Working Group 2017. The study area ranged from 34.33°N to 45.67°N and 5.22°W to 34.09°E and was divided into 17 operative Geographical Sub-Areas (GSAs), for which boundaries were established by the General Fisheries Commission for the Mediterranean Sea (management areas, according to the resolution GFCM/33/2009/2, http://www.gfcm. org, see Fig. 1 for the correspondence between GSA numbers and their names, and Table S1 for number of hauls per GSA).
The sampling procedures were standardized according to a common protocol over GSAs and years. The standard device is a bottom trawl GOC-73 with 20-mm cod-end mesh size. The average vertical opening of the gear is 2 m and its wing-span 18 m (Bertrand et al. 2002b). All the tows were performed during daylight hours. Speed on the ground and duration of the tows were standardized to 3 knots and 30 min for shelf stations. For the slope this length of time was doubled (60 min) in order to overcome potential difficulties in the deep hauls (see Bertrand et al. 2002a for full details). For each haul, we collected raw fish abundance data, bathymetry, swept area (i.e. trawled surface) and haul coordinates.
Information was recorded by an underwater Scanmar system, to control the trawl geometry (horizontal and vertical openings, contact with the bottom), as well as to record swept area and trawling time, which enabled us to select the tows that were properly carried out. Among a total of 394 fish species sampled during MEDITS surveys (Relini 2015), we considered 154 species efficiently sampled by the bottom trawl system used, and properly identified by teams involved in this large scale program , Granger et al. 2015. Indeed, i) only demersal fishes (i.e. benthic and bentho-pelagic fishes) were considered, ii) individuals for which the identification at species level was complex were gathered at genus level (it was only the case for individuals of Lophius spp., see Table S2). This subset was thus designed in such a way as to minimize potential bias caused by gear-selectivity (e.g. Lefcheck et al. 2014), and to strictly limit the risk of a variability of accuracy in sampling identification between the different teams. The list of 154 species considered was set according to those used in previous MEDITS data based studies , Granger et al. 2015; see Table S2 for the list of species considered). Raw abundance data was normalized by the swept area to compute species densities.

Environmental and trawling data
A total of eight explanatory variables were considered: one related to trawling pressure and seven to environment. Among the seven environmental variables, the mean depth of operated haul during surveys was computed. The six other variables (i.e. temperature, salinity, chlorophyll-a concentration, phosphate, nitrate and phytoplankton carbon biomass) were extracted from the Copernicus portal (http://marine.copernicus. eu/) when available for the years of the MEDITS survey (i.e. 1999 to 2015, thus leading to the exclusion of hauls sampled from 1994 to 1998). According to the variable, two models were considered to deliver their values. Both models have a horizontal grid resolution of 1/16×1/16° with 72 unevenly spaced depth levels (Oddo et al. 2009). The values of the variables were extracted according to the depth level that embeds the mean depth of the haul performed. Temperature (°C) and salinity (psu) were derived from the 3D Mediterranean Sea Physics Reanalysis model, whereas chlorophyll-a concentration (mg m -3 ), nutrients (phosphate and nitrate (mmol m -3 ) and phytoplankton carbon biomass (mmol m -3 ) were computed from the 3D Mediterranean Sea Biogeochemistry Reanalysis model (details available at http://marine.copernicus.eu/). Two values for each variable were prepared for each haul: the monthly-mean value over the year before the haul had been performed, and its standard deviation (sd). During a preliminary analysis, we investigated the pairwise redundancy/complementarity of the eight explanatory variables (mean and sd), by using a Draftsman's plot and Spearman's correlation coefficient. Because we found that mean and sd values were correlated, only means were used as explanatory variables in the statistical analyses presented below.
Given that data on fishing vessel locations (e.g. VMS/AIS data) were not available for all spatial and temporal scales of the MEDITS survey, a Multi-Criteria Decision Analysis (MCDA) has been employed to estimate a trawling pressure index for commercial bottom trawlers. As a rule, MCDA provides successive complementary techniques and procedures for structuring decision problems and deriving a way for a compromise in a transparent process (Malczewski 2006). In our case, a GIS-MCDA process that combines influential geographical components affecting bottom trawling in terms of its distribution and intensity (such as fishing capacity, bathymetry etc. see more details in Fig. S1), as well as expert judgement opinions (Table S3), have been used based on a methodological approach proposed by Kavadas et al. (2015). The method was applied to bottom trawlers classified by total length (LOA) segmentation: <12 m, 12-24 m and 24-40 m. The generated trawling pressure index was calculated as the average of trawling pressure indices per LOA on a spatial resolution of 0.01×0.01 decimal degrees. The index resulted from the fuzzy product of an activity index and a suitability index (see details in Fig. S1).
The activity index was based on trawling vessel activities estimated by registration port for 1994, 2004 and 2014, considering the vessel length, gross tonnage and the number of vessels (source: European Fleet Register). An Inverse Distance Weighted (IDW) interpolation method was used to spatialize this index. Then a standardization with a linear Fuzzy Membership function was applied to scale the index between 0-1.
The suitability index reflected the suitability of trawling activity by combining information on several criteria: bathymetry (source: EMODNET), distance from the coast (estimated by ArcGIS ESRI, 2011), annual chlorophyll a (https://oceancolor.gsfc.nasa.gov) and fisheries restricted areas (legislation, source: MEDISEH, Colloca et al. 2013). Each criterion was assigned a grading value by expert judgment (i.e. a rank of order of importance from 0 to 5, Table S3).
The setting of this suitability index from the investigated criteria proceeded in the following steps: (i) creation of spatial information and calibration of each criterion according to a scale of evaluation and formation of the hierarchical structure of the multiple criteria problem; (ii) implementation of the Analytic Hierarchy Process (AHP, Saaty 1980) to estimate the relative importance of the evaluation criteria (Table S4); (iii) application of the Weighted Linear Combination method using the weights (priority vectors) to estimate the suitability index (Table S4); (iv) standardization on a scale from 0 to 1 with linear Fuzzy Membership function (see Kavadas et al. 2015 for details of the overall methodology).
Finally, a value of the trawling pressure index was assigned to each MEDITS operated haul according to its spatial location (defined as mean position) and year. Trawling pressure index values of 1994, 2004 and 2014 were assigned to hauls performed in 1994 to 2003, 2004 to 2013 and 2014 to 2015, respectively. It is thus assumed that the index remained constant over each of these periods. Note that even if generally the index value remained relatively similar for a given location between 1994, 2004 and 2014, it varied in time within some GSAs between these periods (Figs S2A and S2B). In addition, the index showed contrasted spatial patterns among GSAs (Fig. S2A).

Statistical analysis
We investigated the stable patterns and spatiotemporal variability of the relationships between species composition of demersal fish assemblages and the environmental/trawling variables over the period under study at the scale of the northern Mediterranean Sea. To that end, we carried out a multi-table analysis, STA-TICO-CoA (Simier et al. 1999, Thioulouse et al. 2004, Thioulouse 2011) that enabled us to investigate both the stable patterns and the spatio-temporal changes of these relationships. STATICO is a multivariate, multi-table method that allows to access the relationships between two kinds of information. Here, data were organized in a series of pairs of species and environmental-trawling tables associated with each of the 17 GSAs. STATICO, a co-inertia analysis, was used to assess not solely the relationship between species distribution and explanatory variables but also the degree of spatial-temporal reproducibility of this relationship. It is thus a very useful method to investigate key aims in ecology with data of large spatio-temporal scale compared to traditional analysis to study this relationship (Mendes et al. 2012, Le Fur et al. 2019). STATICO has been successfully applied on fauna and flora of various ecosystems, such as demersal fishes (Kidé et al. 2015), diatoms (Mendes et al. 2012), estuarine fishes (Simier et al. 2006 More in details, as recalled by Kidé et al. (2015), STATICO (Simier et al. 1999) is an application of the STATIS method (Lavit et al. 1994) to co-inertia operators (Dolédec andChessel 1994, Dray et al. 2003). In other words, STATICO "combines the objectives of STATIS (finding the stable part of the structure of a series of tables) and the objectives of co-inertia analysis (finding the common structure of two data tables)", as stated by Thioulouse et al. 2004 (see also Kidé et al. (2015) for a reminder of the description of the vectorial approach of STATICO). In comparison with the classical version of STATICO, STATICO-CoA (Thioulouse et al. 2004) used in the present study is based on the logic of correspondence analysis (CoA, based on the Chi 2 distance Hill 1973) to investigate the organization of species assemblages (Le Fur et al. 2019, or Gaertner et al. 2002 for a STATIS-CoA).
The data analyzed are a sequence of pairs of tables with the environmental/trawling variables and, separately, the species density, for each GSA. Each table is first analyzed using a CoA for the species tables, and a PCA for the explanatory tables. Species densities (n) were log (n+1) transformed to reduce the influence of dominant species. Environment-trawling data were centered and reduced (i.e. normed) in order to consider their different units. Each pair of tables is then linked by a co-inertia analysis (Dolédec andChessel 1994, Dray et al. 2003). Co-inertia analysis is a twotable coupling method, which allows a cross-table to be computed between the variables of the two tables (here between species and environmental/trawling variables) for each GSA. The resulting series of species and environmental/trawling variables cross-tables is then analyzed by the STATIS method (Lavit et al. 1994). This analysis is based on three main steps, each of them offering different kinds of results (see a STATICO flow chart in Kidé et al. 2015): (1) the inter-structure step identifying the proximity between each of the K tables (here the 17 tables assessing the link between species composition and potential drivers for each of the GSAs); this step aims to compute the weighs of each of these K tables for the construction of a mean  table of maximum inertia (compromise table). In other words, weights are used to build a linear combination of the series of tables, called «the compromise». These weights are the components of the first eigenvector of the RV coefficients matrix (see Lavit et al. 1994). The fit of each of the tables to the compromise table is given by the Cos 2 value (absolute value ranges between 0 and 1, 1 being the best fit); (2) the compromise analysis, that defines factorial axes which express the stable part of the structures of the K tables, giving an ordination of the reproducible part of the relationships between species composition and environment/trawling that are common for the 17 GSAs across the years. If a/some GSAs contribute less or poorly to the compromise, it means that environmental conditions/trawling pressure were different and/or acted differently on assemblage composition than in the other GSAs; (3) the trajectories step, where species and environmental/trawling variables for each year can be projected as additional elements on the compromise axes in order to summarize their variability around the common structure. For the sake of clarity, a clustering analysis was added at step (2) for identifying on the compromise assemblages of species (see Mazzocchi et al. 2012 andKidé et al. 2015 for a similar approach). In this way, groups of species were identified with a hierarchical classification based on Euclidean matrix of pairwise distances between species on the first factorial plane (i.e. 1 st and 2 nd axes) of the compromise analysis. The UPGMA aggregation criterion, chosen on the basis of an objective approach among the main available aggregation criterion (see Mérigot et al. 2010), was used to compute this classification. The optimum number of species groups was identified with the R package NbClust that computes and compares thirty indices to identify the optimum number of groups (Charrad et al. 2014). All statistical analyses were performed using R 3.2.2 software (R Core Team 2015).

Reproducible part in the relationships of assemblages and environmental/trawling variables among GSAs
The common structure associated with the compromise (first axis of the inter-structure, Fig. S3) explains 52.51% of the total variability. In the construction of the compromise table, greater importance is given to the tables which have similar structures and limited importance to the other tables. Here, we found that the contributions (weights) of the GSAs to the construction of the compromise were relatively similar, ranging mainly from 0.19 to 0.29, except GSA 17 (Northern Adriatic Sea) weighing 0.14 ( Fig. 2), with a standard deviation of 0.04. It means that the nature of the link between species assemblages and environmental/trawling drivers for GSA 17 will be less taken into account in the construction of the compromise table. Moreover, the fit of each of the tables to the compromise (see cos 2 values, Fig. 2) showed both a common and a specific part of the assemblage organization in each of the seventeen GSA, with cos 2 ranging mainly between 0.61 and 0.87 (except GSA 17 and 25 (Cyprus), with 0.38 and 0.56, respectively), with a standard deviation of 0.12. These weaker values observed for GSA 17 and 25 indicated that the structure of these GSAs contributed less to the construction of the compromise table. These results are confirmed by the RV correlation coefficients values among GSAs (Table S5).
The use of a separate Correspondence Analysis (CoA) carried out independently in each GSA and the projection of the first axes of these separate CoAs on the first two axes of the STATICO-CoA compromise confirmed these results (Fig. S6). Except for GSA 17, the first axis of the compromise appeared strongly correlated to the first axis of each of the separate CoAs (Fig. S6). This result showed that the main assemblage structure was strongly similar in each of the sixteen GSA. Projection of the axis 2 of each separate CoA showed a stronger variability in the organization of the assemblages in each GSA, especially for GSA 17, 22 and 25. In short, these results (Figs 2 and S6) suggested a strong inter-GSA stability in the relationships between demersal species composition and environmental/trawling variables in the northern Mediterranean Sea. It also showed that a limited part of this relationship was specific to some GSAs. Along these lines, GSA 17 (Northern Adriatic Sea) and 25 (Cyprus) appear to be the most specific of the seventeen studied GSAs.

Common pattern in the environmental/trawling variables among GSAs
The first two axes of the compromise accounted for 98.41% of the total variability of the common structure on environmental/trawling variables and species with an overwhelming 93.74% for axis 1 (Figs 3 and  4). Thus, the first axis alone very well summarizes the stable part of the link between species composition and environmental/trawling drivers for all studied GSAs (although it is more limited for GSA 17, see Fig. S4). It shows that the dominant pattern of the environmental variables associated with the demersal assemblages variability in the different studied areas consisted in a chlorophyll a, phytoplancton carbon biomass and temperature gradient, inversely correlated with depth, salinity and nutrient gradients (see axis 1 in Fig. 3; pairwise plot and correlation between variables are provided in Fig. S5). The estimate of trawling pressure index only played a role on axis 2 (Fig. 3) which took into account a very small part of the compromise variability (4.67%).

Common pattern in species composition among GSAs
The projection of the species onto the first two axes of the compromise allowed identification of the species mainly contributing to the common structure of the demersal assemblages at the scale of the northern Mediterranean Sea (Fig. 4). It gave a typology of the species showing stable trends over the years for the 17 GSAs. There was a clear gradient in the distribution of species (Fig. 4). Two main fish assemblages were identified on the basis of the hierarchical clustering method and optimum number of groups approach computed on the species coordinates on the two first axes of the compromise analysis (Fig. 4, see Materials and Methods section).
As mentioned above, most of the variability is caught in the first axis, representing 93.74% of the total variability. A first group of species was associated with negative values of coordinates on the left part of this axis (i.e. higher temperature, phytoplanc-   Table S2. ton carbon biomass and chlorophyll a), while a second group was located on the positive and right part of axis 1 (i.e. higher depth, nutrients concentration and salinity , Figs 3 and 4). The relevance of analysing axis 2 which represented 4.67% of the variability is limited. Only some species of the two main assemblages could be distributed according to the estimate of trawling pressure (upper left part) and depth (bottom right part). More precisely, species associated with higher chlorophyll a, phytoplancton carbon biomass and temperature (i.e. left bottom part of the factorial plane, Fig. 3) were mainly Bothus podas, Synodus saurus, Arnoglossus imperialis, Dasyatis pastinaca and Trachinus radiatus (Fig. 4). Species collected in zones where depth, salinity and nutrient concentration were higher (i.e. right bottom part of the factorial plane, Fig. 3) included Alepocephalus rostratus, Centrophorus uyato, Heptranchias perlo, Bellottia apoda and Paralepis coregonoides (Fig. 4).
For each GSA, the projection on the compromise axes of the explanatory variables and species is provided in detail in Figure S4. It shows which GSAs contributed the most to the stable part of the speciesenvironment dynamics (Figs 3 and 4) in the compromise analysis. Except for GSA 17, the stable part of the environmental/trawling variables on axis 1 of the compromise was clearly found for each separate GSA (i.e. similar pattern of explanatory variables in Figures  3 and S4). In contrast, projections of species showed more variability according to the GSAs considered (Figs 3 and S6).

Trajectories of environmental-trawling variables and species assemblages
The main trends of the temporal variability of environmental-trawling conditions and demersal assemblages are provided by their respective trajectories on the first two axes of the compromise for each GSA (Fig. 5). A year point corresponds to the mean coordinates on each axis of hauls performed during that year. The interpretation of the trajectories plots is done by analyzing in each GSA the contraction, the stretching of the cloud of points, and the proximity or remoteness along the axes. All these elements provide information on the relationships between the environment (Fig. 3) and fish assemblage variables (Fig. 4) along the studied years in a given GSA (Fig. 5). Temporal similarity in the trajectory patterns of fish assemblages and environment/trawling suggests that changes in assemblages are related to the environmental-trawling conditions. Overall, the temporal trajectories between environment, trawling features and assemblages did not match along the years analyzed, i.e. there was an absence of a clear coupling in time between environment and trawling features and species trajectories (Fig. 5,Table S6). A link is highlighted from the Kendall correlation coefficient τ computed among the projections of environmental/ trawling variables and species with τ value superior to 0.6 only on x-axis (93.74% of the total variability of the common structure) for GSA 1 (τ x =0.75) and GSA 20 (τ x =0.78) ( Table S6). Fig. 4. -Projections of demersal species variables and assemblage groups on the compromise of the STATICO-CoA analysis. It shows the stable part of the species and environment-trawling relationships. Axis 1 and 2 explained 93.74% and 4.67% of the total variability, respectively. On the first factor plane are shown the two main species assemblages (in blue and red) obtained by hierarchical clustering (UPGMA criterion) and optimum number of groups approach (see Materials and Methods section). Species codes are given in Table S2.

DISCUSSION
The results of the STATICO-CoA analysis suggested a strong inter-GSAs stability in the spatio-temporal relationships between fish assemblage composition and environmental conditions/trawling pressure of the demersal assemblages in the northern Mediterranean Sea (Figs 2 and S4). The analysis highlighted that the most stable structuring variables (explaining over 94% of the total variability) were linked to a combined gradient of chlorophyll a, phytoplancton, carbon biomass and temperature, inversely correlated with the depth, salinity and nutrient gradients. Estimate of trawling pressure contributed relatively poorly to the common organization of the assemblages at the spatial scale and in the years investigated (less than 5% of the total variability). A clear gradient in the distribution of species (Fig. 4) linked to these environmental gradients was evidenced as the stable part of the species environment relationships (note that similar results were obtained for analyses performed without rare species, i.e. those present in less than 5% of the hauls, Gaertner et al. 2002, Figs S7 and S8). However, none of the large majority of GSAs considered showed temporal similarity between the temporal trajectories of the environment (including trawling) and the trajectories of fish assemblages, suggesting an absence of a clear coupling in time between environment, trawling pressure and species trajectories (Fig. 5, Table S5). However, our results also showed specificities in the relationships between species composition and environmental/trawling variables in GSA 1 and 20 (Fig. 5, Table S6), as well as GSA 17 Northern Adriatic Sea and GSA 25 Cyprus (Figs 2 and S3, Table S5).
The relative importance of abiotic and biotic factors in shaping assemblage patterns is considered to vary with spatial scale (Menge and Olson 1990). For instance, at a local scale, physical and biotic factors may interact to influence local patterns of assemblage structure. Larger spatial scales are considered to be associated with an increase in the relative influence of variation in environmental or climatic conditions (Menge and Olson 1990). Our results highlight the importance of environmental conditions with regard to species assemblage composition and structure at the scale of the northern Mediterranean Sea. Using environmental and trawling pressure information describing each hauls for the different GSAs, we provided insights on the determinism of the assemblage structure at this large scale. To our knowledge, there are few quantitative studies either in the Mediterranean Sea or worldwide that investigate simultaneously the effects of both environmental and anthropogenic factors, with a focus on the species composition, structure and abundance of exploited demersal assemblages (both target and non-target species) at large spatio-temporal scales. Firstly, in the Mediterranean Sea, studies performed on large-scale spatial and/or temporal distribution of fish species assessed mainly species diversity (e.g. species richness, evenness, taxonomic, phylogenetic and/or functional diversity) (Mouillot et al. 2011, Granger et al. 2015, Brind'Amour et al. 2016. In contrast, there has been less focus on assemblage structure and mainly at a lower spatial scale, i.e. between the Straits of Gibraltar and the Gulf of Lions (Gaertner et al. 2005). Previous studies based on MEDITS data focusing on species diversity did not find an effect of a longitudinal gradient on diversity (Gaertner et al. 2007, Granger et al. 2015, Brind'Amour et al. 2016, in contrast to results provided by studies based on other types of data (presence-absence, coarser spatial grain, etc. see Coll et al. 2010, Mouillot et al. 2011. Notably, using standardized data, Granger et al. (2015) found that species richness showed relatively low values in the Adriatic Sea and Cyprus at both local and regional scales (i.e. alpha and gamma diversity, respectively), similarly to specificities we evidenced in assemblage organization in these two GSAs. They also highlighted that diversity indices remained stable over the last two decades. This result is in line with our findings that trawling pressure, estimated by the MCDA, had a low effect on assemblage composition at the scale of our study. This does not mean that intensive trawling, such as that practised in numerous GSAs of the northern Mediterranean Sea for more than 50 years, had no effect on the benthic ecosystems and the associated fauna (e.g. D'Onghia et al. 2003, Farriols et al. 2017. It may suggest that most of the impact of trawling pressure had probably altered the demersal fish diversity and composition before the beginning of the studied period (1994 for the MEDIT surveys) (Granger et al. 2015, Farriols et al. 2017. Otherwise, possible uncertainty in estimating the trawling pressure spatialization in our study, due to the lack of spatial fishing effort data (such as from VMS data, not available at the spatio-temporal scale of this study), cannot be ruled out. Thus, to be able to establish more explicit links with fishery management perspective, the availability of finer spatio-temporal fishing effort data, such as VMS data, should be broaden and increased for the different GSAs. It may provide more direct assessment of the relationship between demersal assemblages and trawling pressure (e.g. Farriols et al. 2017), and thus facilitate suggestions towards management.
Secondly, when considering other regions of the world, a recent study was based on two decades of scientific trawl surveys focused on exploited demersal fish assemblages in Mauritania (eastern Atlantic Ocean) among different depth strata and latitudinal areas (Kidé et al. 2015). Similarly to our results, they highlighted a stronger effect of environmental conditions than trawling on the structure and composition of demersal fish assemblages. They found that chlorophyll a and sea surface temperature mainly influenced assemblages in areas where upwellings occurred. Effects of trawling pressure on the structure and composition of assemblages were relatively low. Similarly to our results, the temporal trajectories between environmental and trawling conditions and assemblages did not match over the entire time series they analyzed (only for some specific years and areas).
However, it should be noted that in other regions, targeted species were affected more directly by fishing and led to changes in the structure of demersal assemblages (e.g. Ansari et al. 1995, Greenstreet and Hall 1996, Levin et al. 2006. Regarding temporal trajectories among the environmental conditions and the species, the fact that for the large majority of GSAs a mismatch was observed across years (Fig. 5,Table S6) could be linked to (i) lags in response of assemblages in the face of changes in environmental conditions, (ii) approximate way that environmental values were computed (i.e. monthly mean across the year before the haul was performed) and/or (iii) the fact that assemblages could be sustained despite environmental changes. It may also suggest that factors other than those we have investigated may act on fish assemblages, such as biotic interactions.
In conclusion, our study provides a state of the spatio-temporal relationships between organizational patterns of demersal fish assemblages and several environmental and trawling descriptors, on the basis of a large standardized data set of scientific trawl samplings. Our analysis represents a step towards quantitative studies focused on the drivers of demersal fish diversity in the Mediterranean Sea. It suggests a strong inter-zone stability in the organization of fish assemblages, with specificities for some zones. The most stable structuring factors were linked to a combined gradient of chlorophyll a, phytoplancton carbon biomass and temperature, inversely correlated with a depth, nutrient and salinity gradient. A clear gradient in the distribution of species linked to these environmental gradients was evidenced as the stable part of their relationships. Investigating explicitly potential biotic factors (such as species aggregation, competition or dispersal limitations) would extend our findings within the framework of complementary analyses dedicated to this aspect.

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/sm04954esm.pdf Table S1. -Number (n) of hauls analyzed per Geographical Sub-Areas (GSA). N=18062 hauls sampled between 1999 and 2015 within the 17 GSAs (see Fig.1 for map of locations and GSA names).  Figure 1. Note that for this kind of multi-tables analysis, only axis 1 of the inter-structure has a meaning for the construction of the compromise, and is thus used to determine the contribution of tables to the compromise. The other axes are not used for that purpose, and should not be interpreted (axis 2 is used for display purpose). the compromise of the STATICO-CoA analysis based on a data set of 71 species (without 83 rare species, i.e. those present in less than 5% of the hauls). It shows the stable part of the species-environment-trawling relationships. Axis 1 and 2 explained 94.63% and 4.43% of the total variability, respectively. Species names are provided in Table S2. Fig. S8. -Projections of demersal species variables and assemblage groups on the compromise of the STATICO-CoA analysis based on data set of 71 species (without 83 rare species, i.e. those present in less than 5% of the hauls). It shows the stable part of the species and environment-fishing relationships. Axis 1 and 2 explained 94.63% and 4.43% of the total variability, respectively. On the first factor plane are shown the two main species assemblages (in blue and red) obtained by hierarchical clustering (UPGMA criterion) and optimum number of groups approach (see Materials and Methods section). Species codes are given in Table S2.    Figure 1. Note that for this kind of multi-tables analysis, only axis 1 of the inter-structure has a meaning for the construction of the compromise, and is thus used to determine the contribution of tables to the compromise. The other axes are not used for that purpose, and should not be interpreted (axis 2 is used for display purpose).   species (without 83 rare species, i.e. those present in less than 5% of the hauls). It shows the stable part of the species-environment-trawling relationships. Axis 1 and 2 explained 94.63% and 4.43% of the total variability, respectively. Species names are provided in Table S2. Fig. S8. -Projections of demersal species variables and assemblage groups on the compromise of the STATICO-CoA analysis based on data set of 71 species (without 83 rare species, i.e. those present in less than 5% of the hauls). It shows the stable part of the species and environment-fishing relationships. Axis 1 and 2 explained 94.63% and 4.43% of the total variability, respectively. On the first factor plane are shown the two main species assemblages (in blue and red) obtained by hierarchical clustering (UPGMA criterion) and optimum number of groups approach (see Materials and Methods section). Species codes are given in Table S2.