Samples were collected during the Mediterranean International Trawl Survey (MEDITS) (Bertrand et al., 2002) between July and August 2018 at nine sites ( Figure 1 ; Table 1 ). We examined an area historically exploited by fisheries: the Strait of Sicily (SoS). The SoS is a biologically important area of the central Mediterranean Sea characterised by high habitat complexity and biodiversity (Di Lorenzo et al., 2018; Russo et al., 2019), which notably suffers from ecological alterations (Agnetta et al., 2019; Lauria et al., 2020).

Sampling sites covered three depth strata (10–50, 51–100, and 101–200 m), three distance ranges from the coast (< 10 Km, 10 Km-20 Km, > 20 Km), and different levels of fishing effort intensity (Russo et al., 2019; De Angelis et al., 2020). According to the EUNIS classification, all sampling sites are classified as having a substratum predominantly consisting of muddy sand. All the details about the sampling and the methodological procedures applied are available in Russo et al. (2021).

Here it is important to recall that, through the amplification of two mitochondrial regions (a ~167 bp fragment of the 12S gene using the Teleo02 primer and a ~313 bp fragment of the COI gene using the universal metazoan Leray-XT primer) (Taberlet et al., 2018; Wangensteen et al., 2018), a list of species was obtained for each of the nine sites considered using the eDNA material in the “slush” (the dense water draining from fishing nets just after the end of hauling operations), together with a corresponding list of species visually identified in the catch.

The list of taxa identified by 12S and COI were merged by sampling site and the Multi Response Permutation Procedure (MRPP) (Berry and Mielke, 1994) implemented in the R package vegan (Oksanen et al., 2018) was then applied to assess whether there was a significant difference between the replicates grouped by sampling sites. MRPP is not constrained to the assumptions of MANOVA (e.g. Normality and Equal Variances for each variable) and allows to compare, through a permutational approach, the inter-sites dissimilarities (i.e. between replicates of different sites) with the intra-site dissimilarities (i.e. between replicates of the same site) in terms of species composition.

To investigate how food-web structure varies across the study area, we collected the available information on the trophic relationships of all the species identified through both eDNA metabarcoding and catch data by conducting an extensive literature review (diet and feeding studies obtained from Fishbase, SeaLifeBase, and Google scholar). Since food webs reflect the trophic functioning of the ecological community, the network position of each species (node) is informative about its trophic role. The node position into the network may be a proxy of functional importance and provide macroscopic descriptors of ecosystems. In this analysis, we considered adult feeding links only and excluded feeding interactions involving juveniles, larvae or eggs. This decision was partially justified by the temporal scale of the sampling, which occurred during summer, some months after the spawning periods of most of the species considered. However, juveniles and larvae are often not considered in this kind of modelling exercise since the scant amount of available information for these stages is likely to impair the food web reconstruction (Frelat et al., 2022). The process of eDNA metabarcoding taxon aggregation assumes that there are nodes in the network that are taxonomically similar enough that we can consider them functionally equivalent. Taxonomic aggregation was performed along the hierarchy genus and species level. In this step, each set of congeneric species was lumped into a single functional group, producing a higher genus-level node for the food webs, following Agnetta et al. (2019) ( Supplementary Table 1 ). These aggregation criteria of species may not always reflect the ecological reality of the group; however, we defined the functional groups by combining the taxonomical similarity and the bathymetric distribution of the sampled demersal species, considering these as the main drivers in the definition of functional groups. Following Agnetta et al. (2019), at the lower trophic level of the local food webs we placed primary consumers: Algae, Phytoplankton, and Seagrasses; together with basic resources including: Benthic detritus, Pelagic and Sediment bacteria, Suspended particulate organic matter (SOM). This information was used to reconstruct the predator-prey relationships among all the taxa identified through the MENA. Then, the R package cheddar (Hudson et al., 2013; Hudson, 2022) was used to generate two models (one based on the catch data and one based on the eDNA metabarcoding data) of the food web and to estimate the trophic level (TL_cheddar) of each node (taxon). Then, the consistency and reliability of these estimates TL_cheddar were evaluated by comparing them with the trophic level information extracted from a number of widely used and popular databases (TL_reference) (see Supplementary Table 2 for details). Values of TL_reference were retrieved for each species, attempting to filter the different values contained in the databases so as to select those obtained through studies conducted in the Mediterranean. If more than one value was available for a given species, the mean was calculated. A linear regression analysis was then used to assess the agreement between the TL_cheddar predicted from heuristic food webs (from both catch and eDNA data) and the respective values of the TL_reference from the reference databases (Compson et al., 2019). A delta was computed as the difference between TL_cheddar and TL_reference and the distributions of their values were inspected for catch and eDNA data. In addition, this analysis was used to evaluate the two parallel reconstruction of the food webs (from both catch and eDNA data) and identify the most reliable and robust one.

For illustrative purposes, we reconstructed nine local food-webs combining the information from the replicates of each site. Nodes (taxa) at each site were defined as the subset of taxa identified in at least two of the replicates of that site. The R package igraph (Nepusz, 2022) was then used to produce a graphical representation of the local food webs.

A set of four ecological indicators was computed for each replicate, including: (1) Species Richness (N); (2) Number of links (L); (3) Direct connectance (C); (4) Generality (G) ( Table 2 ). These indicators were used to summarise the changes in food web structure (topological characteristics) and network complexity for each local food web. Then, GAM (Hastie and Tibshirani, 1990) were used to fit the relationships between each food web indicator and three predictors: sea bottom depth, distance from the coast, and fishing effort intensity exerted from trawlers. Nonparametric regression approaches such as GAMs enable the modeling of variable connections without requiring the underlying regression function to take on a certain form. The use of smooth functions as regressors gives GAMs greater flexibility over parametric (including linear) types of models. GAMs were computed using the R package mgcv (Wood, 2022). The model’s goodness of fit was assessed using the coefficient of determination (adjusted R²) and the amount of deviance explained. Residuals for each model were graphically inspected to detect any departures from the model assumptions or other possible anomalies in the data. All statistical analyses were conducted in the programming environment R ver, 1.4.1103 (www.rproject.org). According to the results (section below) of the comparison between estimated (TL_cheddar) and reference (TL_reference) trophic levels, which was used to compare reconstructions obtained from capture data and those obtained from eDNA data (see the subsection above), the procedure described in this subsection was applied only on the eDNA-based food webs.

The detailed information returned by the literature review and from the two databases FishBase (Froese and Pauly, 2022) and SeaLifeBase (Palomares, 2022) about species interactions (i.e. feeding relationships or trophic links) are available in Supplementary Table 3 ( Supplementary Materials ).

Overall, the comparison between the set of TL_cheddar estimated using cheddar (and thus based on eDNA metabarcoding) and TL_reference extracted from popular databases ( Figure 2 ) showed a significant agreement for both the reconstructions, although the values estimated in the model based on eDNA metabarcoding data were much closer to the values known from the literature ( Figure 2-A1 ) than the values obtained for the catch-based data ( Figure 2-A2 ). However, for some species (labelled points in Figure 2-A1 ), the estimated TL was slightly different from the reference values.

The distribution of the differences (delta) between estimated and reference TL, when inspected by class, showed that Chondrichthyes and Actinopterygii were the groups with the greatest variances in the eDNA-based model, while for some other classes (e.g benthic malacostraca and benthic cephalopoda) the variability was much lower ( Figures 2-B1, C1 ). The corresponding distribution of the difference between baseline and estimated TL values for the model based on the capture data was much more platykurtic ( Figure 2-B2 ). The distribution of the delta by class ( Figure 2-C2 ) showed that Chondrichthyes and Actinopterygii (as for the eDNA), but also Echinoidea and Scyphozoa, were the groups characterised by the more variables.

Overall, these results indicate that the eDNA-based MENAs perform better than the corresponding catch-based ones, placing the taxa in trophic positions consistent with previous knowledge. Accordingly, and for the sake of conciseness, the following results are focused on the analysis of the eDNA-based MENA.

Environmental DNA Metabarcoding allowed for the identification of 101 taxa, of which 7 taxa of basal resources, 82 at species level, nine at genus level, one at family level, one at superorder level and one at subclass level ( Table 3 ). The N values ranged from 59 to 71, while the L values ranged from 306 to 382, spanning over 4 TL. Apical species (TL ~ 4.5) included the silver scabbardfish Lepidopus caudatus, the European hake Merluccius merluccius, the common dentex Dentex dentex, the marbled electric ray Torpedo marmorata, and the John Dory Zeus faber. Species with slightly lower TL were the red scorpionfish Scorpaena scrofa, the European conger Conger conger, the whiting Merlangius merlangus, the bluntnose spiny dogfish Squalus megalops, and the European angler Lophius piscatorius. Other predators with TL = 4 were the common smooth-hound Mustelus mustelus, the Blainville’s dogfish Squalus blainville, the stargazer Uranoscopus scaber, and the blue whiting Micromesistius poutassou ( Table 3 ).

The MRPP established that the differences (in terms of species co-occurrence) between samples belonging to the same site (mean value of the Jaccard distance between replicates: 0.451) were significantly less than those between samples belonging to different sites (mean value of the Jaccard distance between sites: 0.743), supporting statistical consistency among replicates at each sampling site. Network analysis revealed considerable variability in food webs without clear patterns. The Local food-web structure varied across space in relation to the intensity of fishing effort and the considered abiotic variables (distance from the coast, depth layers) ( Figure 3 ; Table 1 ). The most complex food webs, in terms of N, L and G were represented by the food webs A1, B12 and C5, characterised by intermediate and low intensities of fishing effort, respectively. In terms of ecological complexity represented by C; the direct connectance varied inversely with N and L (see the three categories of fishing effort intensity Intermedium, low and no fishing effort in Table 1 ). Apical predators present in all the local food-webs included D. dentex (ID = 99) an active and solitary predator, as well as C. conger (ID = 93), M. poutassou (ID = 91), and M. merluccius (ID = 98) except for the local food web B5. As it may be expected, the most striking differences occurred at the top of the webs. Local food webs with fewer top predators corresponded to intermediate depths and distances (Local food webs B5 and B19, Table 1 ; Figure 3 ), while the level of fishing effort did not seem to have a noticeable effect on this aspect. However, there were no immediate visual differences among the local food webs at different locations. The food webs indicators of N, L and C varied substantially between sites, with the widest range of values occurring at intermediate values of depth (50-100m) and distance from the coast (10-20 km) (N = 58 at site B19, N = 71 at site A1; L = 279 at site B19 and L = 371 at site A1). C ranged between 0.074 (minimum value at site C18 and 0.086 (maximum value) at site B5. G ranges between 4.823 (minimum value) at site B19 to 5.354 (maximum value) at site A1 ( Table 3 ). The generated GAMs appeared, in almost all the cases, well-trained (with Deviance explained > 67% and adjusted R² > 0.60) ( Figure 4 ; Table 4 ). In general, both the environmental gradients and the fishing effort had significant effects on all the indicators. N was significantly correlated with all three predictors and showed increasing trends with respect to depth, distance to the coast and fishing effort intensity. The pattern for L was similar. C showed opposite trends when compared to N and L, as it decreased below 100 m depth, showing maximum values at intermediate fishing effort and Distance from the coast. G had a minimum around 100 m depth, being monotonically and directly related to Distance from the coast and fishing effort.

The reconstructed local food webs of the nine sites portrayed complex systems with a total of 101 taxa associated with more than 4 trophic levels ( Figure 3 ). It included several top predators such as C. conger, M. merlangus, S. megalops, L. piscatorius, L. caudatus, M. merluccius, D. dentex, T. marmorata and Z. faber, typical of the marine ecosystem of the South-central Mediterranean. Estimated TL was used to calibrate, parameterize, and validate food-web models and to calculate metrics and indicators of community structure (Belgrano, 2004). Trophic levels obtained from databases (Fishbase and Sealifebase) and the literature strongly correlated with those inferred from eDNA data using the analysis and visualisation of ecological communities (Yodzis and Innes, 1992) integrated into the cheddar package ( Figures 2-A1, A2 ). This indicates the robustness and accuracy of our approach for reconstructing food webs. Similarly, Agnetta et al. (2019) estimated the maximum trophic level of the SoS food webs at 4.5, which was reached by bluefin tuna Thunnus thynnus, with similar values for mesopelagic species, sharks, rays and demersal predators such as M. merluccius also placed in the upper trophic level ( Figure 3 ; Table 3 ). Furthermore, the trophic level values obtained from our approach strongly align with the study of Karachle and Stergiou (2017), which updated a previous compilation of the feeding habits and trophic level collected from 178 papers, regarding 204 Mediterranean fishes and invertebrates species, with values ranging from 2 to 4.50. Although further work is needed across a variety of different marine ecosystems under different levels of exploitation, our approach of comparing local food webs and trophic level appears promising for the characterisation of small-scale ecosystem responses to the interplay of environmental and exploitation drivers. The local food webs exhibited an association between increasing fishing effort intensity and changes in topological structures ( Figure 3 ; Table 1 ), highlighting the impacts of fishing effort on food web structures. The indicators of species richness, number of trophic links, direct connectance and generality are considered the most sensitive and reliable parameters to study spatial variations in food web structures (Martinez, 1992; Vermaat et al., 2009; Riede et al., 2010; Baiser et al., 2012). In our study, they were used as response variables in GAM models with fishing effort and abiotic variables serving as predictors. Results showed that species richness and number of trophic links seem to be the primary aspects of spatial variation in food web complexity in the SoS ( Figure 4 ; Table 4 ). Sites A1 and B12 (characterised by intermediate fishing effort intensity) showed higher values of both species’ richness and number of trophic links than sites with low/no fishing effort intensity ( Table 1 ). These findings revealed more complex and denser local food webs with highly connected species (higher number of trophic links) due to a high degree of generality (more pronounced trophic generalisms). Based on established theory (MacArthur, 1955; May, 1974), our results indicate that relatively high values of species richness and number of trophic links support high trophic complexity, which is considered fundamental in food webs stability (Pimm, 1980). Also, in line with classic ecological theory and with other studies (Cohen et al., 1985; Winemiller, 1989; Cohen et al., 1990; Olesen and Jordano, 2002; Dunne, 2006; Rossi et al., 2015), direct connectance, a proxy of ecological complexity, is inversely correlated which species richness ( Figure 4 ; Table 1 ), and this indicates that, when the species richness is higher, species are exploiting a growing fraction of the total ‘prey’ available (Warren, 1990). Moreover, the higher number of species detected in sites where fishing effort intensity is intermediate ( Table 1 ; Figure 3 ) can be explained by the dominant presence of mobile generalist fishes and invertebrates, which further mirror the higher values of trophic generalism with respect to the other fishing sites.

Previous studies have shown that the negative effects of fishing effort are not necessarily generalised at food web scale (Pauly et al., 1998; Kuparinen et al., 2016). The effect of fishing effort on the structure of local food webs and related indicators has been examined, using as reference the results published by Russo et al. (2019), which carried out an in-depth assessment of trawling areas in the Italian seas (including the whole SoS in the central Mediterranean). The values of fishing effort for the nine sampling sites of this study should thus be interpreted with reference to the values detected in the SoS (Russo et al., 2019). Accordingly, the values of the fishing effort for the nine sampling sites ranged from absent to intermediate, whereas none of the sampling sites considered in this study is characterised by truly high levels of fishing effort. Future works should overcome these limitations and investigate food webs and related indicators in highly-exploited locations. Hence, local food webs indicators we assessed showed significant spatial variations ( Table 1 ) and, in this sense, our analysis could be viewed as a quantitative exploration of how trophic re-organisation of the food webs could contribute to the understanding of marine exploited ecosystems. A major result of this study is that eDNA-based reconstructions of food web’s structure (i.e. species co-occurrence) were consistent over the nine sampling sites. However, the local food webs of different sites differed according to depth, distance from the coast and fishing effort intensity. The role of depth and distance from the coast (which are strongly related descriptors) was not surprising given that we are dealing with food webs of demersal communities. The effect of fishing effort intensity was less obvious and should be interpreted cautiously referring to the available literature. Indeed, fishing is expected to have adverse impacts on marine ecosystems (Pauly et al., 1998; Kuparinen et al., 2016), especially on food web structure and total biomass (Hinz et al., 2017). Specifically, fished areas are characterised by simplified food webs structures (Coll et al., 2008; Good et al., 2022). Lauria et al. (2020) using time series (1994–2016) of MEDITS trawl survey in the SoS, identified four main assemblages sustaining demersal fisheries: coastal (10-70m), outer shelf (70-250m), upper slope (250-450m) and medium slope (400-800m). These assemblages were stable throughout time and showed a clear spatial distribution, different species composition and biodiversity values driven by depth and, to a lesser extent, surface salinity. Our eDNA-based list of species is highly comparable with species found in two (coastal and outer shelfs) of the four assemblages described by Lauria et al. (2020). Moreover, conforming to these authors’ results, we observed that higher values of species richness can be detected in the coastal assemblage, despite being characterised by a relatively high fishing effort intensity. This apparent contradiction could be driven by the fact that the effect of fishing effort is masked by that of other features such as oceanographic features (fronts, upwellings and gyres) and habitat complexity of the coastal area of the SoS (Garofalo et al., 2007). In general, different studies have shown that increased perturbations tend to decrease diversity, and this has been observed in many demersal communities exposed to intense fishing effort (Rijnsdorp, 1996; Jennings and Kaiser, 1998; Gristina et al., 2006). Indeed, the action of bottom trawling is strongly related to the reduction or functional removal and/or loss of fishes and benthic invertebrates from the community, altering the ecological interactions such as predator-prey relationships (Hinz et al., 2017), competition strengths of ecosystems, and bentho-pelagic coupling (Agnetta et al., 2019). In contrast, other authors found that the relationships between human impact (quantified as cumulative removals) and species richness are complex, non-linear and ocean-basin specific (Trebilco et al., 2011). In particular, these authors documented that a higher fishing effort intensity is associated with higher species richness in the Indian and Pacific Oceans, and this overlap is particularly prominent in the central Pacific, and in the Indian Ocean around Sri Lanka (Trebilco et al., 2011). Through this study, we have provided the first indication of a significant positive relationship between rising fishing effort intensity and increasing species richness and number of trophic links in the Mediterranean Sea. If we consider that fishing effort levels corresponding to the sites of this study range from zero (no fishing effort) to intermediate, it is possible to explain the positive relationship between effort level and α diversity by referring to the Intermediate disturbance hypothesis (IDH - Dial and Roughgarden, 1998). Although it will be necessary to supplement and expand the observations of the present study to a wider range of environmental variables and fishing effort, it is certainly possible that the complex system of factors, including fishing effort, that determines α diversity in the SoS at the local level, allows for some level of biodiversity to be maintained even in areas moderately exploited by trawling.

Unfortunately, food web analysis requires a great deal of information, first and foremost the species composition of the communities analysed. Environmental DNA has largely been used for gathering indirect evidence about the presence of a broad spectrum of species, from single-celled organisms to whales (Jeunen et al., 2019). However, food webs reconstruction requires considerable knowledge about community organisms and their feeding habits. Leaving aside some largely acknowledged public databases such as FishBase and SeaLifeBase, accessing this knowledge from journals and reports is painstaking and time-consuming, and generating new knowledge incrementally using traditional sampling and observation is impractical for scaling up. Environmental DNA metabarcoding can boost biodiversity data acquisition to unprecedented levels revealing community compositions at high resolution, by making such efforts less reliant on increasingly scarce taxonomic expertise. Nevertheless, eDNA metabarcoding is not without drawbacks in the perspective of the ecological networks reconstruction. In accordance with basic ecological theories and concepts, including Adaptive Trophic Behaviour (ATB, Valdovinos et al., 2010) and Niche Theory, which can be associated with the concepts of ‘adaptive’ foraging (Kondoh, 2003) and foraging ‘optimisation’ (Pyke et al., 1977), consumers change their food sources over space and time, and thus their linkages in the food web, in response to resource availability, competition and presence of predators. In accordance with theories, it has been demonstrated that consumers will select a subset of prey among those available and an omnivore can feed at different trophic levels depending on prey abundance/productivity/season. The approach described in this paper is not able to capture these kinds of changes in food web structures and, more in general, eDNA metabarcoding should be groud-truthed with empirical studies about the actual feeding choices of the different species (Riccioni et al., 2018). However, we argue that: 1) the aggregation of taxa applied in this study can significantly reduce the distortion determined by the local and temporal ATB, since a given predator could “ignore” a given prey in some conditions, but it is unlikely that would ignore a full set of prey (e.g. copepods, malacostraca or cephalopods); 2) more important, none of the ecological indicators considered in this study is affected by the weights of the trophic links. Another potential issue is represented by the fact that the trophic position of omnivorous species depends on the relative importance of resources on different trophic levels in their diet but information on the strength of trophic linkages in specific locations at specific times cannot be obtained from the literature. To address this, we assumed that omnivorous predators consume prey from different trophic levels in the same proportions. Our study found that the values of TL estimated through food web analysis (TL_cheddar) are significantly coherent with the ones collected from the literature (TL_reference), which is per se a demonstration that the approach presented in this study can return a suitable estimation of the trophic positions. In summary, our methodological framework captured only changes in the food webs linked to changes in species presence/absence; it did not consider other sources of structural changes, such as changes in interaction strength. In this optic, the unique method to capture structural similarities or dissimilarities between food webs is to reconstruct heuristic food webs, based only on species co-occurrence. These latter are based on the assumption that the trophic links, unequivocally determined by the species feeding behaviour, are depicted by unweighted links in a binary network (binary matrix), where 0 represent no trophic relationship and 1 represent the trophic relationship of prey-predator, creating in this way a topological network. There is currently no complementary method to measure or otherwise generate the species abundance and their biomass data needed for FWA beyond quantified DNA sequence proxies. Furthermore, a main limitation of this approach results from the scarce knowledge of biological trophic links in marine communities, especially of less known environments (i.e. mesopelagic and bathypelagic zones). Despite all these considerations, to characterise the spatial variations of the SoS food webs we computed diverse topology properties because they cover diverse and ecologically relevant aspects of food web structures such as the vertical (trophic level), horizontal (generality) dimensions of food webs and their complexity (species richness, number of trophic links and direct connectance). Due to the complexity of our methodological framework, the interpretation of the results must be sequential and can rapidly become complicated if spatial components are entangled. By combining food web indicators with GAMs, we propose a framework to quantify and relate spatial changes in food web structure to changes in marine community composition. Our methodological approach can be readily used in areas where long-term monitoring of multiple trophic groups is carried out synchronously to support effective management strategies aimed at conserving the structure and functioning of ecological communities in times of environmental changes and species’ distributional shifts.