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 (Gaertner et al. 2010, 2013, 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 (Gaertner et al. 2010, 2013, 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.

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).

We investigated the stable patterns and spatio-temporal 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, STATICO-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), heteroptera (Slimani et al. 2017), macrophytes (Le Fur et al. 2019), pelagic fishes and top predators (Certain et al. 2011), pelagic larval and juvenile fishes (Carassou and Ponton 2007), zooplanctons (Mazzocchi et al. 2012), as well as indicators of sustainability (Gallego-Alvarez et al. 2015).

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 and Chessel 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² 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 and Chessel 1994, Dray et al. 2003). Co-inertia analysis is a two-table 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² 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 and Kidé 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).

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² values, Fig. 2) showed both a common and a specific part of the assemblage organization in each of the seventeen GSA, with cos² 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.

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%).

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, phytoplancton 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 species-environment 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).

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).