Anthropogenic impacts, and especially the impacts of climate change and fisheries, are topics of great interest among marine scientists, as they can potentially affect the structure and functioning of marine ecosystems (e.g. Odum 1985, Pauly et al. 1998). Sea warming caused by the greenhouse effect has been evidenced in many seas, including the Mediterranean Sea (Belkin 2009, Skliris et al. 2011, Macias et al. 2013), and it is believed that it causes changes in the marine community structure resulting from changes in the expansion (Cheung et al. 2010, Tsikliras and Stergiou 2014, Tsikliras et al. 2015) or in the abundance of thermopile marine species (Vasilakopoulos et al. 2014, Givan et al. 2018). In addition, the impact of fisheries on the marine environment has been evidenced by numerous authors (e.g. Dayton 1995). Recent studies, mainly using ecological models, have revealed that ecological degradation is often reflected in changes in the marine food webs and particularly in a shift from higher trophic level predators to species at the base of the trophic pyramid, mainly due to the harvesting of species of higher tropic levels by the world fisheries (e.g. Pauly et al. 1998). However, other authors question the existence of a general shift to lower tropic levels in the marine ecosystems, and in some cases even document the opposite trend (Caddy et al. 1998, Branch et al. 2010, Branch 2012).

In the present study we examined the spatiotemporal trends in two biodiversity indicators in several Mediterranean areas based on standardized time series derived from the Mediterranean International Trawl Survey (MEDITS). We attempted to test two working hypotheses related to the impact of sea warming and fishing on the demersal communities of the Mediterranean Sea. The first hypothesis is based on a study by Cheung et al. (2013), who use the term “mean temperature of the catch” (MTC), defined as the average inferred temperature preference of the exploited species weighted by their annual catch, to evaluate the effect of sea warming on fishery catches and marine ecosystems. According to this hypothesis, the global mean temperature of fisheries catches increases positively related to the increase in sea surface temperature as a result of climate change. Past studies have examined the effects of sea warming on the distribution and abundance of marine fish populations in various areas (e.g. Cheung et al. 2010, 2013, Tzanatos et al. 2014), but few have focused on the impact of sea warming on the community structure in the Mediterranean (Keskin and Pauly 2014, Tsikliras and Stergiou 2014, Tsikliras et al. 2015), and most of them are based mainly on landings statistics. A similarly constructed index, the “mean trophic level of the catches” (MTL) is used to test the “fishing down the food web” hypothesis, formulated by Pauly et al. (1998), for the demersal community of the Mediterranean Sea. The authors, based on the decrease in the MTL of global fisheries landings between 1950 and 1994 suggest that the MTL in marine ecosystems is decreasing due to the overexploitation of large fish/predators belonging to a higher tropic level. Several past studies have examined the trophic habits, trophic relations and tropic level of individual species in the Mediterranean Sea (for a review see Karachle and Stergiou 2017). Ecopath models developed for Mediterranean areas address the topic of food web properties but, apart from a recent study examining the spatiotemporal changes of functional-group biomass (Brind’Amour et al. 2016), studies of spatiotemporal changes in the trophic structure of demersal megafauna community on a broad Mediterranean level are lacking.

The present study was based on data obtained from the MEDITS experimental surveys conducted in the northern Mediterranean Sea. The MEDITS survey is carried out annually following a standardized protocol in various areas of the Mediterranean Sea, and its primary goal is to monitor changes in the abundance of demersal megafauna species. The survey design is based on a depth-stratified random sampling scheme that considers five depth strata: 10-50, 50-100, 100-200, 200-500 and 500-800 m. At each station (haul), the total number and biomass by species were estimated (see Bertrand et al. 2002 and the MEDITS Handbook [Anonymous 2017] for further details on survey design and sampling methodology). Data on species biomass by haul used in the present analysis cover all surveys carried out during the 1994-2015 period in 13 Geographic Statistical Areas (GSAs) of the Mediterranean Sea (5, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20, 22, 23), as determined by General Fisheries Commission for the Mediterranean (Fig. 1). In GSAs 5, 20, 22 and 23 the survey was not conducted every year, leaving gaps in the data time series.

The annual standardized biomass index by species (Souplet 1996), expressed in kg per square km of swept area (kg km^–2), was used to calculate the MTC and MTL of the catch per GSA per year. Species with a negligible biomass index (representing in total 1% of the total biomass per GSA per year) were excluded from the analysis. The preferred temperature of each species was acquired from the supplementary material available in Cheung et al. (2013). Species lacking this information were excluded from the MTC analysis. The tropic level of each species was obtained from several dietary studies of the Mediterranean Sea (Valls et al. 2014, Albo-Puigserver et al. 2016, Karachle and Stergiou 2017). Only species with an available tropic level were included in the analysis. The MTC and MTL per GSA per year were estimated by averaging the corresponding species rates weighted by the species biomass indexes. These indexes was preferred against the total catch, as they are standardized to the same unit (km²) and account for minor changes in the sampling (slight differences in the number or location of hauls, etc.) during the sampling period. The species included in the MTC and MTL analyses accounted for more than 76% and 78% of the total catch in terms of biomass, respectively.

The estimated MTC values for each GSA were examined for temporal autocorrelation by means of the R function acf. The same was also accomplished for MTL. The computed autocorrelations were plotted against lag in years and no significant correlations were detected. Given the lack of temporal correlation, a generalized additive modelling (GAM) approach (Hastie and Tibshirani 1990) was used to examine the effects of year (by GSA) and GSA on MTC and MTL separately. Hence, two GAM models were built, having the general form:

Index = a + s(Year, by = GSA) + GSA + e

where Index is either MTC or MTL, a is the intercept, s indicates the smoothing function of the corresponding independent variable and e is a random error term. The smoothing function used for the continuous variable (year) 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 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 was set to 4 (i.e. k=5 in the GAM formulation).

A Gaussian distribution model accompanied by an identity link function were used, as preliminary analysis of the diagnostic residual plots of several models indicated that they provided a pertinent fit to the data.

In order to be able to contrast the spatial MTC pattern with temperature variations, data on bottom temperatures per haul recorded during the MEDITS surveys, collected by temperature sensors attached to the net, were subsequently analysed. After a preliminary analysis of the diagnostic residual plots of the models, the effects of year, depth and GSA on bottom temperature were examined by means of a GAM assuming a Gamma distribution model and a log link function. Thus the form of the GAM model was:

Bottom temperature = a + s(Year) + s(Depth) + GSA + e

where a is the intercept, s indicates the smoothing function of the corresponding independent variable and e is a random error term. Due to important data gaps, interactions among the predictor variables, which would provide information on temperature variations at the GSA level, were not included in the model. The smoothing functions used for the continuous variables (year, depth) were similar to those previously mentioned.

Both GAM models applied to the MTC and MTL indexes explained a high proportion of the total variance (>67%). Differences among GSAs were significant. The models also revealed differences among GSAs with respect to the year effect (Table 1). For MTC, the year was significant only in the Gulf of Lions, northern Adriatic, western Ionian and Cretan seas (GSAs 7, 17, 19, 23), while for MTL it was significant in Sardinia, the Strait of Sicily, and the northern and southern Adriatic and Cretan seas (GSAs 11, 16, 17, 18, 23). The effect of year (when significant) on MTC and MTL is presented in Figures 2 and 3, respectively. The MTC and MTL estimates by GSA are presented in Figures S1 and S2 of the Supplementary Material. In most cases there were no great differences among years.

In the Gulf of Lions (GSA 7) and western Ionian (GSA 19), MTC showed a slight decreasing trend (Fig. 2). In the northern Adriatic (GSA 17) and Cretan seas (GSA 23), MTC showed a decreasing trend during the first years of the time-series, but in GSA 17 it increased in the last two years, while in GSA 23 an increase was observed until 2010 and a decrease thereafter (Fig. 2). However, the decrease observed during recent years in GSA 23 may be an artefact due to the lack of data between 2009 and 2013 in this area. For the area effect, a general increasing trend of MTC was observed from west to east and from north to south. The highest MTC values were observed in the Cretan and western Ionian seas and the Strait of Sicily (GSAs 23, 19, 16, respectively) and the lowest in the Gulf of Lions, Corsica and northern Adriatic seas (GSAs 7, 8, 17) (Fig. 4).

Regarding the significant effect of year on MTL, variations were also observed among areas. In Sardinia and the Strait of Sicily (GSAs 11, 16) an increasing trend was observed, while decreasing trends were observed for the northern and southern Adriatic (GSAs 17, 18). In the Cretan Sea (GSA 23), a decreasing trend was observed until 2005, followed by a slight increase (Fig. 3). For the area effect, the highest MTL values were observed in the Gulf of Lions and Ligurian Sea (GSA 7, 9), followed by the values estimated in Sardinia, the Strait of Sicily and the southern Adriatic (GSAs 11, 16, 18). The lowest values were observed in Corsica, followed by the southern and central Tyrrhenian and the northern Adriatic (GSAs 8, 10, 17) (Fig. 5).

The analysis of deviance of the GAM model applied to the recorded bottom temperatures by haul indicated that the effect of GSA, year and depth was always significant (Table 2). Bottom temperature showed a prominent increasing trend after 2010 (Fig. 6), decreasing with depth down to 200 m and showing small variations beyond that (Fig. 6). Finally, a clearly increasing temperature trend was observed from west to east and from north to south (Fig. 7).

The Mediterranean Sea is a semi-enclosed marine area characterized by a diversity of environmental conditions and a high diversity of marine organisms (Durrieu de Madron et al. 2011). Due to its size and limited water exchange at the Gibraltar Strait, Mediterranean Sea dynamics are mainly linked to the local climate and are therefore particularly sensitive to climatic/anthropogenic perturbations. Such perturbations are caused by global warming, fishing activities, and other human interventions such as the opening of the Suez Canal and the construction of the Aswan Dam, which led to serious changes in the environment and the biodiversity of the Levantine basin (Galil 2006).

Intense fishing pressure and overexploitation of fish stocks in the Mediterranean has been reported by Botsford et al. (1997) and Vasilakopoulos et al. (2014). Mediterranean Sea warming has also been documented by numerous studies. Rixen et al. (2005) found an increase of about 0.5°C in the average temperature of the upper layer (0-150 m) of the Mediterranean Sea during the years 1980-2000, Belkin (2009) estimated a much higher rise of 1.4°C from 1978 to 2003, and Skliris et al. (2011) found an annual warming rate of 0.037°C per year for the entire Mediterranean over the 1985-2008 period, with higher estimates for the eastern basin (0.042°C per year). The results of the current study indicate a slower increase in the bottom temperature than in the surface and upper water layers. However, the observed spatial bottom temperature trends are similar to the mean sea surface temperature and salinity trends, showing increases along a longitudinal gradient from west to east, and along a latitudinal gradient from north to south (MyOcean 2014). A similar trend has been found for primary production (Bosc et al. 2004, Moutin and Raimbault 2001, Ignatiades et al. 2009). The present spatial analysis of the MTC values revealed a similar pattern. For the spatial trends of MTL, the highest values were observed in the northern-western and northern-central areas, which also undergo higher fishing pressure than the eastern Mediterranean (Vasilakopoulos et al. 2014). Lower MTL values were observed in areas with lower fishing pressure, such as Corsica and Crete (STECF 2015, Sylaios et al. 2010). The above findings disagree with the “fishing down food web” theory, which assumes that intense fishing results in depletion of larger predators and subsequent dominance of low-trophic-level species (Pauly et al. 1998). Apart from fishing pressure, primary productivity also affects the marine community structure (Piroddi et al. 2017), and combinations of those parameters may explain the spatial variation of the observed MTL values. Nevertheless, further information on the spatiotemporal patterns of fishing effort on a broad Mediterranean scale is needed to identify the effects of fishing on the community structure. The simple use of landing rates as a fishing pressure indicator ignores the importance of biotic and environmental factors on stock productivity. Also, recent changes in fisheries legislation (spatial closures, temporal fishery restrictions etc.) may affect the species composition of the landings.

For MTC and MTL, in the vast majority of areas examined no significant temporal trends were found, suggesting that fishing pressure and sea warming did not generally affect demersal communities during the study period (1994-2016). In only two of the studied areas was a temporal increasing trend of MTC observed, in line with past findings in various areas (Cheung et al. 2013, Tsikliras et al. 2015). However, past studies in the Mediterranean were based on landings statistics, largely comprised of pelagic species that are more susceptible to temperature changes (Alheit et al. 2014). By contrast, the present study takes advantage of a relatively long (22 years) fishery-independent data series obtained through a standardized sampling protocol specifically designed for demersal communities, and is thus more representative of demersal communities’ structures. Rather than the MTC increase expected according to the hypothesis of Cheung et al. (2013), MTC estimates showed a decreasing trend in two of the study areas.

For MTL, significant decreasing trends were observed only in 2 out of the 13 areas studied, while in 2 others increasing trends were observed. These results conflict with the findings of past and recent studies that were based on landings statistics (e.g. Pauly et al. 1998, Keskin and Pauly 2018). However, because landings data are susceptible to market demand, discarding practices and differentiation of fishing strategies (Caddy et al. 1998, Branch et al. 2010, Branch 2012), it is likely that they are not fully representative of the marine communities and may result in biased estimates.

In the few areas which show significant temporal trends in the MTC and MTL values, it is worth examining the underlying causes of these changes. The increasing MTC trends in the Cretan Sea may be due to the particularly high increasing rates of sea temperature (from 0 to 300 m depth) observed in the eastern Mediterranean since 1985 (Skliris et al. 2011, Soto-Navarro and Criado-Aldeanueva 2012). These increases favour the expansion of alien thermophile marine species of Indo-Pacific origin, resulting in changes in the marine faunal composition (Golani et al. 2002). This is particularly evident in the case of Crete, as has been demonstrated by various studies (e.g. Kasapides et al. 2007, Peristeraki et al. 2007, 2017). In the same area, the significant temporal variations in MTL values could be, at least partly, attributed to changes in the fishing pattern. The implementation of additional management measures in the last decade, such as the extension of the banned area to trawlers from 1 to 1.5 miles from the coast (EC Regulation 1967/2006) and additional temporal closures after 2012 (Greek Management Plan 2012), may have a serious impact on the coastal marine communities in the area, as Crete has an extremely narrow continental self and a steep slope (Lykousis et al. 2002) with a restricted trawlable area (Tserpes et al. 2011). Thus, it is likely that a drastic reduction of the fishing effort on the continental self has occurred in the recent years. However, lack of data in this area for the period 2009-2013 increases the uncertainty of model predictions for the more recent years.

In the northern Adriatic, MTC values showed an increasing trend after 2010, similarly to the bottom temperature pattern, while decreasing trends of MTL values were observed in both the northern and southern Adriatic. This response of the northern Adriatic demersal communities to sea warming may be attributed to temperature increases in the Mediterranean upper layer observed since the 1980s (Rixen et al. 2005, Belkin 2009). Due to the prevailing shallow depths, these increases may have a relatively direct impact on the benthic communities of the area. The decreasing trend in the MTL values in both the northern and southern Adriatic can be attributed to the increase in low-trophic-level species in the survey catches, such as M. barbatus, B. boops, E. encrasicolus and S. pilchardus. The decreasing MTC trends in the Gulf of Lions and the western Ionian are inconsistent with the hypothesis of Cheung et al. (2013). Taking into account that bottom temperatures in the Gulf of Lions were the lowest of the examined western areas and those in the western Ionian were the lowest of the southeastern areas, it could be speculated that, due to the general sea warming, these areas act as a “refuge” for psychrophilic species inhabiting adjacent areas with higher temperatures. Hence, they show increasing abundances of psychrophilic species. A similar hypothesis has been put forward by Lasram et al. (2010), but only as a prediction for the remote future, i.e. the middle of the 21st century. Examining the data sets used in the present analysis, we found that, indeed, psychrophilic species showed higher biomass indexes in recent years in the aforementioned areas. More specifically, in the Gulf of Lions, the psychrophilic species Sprattus sprattus showed a huge increase in its biomass index after 2009, and Engraulis engrasicolus has also shown generally higher biomass indexes during recent years. Though small pelagic species are usually not fully available to the bottom trawl, their increase in the survey catches may indicate a high concentration of these species in the aforementioned areas or a differentiation of their distribution in the water column, possibly due to sea warming. Additionally, in the case of the western Ionian, psychrophilic species such as Phycis blennoides showed biomass increases from 2005 onwards, while thermophilic species such as Lepidopus caudatus have decreased in the recent years. Certainly, a detailed investigation of the spatiotemporal abundance changes at the species level would better clarify the species-specific changes in the community structure of these areas. Further, Lasram et al. (2010) highlighted the importance of considering habitat fragmentation (e.g. semi-isolated basins, islands) in the Mediterranean Sea when exploring species’ spatial responses to climate change.

In Sardinia and the Strait of Sicily, a general increasing trend of the MTL values was observed. In the Strait of Sicily this finding could be explained by the combined effect of two main driving forces: a drastic decrease in fishing effort on the demersal resources followed by an increase in the density and biomass of many elasmobranch species. Since 2005, both fleet capacity and fishing effort of the trawl fleet operating in the Strait of Sicily has shown a slow, progressive decline (Mannini and Sabatella 2015). During the same period, some taxonomic groups have shown an increase in density and biomass indexes. During the study period, there was a significant increase in the abundance of cartilaginous fish (Gancitano et al. 2015, STECF 2013), mainly represented by top predators of a high tropic level. In our data set, the elasmobranch species Galeus melastomus, Raja clavata and Squalus blainvillei, usually abundant in trawl catches, show increasing biomass indexes during recent years. Similarly, in Sardinia, a decrease in the number of small trawlers, which mainly exploit the shallow resources, occurred during the first decade of the study period. Additionally, the total number of fishing boats has been progressively decreasing during the last two decades (Mannini and Sabatella 2015). These conditions could have positively affected the marine communities in the area and explain the observed increasing trend of MTL, in the light of the theory of Pauly et al. (1998).

The absence of temporal MTC trends in most of the areas suggests that the recent bottom temperature increases have not caused an immediate response to the structure of demersal marine communities, so it can be deduced that there is a delay in the response of the demersal communities to sea warming. In areas where there is already a prolonged temperature increase in the upper and intermediate waters (e.g. Crete), a shift in the demersal communities regime to more thermophile species has occurred. The same is valid for relatively shallow areas (e.g. the northern Adriatic).

Nevertheless, sea warming is a phenomenon in evolution, so a constant effort is needed to investigate its impact on the marine environment. Areas more susceptible to sea-warming effects could be advantageous for scientific monitoring aimed at identifying the impact of the greenhouse effect on benthic communities.

The lack of significant temporal changes in MTL rates in most of the areas suggested the absence of any general shift in the Mediterranean demersal food webs during the study period. The current findings indicate that MTL trends based on the landings of commercial fisheries may diverge from ecosystem MTL trends obtained from surveys and assessments.

Although temporal MTL trends could be related to fishing pressure changes in a few areas, the present results do not show a clear connection between MTL and fishing pressure, so they do not generally support the validity of the “fishing down marine food web” theory in the Mediterranean.

The variability among GSAs regarding the observed MTC and MTL temporal patterns reflects differences in both environmental conditions and exploitation patterns that have implications for the ecology and dynamics of the stocks.