# A trophic latitudinal gradient revealed in anchovy and sardine from the Western Mediterranean Sea using a multi-proxy approach

Eneko Bachiller, Marta Albo-Puigserver, Joan Gimenez, Maria Grazia

Pennino, Neus Mari-Mena, Antonio Esteban, Elena Lloret-Lloret, Angelique

Jadaud, Belen Carro, Jose Maria Bellido, et al.

► To cite this version:

Eneko Bachiller, Marta Albo-Puigserver, Joan Gimenez, Maria Grazia Pennino, Neus Mari-Mena, et al. A trophic latitudinal gradient revealed in anchovy and sardine from the Western Mediterranean Sea using a multi-proxy approach. Scientific Reports, 2020, 10 (1), pp.17598. 10.1038/s41598-020-74602-y. http://dx.doi.org/10.1038/s41598-020-74602-y

HAL Id: hal-03405006

https://hal.umontpellier.fr/hal-03405006

Submitted on 25 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# A trophic latitudinal gradient revealed in anchovy and sardine from the Western Mediterranean Sea using a multi-proxy approach

Eneko Bachiller<fnref n="1" />,⁶✉, Marta Albo-Puigserver<fnref n="1" />, Joan Giménez<fnref n="1" />, Maria Grazia Pennino², Neus Marí-Mena³, Antonio Esteban⁴, Elena Lloret-Lloret<fnref n="1" />, Angelique Jadaud⁵, Belén Carro³, José María Bellido⁴ & Marta Coll<fnref n="1" />

This work combines state-of-the-art methods (DNA metabarcoding) with classic approaches (visual stomach content characterization and stable isotope analyses of nitrogen (δ¹⁵N) and carbon (δ¹³C)) to investigate the trophic ecology of anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus) at high taxonomic and spatial resolution in the Western Mediterranean Sea. Gut contents observed are in accordance with the dietary plasticity generally described for anchovy and sardine, suggesting a diet related to the opportunistic ingestion of available prey in a certain area and/or time. Genetic tools also showed modest inter-specific differences regarding ingested species. However, inter-specific and intra-specific differences in ingested prey frequencies and prey biomass reflected a latitudinal signal that could indicate a more effective predation on large prey like krill by anchovy versus sardine, as well as a generalized higher large prey ingestion by both species southwards. In fact, both species presented lower δ¹⁵N in the northernmost area. This latitudinal gradient indicates changes in the trophic ecology of anchovy and sardine that coincide with previously described better biological conditions for fish in the southern part of the study area as well as higher landings of both species in recent years.

European sardine (*Sardina pilchardus*) and European anchovy (*Engraulis encrasicolus*) are two of the most exploited small pelagic fish in the Mediterranean Sea, representing ca. 50% of the total Mediterranean fish landings<fnref n="1" />. In addition to their commercial importance, their high abundance and pivotal trophic position in the marine food web highlights their relevance for the Mediterranean ecosystem²⁻⁵.

During the last decade, declines in stock biomass of these species have concerned the fishing industry²,⁶. Historical observations over the past two decades in the Western Mediterranean Sea have shown decreases in stock biomass (and landings) for sardine in the north and fluctuations in the southern part⁷⁻¹⁰, whereas anchovy biomass has decreased or fluctuated in the north and increased in the south⁷,⁸,¹⁰. In addition, changes in body condition, growth, size at first maturity and disappearance of older ages have been observed for both species, also shows a latitudinal trend with higher incidence of change in the northern versus the southern part of the Western Mediterranean Sea⁷⁻⁹,¹¹.

In this sense, understanding ecological processes affecting fluctuations of such fisheries stocks has been a key issue in recent studies²,¹². In fact, these species are planktivorous consuming a wide range of prey during their life cycle¹³⁻¹⁹. Recent studies have suggested a relatively high niche overlap since 2010 and a reduction in prey diversity in the northern part of the Western Mediterranean²⁰, that might have caused fluctuations in stocks partially affected by inter-specific feeding interactions. However, most studies have been focused on the northwesternmost Mediterranean area, the Gulf of Lion²⁰, whereas the ecological status and trophic ecology of small pelagic fish in the whole Western Mediterranean Sea remains uncertain²¹.

<footnote n="1">¹Institut de Ciències del Mar (ICM-CSIC), Passeig Marítim de la Barceloneta, 37-49, 08003 Barcelona, Spain. ²Instituto Español de Oceanografía, Centro Oceanográfico de Vigo, Subida Radio Faro, 50, 36390 Vigo, Spain. ³AllGenetics and Biology SL, Edificio CICA, Campus de Elviña, 15008 A Coruña, Spain. ⁴Instituto Español de Oceanografía, Centro Oceanográfico de Murcia, Varadero 1 Apdo 22, 30740 San Pedro del Pinatar, Murcia, Spain. ⁵Marine Biodiversity, Exploitation and Conservation (MARBEC), Ifremer, University Montpellier, CNRS, IRD, Sète, France. ⁶Present address: AZTI, Sustainable Fisheries Management (Data), Basque Research and Technology Alliance (BRTA), Txatxarramendi uhartea z/g, 48395 Sukarrieta, Bizkaia (Basque Country), Spain. ⁷email: ebachiller@mail.com</footnote>
<table><caption>**Table 1.** Comparative methodological framework applied in this study to characterize the trophic ecology of anchovy (*Engraulis encrasicolus*) and sardine (*Sardina pilchardus*) in the Western Mediterranean Sea, using four different aspects (1–4), the corresponding comparable metrics and three sources of information (A–C).</caption><thead><tr><th>Variable</th><th>A. Gut content characterization</th><th>B. Stable Isotope Analysis</th><th>C. DNA metabarcoding</th></tr></thead><tbody><tr><td>1. Diet characterization</td><td>Feeding Intensity (SFD)<br/>Numerical frequency (%N)<br/>Prey biomass (%B)<br/>Frequency of occurrence (%FO)</td><td>δ<sup>13</sup>C and δ<sup>15</sup>N</td><td>%FO mOTUs/OTUs</td></tr><tr><td>2. Prey diversity/niche width</td><td>Species Richness<br/>Rarefaction curves (species richness)<br/>Shannon–Wiener diversity index</td><td>Ellipses width (SEA<sub>B</sub>)</td><td>Species Richness<br/>Rarefaction curves (species richness)</td></tr><tr><td>3. Diet similarity</td><td>Beta-diversity</td><td>–</td><td>Beta-diversity</td></tr><tr><td>4. Niche overlap</td><td>Pianka Niche Overlap</td><td>Ellipses area overlap (SEA<sub>B</sub>)</td><td>Pianka Niche Overlap</td></tr></tbody></table>

Recently, DNA metabarcoding has emerged as a useful technique to study the trophic ecology of several organisms<sup>22–24</sup>. However, these techniques are not so often combined with other methods to enhance their utility. In fact, methods like the high-throughput DNA sequencing (HTS) of diets<sup>25–27</sup> can be useful to parameterize food webs at enhanced taxonomic resolution, potentially improving the parametrization of ecosystem models, although they cannot quantify the amount of ingested prey<sup>28</sup>.

The direct quantification of ingested prey can be achieved by analysing gut contents under the microscope, which is time-consuming and requires high taxonomic expertise<sup>29</sup> or through mass-balanced stable isotope mixing models<sup>30</sup>, which require a proper isotopic reference collection of potential food sources and normally gives poor taxonomic resolution<sup>31,32</sup>. Furthermore, gut content analysis by visual inspection and DNA metabarcoding only provides a snapshot of what the fish has ingested in the last few hours<sup>29,33</sup>, which can lead to some bias if we are interested in the long term trophic ecology of the species. In addition, the gut content characterization technique also presents some limitations, as for example some gelatinous species and certain fish eggs that are difficult to detect due to their high vulnerability to the digestion processes<sup>34</sup>.

Complementary, stable isotope analysis provides information on the assimilated diet rather than the ingested one over a longer time period<sup>35</sup>. In the absence of potential prey stable isotope signatures (imperative to estimate the assimilated diet), they still provide useful information about trophic position, isotopic width and overlap, as a proxy of trophic niche width and overlap between species<sup>36,37</sup>. Thus, the combination of different techniques can provide new insights and show a general overview of the feeding ecology of organisms<sup>22</sup> (Table 1).

In this study, we characterize the trophic ecology of sardine and anchovy in the Western Mediterranean Sea, and explore potential geographical (i.e. latitudinal) variations in their diet, combining the three different methodologies mentioned above: microscope analysis and DNA metabarcoding of gut contents, and stable isotope analysis of muscle samples. In addition, we perform generalized additive models (GAMs) to test the influence of latitude, area, bathymetry, total length of fish on gut content indices and stable isotope variables for both species.

Our hypothesis is that due to an uneven stock status and population dynamics of the two species during the last decade, latitudinal changes in trophic ecology may partially explain apparent fluctuations in stocks. Previous studies have shown that the genetic structure of both species is similar in the Western Mediterranean Sea<sup>38</sup> except around south-western waters (Alboran Sea), where a transition area between the Mediterranean Sea and the Atlantic Ocean populations has been established<sup>9,38</sup>. Therefore, plausible causes of ecological change need to be investigated.

The present study represents the first approach that combines individual visual diet characterization of small pelagic fish with metabarcoding and stable isotope analyses in the Mediterranean marine ecosystem. This new information is essential to advance on our understanding of the ecological status and processes affecting fluctuations and change in population dynamics of anchovy and sardine in the study area. It also provides useful information about their feeding dynamics for future studies addressing potential trophic interactions between pelagic species sharing the ecosystem.

## Results

Diet characterization under the microscope, DNA metabarcoding and stable isotope analyses were combined in order to assess the trophic ecology of anchovy and sardine in terms of diet composition and diet similarity between both species, ontogenetic stages (juveniles vs adults) of each species and geographical areas by a latitudinal order (Table 1). The GSA07 corresponding with the Gulf of Lion was the northernmost region of the studied area, followed by GSA06-North in the Catalan Sea and Gulf of Valencia, and GSA06-South in the Gulf of Alicante (Fig. 1).

Specifically, the following aspects were investigated (see Methods; Table 1): (1) feeding intensity and prey composition in the diet, based on individual characterization of prey items in gut contents; (2) prey diversity in gut contents and niche width, assessed with species richness, Shannon–Wiener diversity index and isotopic Bayesian standard ellipse areas (SEA<sub>B</sub><sup>37</sup>); (3) similarity of diets, assessed by beta-diversity between anchovy and sardine in the three areas<sup>39,40</sup>; and (4) trophic niche and diet overlap, assessed by Pianka's<sup>41</sup> niche overlap and by isotopic SEA<sub>B</sub> overlap. In order to ease the interpretation of results, juvenile vs adult comparison figures are presented as Supplementary information.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sampling area in the Western Mediterranean Sea where (a) anchovy (Engraulis encrasicolus) and (b) sardine (Sardina pilchardus) were collected. The sampling area was divided in GSA07 (Gulf of Lion), GSA06-North (Catalan Sea and Gulf of Valencia) and GSA06-South (Gulf of Alicante). Locations where adults and juvenile individuals were collected are indicated with a red diamond for adults and a green cross for juveniles. All samples collected in the marked stations were used for gut content characterization (under the microscope and with DNA metabarcoding) and for isotope analysis (see Table 6 for details). Figure generated with QGIS v. 3.2.1-Bonn¹¹⁷ (https://qgis.org/en/site/).</figcaption></figure>

**Feeding intensity.** The best and most parsimonious generalized additive model (GAM) based on significant predictors, low AIC (Akaike Information Criterion) and high D (deviance explained) values included total length of fish, latitude and depth as explanatory variables of the stomach filling degree (SFD, a proxy of feeding intensity) variability, for both species (Table S01 online). The final models explained 66.3% and 76.2% of the deviances in anchovy and sardine, respectively (Table S01a online). Adults of both anchovy and sardine (i.e. larger fish) showed relatively higher SFD values, or higher feeding intensity, than juveniles (i.e. smaller fish) (Figure S01 online; Table S01b online). Significantly higher SFD values were observed for anchovy and sardine at lower latitudes (i.e. southern areas), and also in shallower depths in the case of sardine (Figure S01 online; Table S01b online).

**Diet composition.** Adults of both species ingested a broader range of prey sizes compared to juveniles, the latter demonstrating a diet composition generally dominated by copepods (Figure S02 online). Larger prey such as euphausiids and decapods were the main source of biomass in all species and areas, except for adult sardine and anchovy juveniles in GSA06-North, where copepods still represented more than half of the biomass input in the diet. Anchovy juveniles from GSA06-South also ingested large amounts of molluscs (gastropods and bivalve larvae) (Figure S02 online).

In terms of numerical frequencies, large krill was frequently ingested in GSA06-South, especially by adults of anchovy and sardine. This is in accordance with previous SFD analysis, since higher abundances of relatively larger prey in stomachs can partially explain higher SFD values in southern areas. In this sense, when GSA07 results (considering juvenile and adults together, Fig. 2a,b) were incorporated, although in general anchovy ingested relatively larger prey more frequently than sardine, some latitudinal differences were observed. While in northern areas (GSA07) anchovy preyed mainly on copepods but obtained most biomass from amphipods (i.e. 'Other Malacostraca', Fig. 2a,b), moving southwards, they fed more frequently on molluscs (mainly gastropod and bivalve larvae), euphausiids and decapods. Northern sardine based more than 50% of their diet on copepods (and some fish eggs contributing important amounts of biomass), and moving southwards sardine ingested more euphausiids and decapods, the highest numbers and biomass inputs coming from such large prey in GSA06-South (Fig. 2a,b).

When frequency of occurrence of prey is considered (%FO), results from the microscope analysis and DNA metabarcoding can be compared, obtaining new insights on diet composition analyses. This allows increasing the taxonomic resolution defined under the microscope to the highest (i.e. species) level, but is, at the same time,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Prey group composition per stomach determined under the microscope as percentage of (a) mean prey abundance, (b) mean prey biomass, and (c) occurrence, and (d) prey group occurrence determined with DNA metabarcoding, presented as means averaged across sampling sites within areas, for anchovy (*Engraulis encrasicolus*) and sardine (*Sardina pilchardus*).</figcaption></figure>

a trade-off between comparability of methods, as well as increasing the taxonomic resolution versus loosing relevant quantitative information (Fig. 2c,d, Tables S02-S07 online). Comparing both methodologies, our results showed that for both anchovy and sardine, copepods were the most abundant prey in numbers over the entire sampling area (Fig. 2a), whereas in terms of occurrence, the combination of methods showed that the relative importance of krill (i.e. euphausiids, decapods and amphipods defined as 'other Malacostraca') and other small prey like cladocerans or molluscs was higher than expected from the microscope analysis (Fig. 2c,d). Further, the relative importance of copepods in sardine diet from the northernmost area (GSA07) could be much lower than expected from microscope analysis (Fig. 2c), if we consider the presence of other prey groups only detected with DNA metabarcoding (e.g. decapods and other crustacean and molluscs, Fig. 2d).

Genetic methods (see detailed information in Tables S02-S04 and S06 online) allowed detecting 12 (10 demersal + 1 pelagic + 1 lantern-fish) fish species (as 'Actynopterygii cl., Fig. 2d), which could be presented as potential evidence of intraguild predation (i.e. predation of fish egg and/or larvae) by the two species in all areas.

DNA metabarcoding also indicated that the 'Others' group was composed of cnidarians (15 hydrozoan, 1 scyphozoan and 1 anthozoan groups), other molluscs (5 additional species), annelids (9 Polychaeta groups), echinoderms (7 groups), nemertines (1 group) and chaetognaths (1 group) (Fig. 3a). Results show that medusae were common in all areas and species, as well as annelids (with higher occurrences, at least based on percentages of detections, in GSA07) and echinoderms (especially important for sardine in GSA07) (Fig. 3a; Table S06 online).

Both anchovy and sardine ingested several taxa of diatoms, which were unidentifiable under the microscope but assessed with DNA metabarcoding (Fig. 3b). Results showed that the different species found in gut contents were ingested across the sampling area with no apparent trend. However, sardine from GSA07 seemed to have eaten different taxa from the others, as reflected by the 100% of occurrence of diatom groups that were not common (i.e. smaller occurrence % than the 70th percentile) in the other areas and species. It should also be noted that potentially harmful diatom species (i.e. within the Bacillariaceae and Chaetocerotaceae families; Table S04 online), as well as many rare species (Table S04 online), were observed in guts of anchovy and sardine (occurrence percentages of all identified taxonomic groups are presented in Table S07 online).

On the other hand, stable isotope analyses extended the diet traceability to a longer timescale, and allowed us to elucidate key differences for gut content analyses between areas for both species. Anchovy and sardine showed significant differences between areas in δ¹³C and δ¹⁵N (Table S08 online).

In the case of anchovy, δ¹⁵N changed from the lowest values in the north (i.e. GSA07) to the highest in the south (i.e. GSA06-South), in accordance with relatively larger prey ingested in the latter (observed in microscope analysis) (Figs. 4 and S03 online). Accordingly, the final selected GAM for δ¹⁵N of anchovy included the latitude as the only significant variable, which explained alone the 70.6% of the total variability, highlighting a negative relationship (Figure S04 online, Tables 2a and S09a online). Regarding δ¹³C of anchovy, the final GAM included not only the depth and latitude but also the total length of fish as significant variables. A negative linear
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Occurrence percentage of (a) phyla within 'Others' prey group and (b) diatoms, determined with DNA metabarcoding and presented as means averaged across sampling sites within areas for anchovy (*Engraulis encrasicolus*) and sardine (*Sardina pilchardus*) diet.</figcaption></figure>

relationship was found between latitude, depth and δ¹³C, while total length of fish showed a positive linear relationship. These results suggest that larger fish in shallow waters of the northern areas have higher δ¹³C values (Figure S05 online, Tables 2b and S09b online).

For sardine, a similar latitudinal trend was observed for δ¹⁵N (i.e. higher values in GSA06 than in GSA07; Figs. 4 and S03 online, Tables 2a and S09a online). In the selected GAM for δ¹⁵N, the final significant variables were the total length of fish and the area factor, which jointly explained 67.4% of the total variability (Figure S04 online, Tables 2a and S09a online). In particular, results highlight the lowest δ¹⁵N values in GSA07 with respect to GSA06-North (i.e. reference level, mean = -0.09, sd = 0.01) and GSA06-South (mean = -0.06, sd = 0.01) (Table 2a). Higher values of δ¹⁵N were found in larger (i.e. adult) fish (Figure S06 online, Table S08). The selected GAM model for δ¹³C showed bathymetry as the only significant variable, explaining 28% of the total δ¹³C variability and highlighting a decreasing trend from 80 m depth (Figure S05 online, Tables 2b and S09b online).

*Species richness and diversity.* As expected from the previous section, prey species richness as well as the Shannon–Wiener diversity index in diets were the highest in the southernmost area (GSA06-South, Table 3). A total of 144 prey groups were identified with DNA metabarcoding, 2.6 times more than the 55 groups identified under the microscope. DNA metabarcoding detected higher numbers of prey species (Fig. 5) in all taxonomic groups, e.g. copepods (40 vs 28 taxa), molluscs (10 vs 2), decapods (24 vs 4 taxa), fish (14 vs 3) or groups merged as 'Others' (42 vs 3 taxa).

In addition, DNA metabarcoding detected 41 species of diatoms, not observed under the microscope. This group was excluded from prey richness and diversity comparisons in order to make results comparable.

While differences within copepods were largely due to taxonomic identification (Table S02), most of fish, decapods and groups within 'Others' were only detected with DNA metabarcoding (Fig. 3a; Table S03 online). This resulted in higher prey richness indices than those obtained with microscope analysis, and small inter-specific differences regarding certain prey groups. For example, in GSA07, prey richness in anchovy and sardine were different depending on the method used to detect prey, possibly due to differences in the number of identified
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Corrected standard ellipses areas (SEA<sub>C</sub>; solid lines) in each area, for (a) anchovy and (b) sardine. Individual δ<sup>13</sup>C and δ<sup>15</sup>N values of GSA07 (red dots), GSA06-North (green dots) and GSA06-South (blue dots) are also graphed. Density plots showing Bayesian standard ellipses areas (SEA<sub>B</sub>) of (a) anchovy (*Engraulis encrasicolus*) and (b) sardine (*Sardina pilchardus*) are represented in the bottom-right of each graph as a measure of trophic niche width.</figcaption></figure>

copepod species (Table 3). Nevertheless, relative spatial differences remained similar between DNA metabarcoding and microscope analysis.

Accordingly, inter-specific differences in prey richness and diversity were generally low within each area, whereas a latitudinal decreasing trend in both prey richness and prey diversity was observed, which is partially related to the widening of the prey spectrum in the diet of both anchovy and sardine in the south (Table 3). In fact, the final selected GAMs showed a significant negative relationship between Shannon index and latitude for both anchovy and sardine, as well as a positive relationship with the length of fish (Figure S07 online, Tables 2c and S09c online). The interaction between the length of the fish and the latitude was significant for both species and highlighted that smaller fish sampled in the northernmost areas (GSA07) showed lower prey diversities than the larger fish sampled in the south (GSA06-South), with higher prey diversity in their diet composition (Table S12 online). Overall, the best GAMs explained 29.6% and 57.5% of the deviance for anchovy and sardine prey diversity, respectively (Table S09c online).

Regarding the isotopic niche width (Bayesian standard ellipse area, SEA<sub>B</sub><sup>36</sup>), contrary to what was observed with gut content analyses both anchovy and sardine showed a wider niche in GSA07. Specifically, in the case of anchovy the area that presented the widest SEA<sub>B</sub> was the GSA06-North followed by GSA07. The GSA06-south was the area where anchovy had the narrowest isotopic niche width. In the case of sardine, in the three areas the isotopic niche width was similar with higher values of SEA<sub>B</sub> for GSA07 and GSA06-South than GSA06-North, mainly due to a wider variability in δ<sup>13</sup>C values (Fig. 4; Figure S06 online).

*Similarity of diets, trophic niche and diet overlap.* In order to compare the similarity of diets between anchovy and sardine for each area, determined by microscope and DNA metabarcoding analyses, a beta-diversity approach was applied for anchovy and sardine with prey presence-absence data<sup>39</sup>. Assuming the same diet activity for both species, results obtained with both microscope analysis and DNA metabarcoding showed that
<table><caption>**Table 2.** Numerical summaries of the best (a) δ¹⁵N, (b) δ¹³C, and (c) Shannon diversity index (H') GAMs, obtained for anchovy and sardine (see model comparisons in Table S09 online). Statistics acronyms are: edf = degrees of freedom, Ref.edf = relative degrees of freedom, F = F statistic, Std.Error = Standard error.</caption><thead><tr><th colspan="5">A. Response variable: δ<sup>15</sup>N</th></tr><tr><th>Explanatory variables</th><th>edf</th><th>Ref.edf</th><th>F</th><th>p value</th></tr></thead><tbody><tr><td><strong>Anchovy</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Lat</td><td>2.59</td><td>2.87</td><td>65.05</td><td>&lt;2e-16</td></tr><tr><td><strong>Sardine</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>TL</td><td>2.86</td><td>2.98</td><td>8.43</td><td>9.22- e-05</td></tr><tr><td><strong>Parameter coefficients</strong></td><td><strong>Mean</strong></td><td><strong>Std.Error</strong></td><td><strong>t-value</strong></td><td><strong>p value</strong></td></tr><tr><td>Intercept (Area GSA06-North)</td><td>2.15</td><td>0.00</td><td>245.59</td><td>&lt;2e-16</td></tr><tr><td>Area (GSA07)</td><td>-0.09</td><td>0.01</td><td>-7.17</td><td>2.92e-10</td></tr><tr><td>Area (GSA06-South)</td><td>-0.06</td><td>0.01</td><td>-6.09</td><td>3.42e-10</td></tr><tr><th colspan="5">B. Response variable: δ<sup>13</sup>C</th></tr><tr><th>Explanatory variables</th><th>edf</th><th>Ref.edf</th><th>F</th><th>p value</th></tr><tr><td><strong>Anchovy</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Latitude</td><td>1.00</td><td>1.00</td><td>29.70</td><td>4.49e-07</td></tr><tr><td>Total Length (cm)</td><td>1.26</td><td>1.46</td><td>4.99</td><td>0.027</td></tr><tr><td>Depth (m)</td><td>1.00</td><td>1.00</td><td>4.75</td><td>0.032</td></tr><tr><td><strong>Sardine</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>Depth (m)</td><td>2.70</td><td>2.93</td><td>9.44</td><td>&lt;2e-16</td></tr><tr><th colspan="5">C. Response variable H'</th></tr><tr><th>Explanatory variables</th><th>edf</th><th>Ref.edf</th><th>F</th><th>p value</th></tr><tr><td><strong>Anchovy</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>TL</td><td>1.47</td><td>1.61</td><td>7.16</td><td>0.00222</td></tr><tr><td>Lat</td><td>2.00</td><td>2.00</td><td>6.71</td><td>0.00217</td></tr><tr><td>TL: Lat</td><td>7.87</td><td>8.89</td><td>4.69</td><td>8.85e-05</td></tr><tr><td><strong>Sardine</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>TL</td><td>2.92</td><td>2.99</td><td>3.12</td><td>3.07e-02</td></tr><tr><td>Lat</td><td>1.00</td><td>1.00</td><td>6.55</td><td>12.6e-03</td></tr><tr><td>TL:Lat</td><td>6.80</td><td>6.98</td><td>5.29</td><td>4.17e-05</td></tr></tbody></table>

anchovy and sardine shared the same niche in all areas equally (i.e. beta-diversity values were closer to 0 than to 1, see Methods). The highest differences between methods (i.e. microscope and DNA metabarcoding) were obtained in GSA06-South, probably due to the remarkable increase in numbers of copepods, decapods and mollusc species identified by DNA metabarcoding analyses. Accordingly, the inter-specific diet similarity was the highest in GSA06-South, where higher beta-diversity values were obtained compared to GSA07 (Table 3).

A high diet overlap (O) was found between anchovy and sardine, regardless of the methodology used for the analysis. For the whole sampling area, the O index was 0.98 and 0.96, considering diets based on microscope analysis and DNA metabarcoding, respectively. The O index was > 0.85 when considering each area separately (Table 4). When the overlap indices are compared between areas, the diet overlap was the highest between GSA06-North and GSA06-South for both species (O > 0.95, Table 4). The lowest diet overlap was observed for sardine between GSA07 and GSA06-North based on microscope analysis, but DNA metabarcoding also revealed a high degree of overlap in this case. In contrast, anchovy seemed to show lower niche overlap between these two areas regardless of the methodology used to analyse gut contents (Table 4).

When calculating trophic overlap with stable isotope data, a longer temporal period than the one covered with gut content analysis is integrated, obtaining information not only based on the prey consumed but also on the assimilated food. Thus, such overlap analysis of isotopic data (SEA<sub>B</sub>) between areas (i.e. GSA07, GSA06-North and GSA06-South) presented a general smaller overlap for each species (Fig. 4, Figure S06 online, Table S10 online). Some overlap was detected between GSA06-North and GSA06-South for anchovy (17.96% and 11.07% overlap) and between GSA07 and GSA06-South for sardine (30.16% and 28.17% of overlap). Regarding the overlap between species in each subarea, smaller overlap was detected except for GSA06-South, where a medium level overlap was apparent (27.53% and 54.80% overlap). SEA<sub>B</sub> overlap values and credible intervals are presented in Table S10 online).

Discussion

The gut contents characterized in this study is in accordance with the dietary plasticity generally described for anchovy and sardine<sup>13,17,42-46</sup>, suggesting that the diet composition of both pelagic species is generally related to the opportunistic ingestion of available prey in a certain area and/or at a certain point in time<sup>47,48</sup>. Genetic tools indicated small inter-specific differences regarding ingested species, which also supports the opportunistic
<table><caption>**Table 3.** Species richness and diversity in gut contents of anchovy (*E.enc: Engraulis encrasicolus*) and sardine (*S.pil: Sardina pilchardus*), averaged by prey groups (number of prey species), areas and methods. ‘Total’ denotes results for all the areas together. Shannon diversity index, calculated for gut contents analysed under the microscope, is presented as the mean ± standard error. Pairwise beta-diversity was calculated for anchovy and sardine according to Koleff et al.³⁹. *Diatoms were excluded in diversity calculations.</caption><thead><tr><th rowspan="3"></th><th colspan="8">Microscope analysis</th><th colspan="8">DNA metabarcoding</th></tr><tr><th colspan="2">GSA07</th><th colspan="2">GSA06-North</th><th colspan="2">GSA06-South</th><th colspan="2">Total</th><th colspan="2">GSA07</th><th colspan="2">GSA06-North</th><th colspan="2">GSA06-South</th><th colspan="2">Total</th></tr><tr><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th></tr></thead><tbody><tr><td>Calanoids</td><td>8</td><td>7</td><td>16</td><td>11</td><td>14</td><td>13</td><td>16</td><td>13</td><td>10</td><td>14</td><td>19</td><td>22</td><td>25</td><td>26</td><td>28</td><td>29</td></tr><tr><td>Cyclopoids</td><td>2</td><td>2</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>1</td><td>2</td><td>4</td><td>2</td><td>4</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Harpacti-coids</td><td>3</td><td>4</td><td>5</td><td>5</td><td>3</td><td>5</td><td>5</td><td>6</td><td>1</td><td>2</td><td>4</td><td>3</td><td>4</td><td>4</td><td>5</td><td>4</td></tr><tr><td>Euphausia-cea ord</td><td>3</td><td>0</td><td>5</td><td>3</td><td>5</td><td>6</td><td>5</td><td>6</td><td>0</td><td>0</td><td>3</td><td>0</td><td>2</td><td>2</td><td>3</td><td>2</td></tr><tr><td>Decapoda ord</td><td>1</td><td>0</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td><td>7</td><td>9</td><td>15</td><td>11</td><td>12</td><td>16</td><td>19</td><td>19</td></tr><tr><td>Other Malacos-traca</td><td>3</td><td>0</td><td>2</td><td>1</td><td>2</td><td>0</td><td>4</td><td>1</td><td>3</td><td>1</td><td>1</td><td>1</td><td>2</td><td>3</td><td>5</td><td>4</td></tr><tr><td>Crustacean remains</td><td>0</td><td>1</td><td>1</td><td>0</td><td>2</td><td>2</td><td>2</td><td>3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Mollusca ph</td><td>2</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>0</td><td>1</td><td>3</td><td>5</td><td>9</td><td>6</td><td>10</td><td>7</td></tr><tr><td>Cladocer-ans</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>4</td><td>3</td><td>4</td><td>2</td><td>3</td><td>3</td><td>4</td></tr><tr><td>Actinop-terygii cl</td><td>0</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>4</td><td>4</td><td>4</td><td>5</td><td>7</td><td>5</td><td>12</td><td>9</td></tr><tr><td>Others</td><td>0</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>6</td><td>11</td><td>22</td><td>30</td><td>26</td><td>23</td><td>31</td><td>35</td></tr><tr><td>*Diatoms</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>26</td><td>32</td><td>36</td><td>39</td><td>32</td><td>36</td><td>38</td><td>41</td></tr><tr><td>Species Richness</td><td>23</td><td>17</td><td>41</td><td>33</td><td>41</td><td>41</td><td>48</td><td>46</td><td>33</td><td>48</td><td>78</td><td>83</td><td>93</td><td>93</td><td>121</td><td>118</td></tr><tr><td>Shannon Diversity</td><td>1.27 ± 0.09</td><td>1.29 ± 0.08</td><td>1.57 ± 0.10</td><td>1.50 ± 0.11</td><td>1.72 ± 0.11</td><td>1.77 ± 0.07</td><td>1.52 ± 0.17</td><td>1.56 ± 0.17</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Beta-Diversity</td><td>0.35</td><td></td><td>0.24</td><td></td><td>0.17</td><td></td><td>0.14</td><td></td><td>0.44</td><td></td><td>0.28</td><td></td><td>0.27</td><td></td><td>0.20</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Sample-size-based rarefaction (solid line segments) and extrapolation (dotted line segments) sampling curves with 95% confidence intervals (shaded areas) for prey species richness obtained analyzing gut contents of anchovy and sardine under the microscope and with DNA metabarcoding. The solid dots/triangles represent the reference samples. Graphical representation was made according to Hsieh et al.¹¹⁸.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Table 4.** Pairwise contingency table presenting Pianka index of niche overlap between species (*E.enc*: anchovy; *S.pil*: sardine) and areas. Results in the top-right part of the table correspond to diet information obtained from visual microscope analysis, whereas the results in the bottom-left part of the table (numbers in italics) are based on DNA metabarcoding analyses. Numbers in bold represent species comparison of overlap indices within a certain area, and light numbers, the overlap indices between areas (i.e. considering the species separately).</figcaption></figure>

feeding and the ability of both species to switch between particulate and/or filter feeding strategies according to prey availability. However, inter-specific and intra-specific differences in ingested prey frequencies and the ingested prey biomass reflected a latitudinal signal that could illustrate a more effective predation on large prey like krill by anchovy than by sardine, e.g. especially in the *north* (hereafter referring to GSA07), as well as a generalized higher large prey ingestion by both species as they move *southwards* (hereafter referring to GSA06-South as the southernmost area) (Table 5).

Combining the stomach filling degree (SFD) approach with numerical prey frequency and prey biomass analyses, it seemed clear that predation on large prey is of great interest in terms of energy intake. Not only northern anchovy but also both species in the south get most of the biomass input from decapods, euphausiids and other malacostracans. It might be expected that such opportunistic active predation widens the prey size spectrum as fish grow, and that adult fish ingest higher numbers of larger prey⁴⁹ due to their energetic requirements (e.g. for reproduction⁵⁰) and a larger gape (i.e. mouth) width⁴². This high voraciousness in adults might also explain the higher SFD determined in adult anchovy and sardine. This does not necessarily mean that juveniles cannot get some success feeding on relatively larger prey. This is the case in GSA06-North, where sardine juveniles ingest important numbers of krill, in contrast with anchovy juveniles that get > 50% biomass from copepods; or in the south, where molluscs are often found in guts contents of anchovy juveniles. Such opportunistic predation by small (i.e. juvenile) anchovy and sardine has also been previously reported e.g. in the Bay of Biscay⁴²,⁴³,⁵¹.

As large prey are incorporated in the diet, prey richness and diversity also tend to increase, as reflected especially in the southern area for both species (also observed in the modelling results) (Table 5). This trend is in accordance with previous observations in the Bay of Biscay⁴³, but also contrasts with sardine studies a decade ago in North Aegean Sea that showed higher prey diversities related to higher phytoplankton ingestion by large fish¹⁶. The fish size range in the northernmost area was smaller than in the other areas (e.g. sardines were all juveniles, Table 6; sampled fish were smaller in GSA07, Figure S08 online). This may be caused by a lower growth and poorer feeding condition of northern Mediterranean fish observed in recent years by other authors²⁰,²¹,⁵². However, further research is required in order to determine a hypothetical size-dependent latitudinal feeding behaviour, while considering that sardine can effectively filter even the smallest prey sizes during their whole life-cycle¹³,¹⁶,¹⁷,⁴⁶,⁵³. In any case, in this study the major diversity compound is based on copepod ingestion, with several ingested species (33 and 43 groups based on microscope analysis and DNA metabarcoding, respectively), most likely determined by the available zooplankton in the area¹⁴,⁴⁹,⁵⁴. This results in low inter-specific or spatial differences and a relatively high degree of diet overlap (especially in the south). Then, when the number of ingested prey groups is low (e.g. by microscope analysis in northern anchovy and sardine), the diet dissimilarity (i.e. inter-specific difference) is relatively high. In contrast, incorporating DNA metabarcoding increases the number of detected species and therefore the beta-diversity³⁹, resulting on lower spatial inter-specific differences in such diet dissimilarities. Nevertheless, latitudinal differences remain clear regardless of the method applied, i.e. anchovy and sardine sharing the niche and showing higher diet similarity in the south (Table 5).

On the other hand, the isotopic niche overlap gives contrasting results compared to gut content data. The higher δ¹⁵N values of sardine in the northern and central part in comparison with anchovy values, explain the segregation of both species. Such higher values of δ¹⁵N are usually associated with the ingestion of prey with higher trophic position⁵⁵. In the north, inter-specific differences might be explained by the higher biomass obtained from fish eggs and/or larvae (i.e. Actinopterygii cl.) by sardine as well as by the high biomass input coming from decapods for anchovy. Instead, in the south the niche overlap between anchovy and sardine is higher, similar to the results from gut content analysis. Such an overlap between sardine and anchovy has also been previously reported¹⁴,²⁰,⁴⁷,⁵⁶. When comparing latitudinally, both species present lower δ¹⁵N in the northernmost area, which is in agreement with the observed prey composition²⁰, since a higher presence of larger
<table><caption>**Table 5.** Summary table showing comparative results regarding the latitudinal gradient hypothesis of the study, depending on species. (E.enc: Engraulis encrasicolus; S.pil: Sardina pilchardus), methods (GCC: gut content characterization; DNA-M: DNA metabarcoding; SIA: stable isotope analysis) and metrics (SFD: stomach filling degree; %ABD: percentage of prey number; %BIO: percentage of prey biomass; %FO: percentage of prey occurrence; S: prey richness; H: Shannon diversity index; β<sub>w</sub>: beta-diversity; O: Pianka overlap index; SEA<sub>C</sub> & SEA<sub>B</sub>: corrected and bayesian standard ellipses areas, respectively). “Copepods” refers to Calanoids, Cyclopoids and Harpacticoids; “Krill” refers to “Euphausiacea ord,” “Decapoda ord.” and “Other Malacostraca”; “Fish egg/larvae” refers to “Actinopterygii cl.”; remaining groups are presented as “Other prey”. “NLD” means no latitudinal difference (i.e. < 30% difference). Paired comparisons are made between the northernmost (GSA07, referred as N) vs. the southernmost (GSA06-South, referred as S) areas. “*” indicates comparison between methods. “>” and “<” indicate in which area the magnitude of the comparison is higher.</caption><thead><tr><th rowspan="2">Variable (metrics)</th><th>Method</th><th colspan="2">GCC</th><th colspan="2">SIA</th><th colspan="2">DNA-M</th></tr><tr><th>Species</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th><th>E.enc</th><th>S.pil</th></tr></thead><tbody><tr><td>Diet characterization</td><td>SFD</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>%ABD</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>copepods</td><td>N&gt;S</td><td>N&gt;S</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>krill</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>others</td><td>NLD</td><td>NLD</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>%BIO</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>copepods</td><td>NLD</td><td>N&gt;S</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>krill</td><td>NLD</td><td>S&gt;N</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>fish eggs/larvae</td><td>NLD</td><td>N&gt;S</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>other prey</td><td>NLD</td><td>N&gt;S</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>%FO</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>copepods</td><td>NLD</td><td>N&gt;S</td><td>-</td><td>-</td><td>NLD</td><td>NLD</td></tr><tr><td></td><td>krill</td><td>NLD</td><td>S&gt;N</td><td>-</td><td>-</td><td>NLD</td><td>NLD</td></tr><tr><td></td><td>fish eggs/larvae</td><td>NLD</td><td>NLD</td><td>-</td><td>-</td><td>NLD</td><td>N&gt;S</td></tr><tr><td></td><td>cnidarians</td><td>-</td><td>-</td><td>-</td><td>-</td><td>N&gt;S</td><td>NLD</td></tr><tr><td></td><td>other prey</td><td>NLD</td><td>NLD</td><td>-</td><td>-</td><td>NLD</td><td>NLD</td></tr><tr><td>S (values)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Prey diversity/niche width</td><td>[total]</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>S&gt;N</td><td>S&gt;N</td></tr><tr><td></td><td>copepods</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>S&gt;N</td><td>S&gt;N</td></tr><tr><td></td><td>krill</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>S&gt;N</td><td>S&gt;N</td></tr><tr><td></td><td>fish eggs/larvae</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>S&gt;N</td><td>NLD</td></tr><tr><td></td><td>other prey</td><td>S&gt;N</td><td>S&gt;N</td><td>-</td><td>-</td><td>S&gt;N</td><td>S&gt;N</td></tr><tr><td></td><td>diatoms</td><td>-</td><td>-</td><td>-</td><td>-</td><td>NLD</td><td>NLD</td></tr><tr><td></td><td>S (rarefaction curves)*</td><td>DNA&gt;GCC</td><td>-</td><td>DNA&gt;GCC</td><td>-</td><td>-</td><td>-</td></tr><tr><td></td><td>H and SEA<sub>B</sub> area</td><td>NLD</td><td>NLD</td><td>N&gt;S</td><td>NLD</td><td>-</td><td>-</td></tr><tr><td>Diet similarity</td><td>β<sub>w</sub></td><td>N&lt;S</td><td>-</td><td>N&lt;S</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Niche overlap</td><td>O and SEA<sub>B</sub> overlap</td><td>NLD</td><td>S&gt;N</td><td>NLD</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>**Table 6.** Number of stomach samples (N) and total length (TL) range (i.e. minimum-maximum length in cm) of anchovy. (E.enc: Engraulis encrasicolus) and sardine (S.pil: Sardina pilchardus). Results are presented by area, ontogenetic stage (juv: juveniles; ad: adults) and method applied for the analyses: gut content characterization under the microscope and DNA metabarcoding (i.e. the same fish were used for these two methods), and Stable Isotopes</caption><thead><tr><th></th><th></th><th>TL (cm)</th><th>Gut content characterization (N)</th><th>DNA metabarcoding (N)</th><th>TL (cm)</th><th>Stable Isotope Analysis (N)</th></tr></thead><tbody><tr><td rowspan="2">GSA07</td><td>E.enc</td><td>8.50-11.00</td><td>40 (juv)</td><td>10 (juv)</td><td>8.50-12.40</td><td>26 (juv) + 13 (ad)</td></tr><tr><td>S.pil</td><td>9.60-12.80</td><td>40 (juv)</td><td>10 (juv)</td><td>9.60-12.90</td><td>38 (juv)</td></tr><tr><td rowspan="2">GSA06-North</td><td>E.enc</td><td>8.90-14.80</td><td>15 (juv) + 20 (ad)</td><td>15 (juv) + 20 (ad)</td><td>8.90-14.50</td><td>14 (juv) + 14 (ad)</td></tr><tr><td>S.pil</td><td>7.80-16.90</td><td>15 (juv) + 15 (ad)</td><td>15 (juv) + 15 (ad)</td><td>7.80-16.90</td><td>12 (juv) + 12 (ad)</td></tr><tr><td rowspan="2">GSA06-South</td><td>E.enc</td><td>9.20-16.20</td><td>6 (juv) + 22 (ad)</td><td>6 (juv) + 22 (ad)</td><td>9.20-16.20</td><td>6 (juv) + 16 (ad)</td></tr><tr><td>S.pil</td><td>8.50-17.80</td><td>15 (juv) + 19 (ad)</td><td>15 (juv) + 15 (ad)</td><td>8.50-17.70</td><td>12 (juv) + 15 (ad)</td></tr><tr><td rowspan="2">All samples</td><td>E.enc</td><td>8.50-16.20</td><td>103</td><td>73</td><td>8.50-16.20</td><td>89</td></tr><tr><td>S.pil</td><td>7.80-17.80</td><td>104</td><td>70</td><td>7.80-17.70</td><td>89</td></tr></tbody></table>
prey in the south than in the north might also explain differences in the δ¹⁵N signal. Also, sardine and anchovy individuals were smaller in GSA07 than in the other areas, which could partially explain the lower δ¹⁵N values in the northernmost area. However, the isotopic niche width is similar across areas for both species, with slightly lower values in the south for anchovy and in the centre for sardine. The longer time-integrated period (e.g. weeks, months) reflected in stable isotopes in contrasting with the snapshot represented by gut content analyses could explain such differences in these results. Differences in the information provided from metrics on gut content and stable isotopes have also been highlighted in previous studies⁵⁷.

In addition, information obtained from certain prey species highlights differences both in terms of feeding efficiency and potential (i.e. energetic) interest in the diet. For instance, DNA metabarcoding allows determining several (previously underestimated) decapod species, such as *Solenocera membranacea*, indicator of muddy or sandy-muddy bottoms within the continental shelf and/or slope⁵⁸, which occurs in higher percentage of sardine gut contents in the north, as well as in anchovy gut contents in the south (Table S07 online). Similarly, DNA remains of other species like the annelid *Magelona* sp. and starfish *Paracentrotus lividus* were also identified in both fish species all over the sampling area, suggesting that the feeding activity of anchovy and sardine might also affect to (most likely) pelagic early-life stages of demersal species. Moreover, when samples were collected during the daytime (see Methods), both anchovy and sardine could also be feeding in deeper areas¹⁷,⁵⁹, which highlights the importance of these species not only for pelagic but also for the demersal-pelagic energy transfer⁶⁰. On the other hand, DNA metabarcoding reveals an important presence of fish eggs and/or larvae in (mainly juvenile fish) gut contents, especially of demersal species, such as *Trisopterus capelanus* and *Spicara maena*, but also of pelagic species like *Sprattus sprattus* (mainly sardine predation). However, just with occurrence data it is difficult to determine whether such predation might cause potential (negative) effects in survival indices of early life stages of potential prey and competitors⁶¹⁻⁶³.

In any case, our study showed that both anchovy and sardine took advantage of whatever prey species was available for them; and fish predation, as well as the predation on relatively large prey groups, may occur especially during late hours of the day and/or night, when such prey are more vulnerable to active predation¹⁸. The same may apply for the predation on jellyfish, which have not been detected in the diet of anchovy and sardine before this study. Such jelly organisms might be considered as especially vulnerable during the night, and were mainly identified by DNA metabarcoding probably because they were mostly digested by the time of fish sample collection (i.e. daytime). Determining jelly organisms ingestion by fish is of special interest, since such expanding organisms could be considered indicators of global warming and changes in the pelagic ecosystem⁶⁴,⁶⁵. However, it is not clear whether hydrozoans have not been ingested in the past, or the lack of such information is related to methodological issues. Recently, other jelly organisms like salps have been detected as part of the diet in *Sardinella aurita* within the same area⁶⁶, which is in accordance with the observed ingestion of siphonophores such as *Nanomia bijuga* and *Muggiaea atlantica*, or hydrozoans like *Clytia hemisphaerica*, and medusae *Lizzia blondina* and *Aglaura hemistoma*, all of them often detected as DNA remains of gut contents in anchovy and especially in sardine in this study (Table S07 online). Further research is needed to determine potential consequences of the increase in numbers of cnidarians that might be expected in the diet of anchovy and sardine in a near future due to environmental changes.

Genetics also determined that both anchovy and sardine ingest diatoms. It is generally known that phytoplankton is more important in the diet of sardine due to their more effective filtering apparatus (i.e. gill-raker size⁴⁶,⁶⁷), and the occurrence and diversity of phytoplankton in gut contents is higher in sardine. According to our study, the ingested diatom species by anchovy and sardine are most likely determined by the phytoplankton diversity in the feeding area during the sampling period; such diatoms might have been directly ingested through filter feeding, or even as part of the zooplankton ingested by fish prey, such as copepods and other phytophagous crustaceans or molluscs. That said, for instance, *Chaetoceros* spp. and *Thalassiosira* spp., known as common spring bloom species in the Mediterranean⁶⁸⁻⁷⁰, are found in guts from both anchovy and sardine. Other diatoms such as *Minidiscus trioculatus*, reported as 'extremely rare' in previous work²¹ are commonly found in samples of both anchovy and sardine in the present study (Tables S05 and S07 online), indicating a potential expansion of such species distribution in the Mediterranean during recent years. The random phytoplankton ingestion is also confirmed with the observed harmful diatom species⁷⁰, such as *Pseudo-nitzschia* spp., *Chaetoceros socialis* and *Cerataolina pelagica* (Tables S05 and S08 online), which might not cause further problems in anchovy or sardine population, unless blooms of this kind of algae occur.

It should be noted that DNA metabarcoding have some limitations when is used in diet-related studies, such as DNA degradation from the time of feeding until the sample collection, DNA extraction efficiency, selection of molecular markers, PCR inhibition or uncompleted reference databases⁷²⁻⁷⁴. In addition, DNA metabarcoding without microscope analysis does not determine the trophic position of each detected prey species in the food web. In the same way, further development is needed in this method in order to quantify the detected prey, which would increase the sampling coverage and the diet characterization in the laboratory more cost-effective. However, at this stage our DNA metabarcoding approach has increased the detection success of prey species and, complementary to the results obtained from the other two methods, provides relevant information of anchovy and sardine diet composition.

In summary, this study reports a latitudinal gradient in the diet composition of both anchovy and sardine in the Western Mediterranean Sea, which may reflect a widening of the prey spectrum for both anchovy and sardine as they move southwards (Table 5). All applied methods show that the northern area is the most different one, where both species showed a higher inter-specific difference in the diet probably related to poorer feeding conditions (i.e. less numbers or relatively large prey²⁰) that lead to higher phytoplankton and small copepod ingestion especially by sardine, and a more effective opportunistic (active) predation on large prey by anchovy³, 13, 16, 17, 20, 21, 53, 54. The largest ingestion of biomass in the diet of both species came from large prey, such as fish eggs and cnidarians in the north for sardine and anchovy, respectively, or krill, ingested by all fish all over
the sampling area but with higher apparent success in the south. In fact, according to this study, in the south, inter-specific feeding differences were smaller, most likely because anchovy and sardine took advantage of more abundant and diverse prey groups, including a wider range of large prey. The generalist feeding behaviour of anchovy and sardine is known to result in intra-specific individual variation in the diet composition⁴²,⁴³, as observed also in the present study.

Overall, every method used in this study provided specific diet information of anchovy and sardine, and methodological differences were mainly reflected in the obtained absolute numbers (e.g. higher diversity and richness values obtained with DNA metabarcoding, Fig. 5)²². Only when all of the methods were combined, a global insight of the feeding ecology could be obtained. This is the very first time these three different methods are combined to investigate feeding ecology of pelagic fish in a broad spatial scale. This is indeed the greatest contribution of this work as, for instance, when considering different time scales, the whole prey size range and other potential interactions could be often underestimated when using methods separately. It should be noted also that despite the use of different methods, the latitudinal trends were clearly detected by all methods, illustrating changes of the trophic ecology of anchovy and sardine in the Western Mediterranean Sea. This latitudinal gradient matches with the ecological gradient previously described and the better conditions of both species in the southern area in recent years⁹. The evidenced inter-specific latitudinal difference is fundamental for further multidisciplinary approaches that could integrate different aspects such as fish distribution²,⁷⁵, plankton availability⁷⁶, fish growth²⁰,⁵² or the potential ecological effects of microplastic ingestion by fish⁷⁷,⁷⁸.

## Methods

**Sample collection.** Fish samples of anchovy and sardine were collected in the Western Mediterranean area (Fig. 1) in Geographical Sub-Areas 6 and 7 (GSAs) defined by the General Fisheries Commission for the Mediterranean (GFCM¹), during the MEDITS 2018 survey (Mediterranean International bottom Trawl Survey). MEDITS was considered an ideal survey for obtaining standardized information (i.e. abundance index by swept area method), as researchers use the same sampling protocol throughout the Mediterranean Sea⁷⁹. Therefore, the sampling procedures were standardized according to a common protocol over GSAs and years. The standard fishing device was a bottom trawl GO73 with 20 mm cod-end mesh size net⁸⁰. The average vertical opening of the gear was 2 m and its wing-span was 18 m. All the tows were performed during daylight hours. The samples of the study were collected during the months of May and June, 2018.

**Gut content characterization under the microscope.** Gut contents of 103 anchovy and 104 sardine samples (Table 6) were analysed individually, with no subsampling, under a NIKON SMZ1270 stereomicroscope with 20–80× amplification. The laboratory practice followed the standard procedures to avoid contamination during sample processing. Accordingly, microscope analysis was conducted in a ‘clean room’ and with an air extractor placed 20–30 cm above the petri plate with the gut samples. To avoid contamination between samples, glassware, bench, microscope slide and dissection equipment (i.e., stainless-steel scissors, scalpel and lancet) were rinsed with 96% ethanol prior to each gut content analysis⁸¹.

In order to exclude bias caused by different rates of digestion and cod-end feeding⁸² only material contained in the stomachs was considered, whereas the contents of the intestine and oesophagus were discarded. During processing, stomach contents were carefully taken apart and all identifiable prey counted and specified to the lowest possible taxonomic group, not including broken parts of appendixes when quantifying, and categorized into 59 groups. Parasitic organisms found in stomachs (e.g. Trematoda and Nematoda larvae) did not show any relationship either with total prey abundance in gut contents or with stomach fullness, so they were excluded from the diet analysis. After the microscope analysis, stomach contents were preserved on 96% ethanol for later DNA metabarcoding analysis.

The feeding intensity was assessed calculating the stomach filling degree (SFD). This parameter is a useful qualitative metric that allows determining if feeding intensity (or efficiency) is relatively higher in a certain area and/or time⁴³,⁸³. In order to exclude the effect of fish size, the response variable SFD was defined as the sum of the weights of all the prey in a stomach (mg) divided by the total length of the fish (mm).

The diet composition was first explored using numerical and weight percentages of prey groups relative to total prey consumption. To determine the weight of each prey group, length–weight conversion equations were used based on literature. Average total length of prey species (mm) was obtained from Bachiller & Irigoien⁴², where direct length measurements were made for the first 50 prey items in each stomach. Since most of prey species observed in the present study were also observed in the Bay of Biscay, conversion of prey counts into biomass was made to the detailed species level. For missing prey species or groups (i.e. observed only in the Mediterranean), biomass of the same genera or the corresponding upper taxonomic level was assigned (Table S03 online). To exclude the effect of the sample size (number of fish per station) on the identified prey abundance, the biomass of each prey group was weighted by the number of fish per predator species and area.

Diet composition was presented as percentages of the total prey consumption for the three areas when considering all sampled fish, and for GSA06-North and GSA06-South also by ontogenetic stages (i.e. adults vs. juveniles); fish < 11 cm and < 13 cm were considered as juveniles for anchovy⁸⁴,⁸⁵ and sardine⁸⁶, respectively.

To ease later interpretation of the figures, prey groups were categorized into 11 groups (see Table S04 online): Calanoids, Cyclopoids, Harpacticoids, Euphausiacea ord., Decapoda ord., Other Malacostraca, Crustacean remains (i.e. undefined taxa and/or broken parts of organisms within Crustacea subph.), Mollusca ph., cladocerans, Actinopterygii cl., and Others (including the rest of the groups with a frequency in number < 5% of the total prey consumption observed under the microscope). This ‘Others’ group was then broken down in detail, based on results from the DNA metabarcoding and presenting 38 taxonomic groups, merged into 8 phyla (Table S04). Finally, the diatom (Bacillariophyta ph.) content in gut contents of anchovy and sardine was also assessed by the
DNA metabarcoding (see procedures below), detecting 41 different algae taxa, merged into 11 groups (i.e. based on family taxonomic level, and presenting algae with < 7% of the total algae occurrence as 'Diatom remains' for graphical representation; Table S05 online).

Prey species richness was defined as the number of different prey groups found in gut contents. The diet diversity was expressed by the Shannon-Wiener diversity index (H'), calculated from prey abundance composition based on the extended prey species list obtained from microscope analysis (i.e. 59 groups).

Beta-diversity (β<sub>w</sub>) for anchovy and for sardine was calculated using prey presence-absence data (considering each area separately), according to the following equation 40:

$$ \beta_w = \frac{a + b + c}{(2a + b + c)/2} - 1 \qquad (1) $$

where component a comprises the total number of species that occur in both anchovy and sardine; component b comprises the total number of species that occur in anchovy but not in sardine; and component c comprises the total number of species that occur in sardine but not in anchovy. Regarding this pairwise comparison and assuming that the diet composition is similar for both species, a minimum value of 0 beta-diversity (similarity measures) would mean that both species share the niche equally, whereas a maximum value of 1 (dissimilarity measures) might be expected when one community dominates the ecosystem<sup>39</sup>.

The overlap in resource use between species (i.e. anchovy and sardine) and between areas was assessed using Pianka's<sup>41</sup> index of niche overlap:

$$ O = \frac{\sum p_{i,j}p_{i,k}}{\sqrt{[(\sum p_{i,j}^2)(\sum p_{i,k}^2)]}} \qquad (2) $$

where O is the overlap index between the two species j and k expressed as a value between 0 and 1, where 0 means no overlap and 1 complete overlap in diets. P<sub>i,j</sub> and P<sub>i,k</sub> are the proportions of presence of the i<sup>th</sup> prey group in the diets of species j and k, respectively. To assess the overlap index between areas (considering the species separately), 'species' were replaced with 'areas' in the same equation. For the diet overlap comparisons, the gut contents were categorized into the 11 prey groups mentioned before and used for graphical representations. To test for significance, the presence-proportion of a given prey group in a given diet was randomized according to the Randomization Algorithm (RA2) defined by Lawlor<sup>87</sup> and iterated 1000 times for each comparison of diet overlap. Lawlor<sup>87</sup> described four randomization algorithms (RA1-RA4) for niche overlap, in which the zero states (the empty prey groups) and the niche breadth (the degree of utilization of a prey group) can be either relaxed or retained. Under RA2, the zero states are retained (i.e. empty prey groups from the stomach samples remain empty in the simulations), while niche breadth is relaxed (i.e. the proportion in the diet of each non-empty prey group is replaced by a uniform value between 0 and 1). As in the case of the Bay of Biscay<sup>43</sup>, RA2 was considered to give the most realistic reflection of the Mediterranean Sea pelagic system because some of the prey groups would be unavailable to fish in certain areas, due to the patchy distribution of the plankton prey<sup>88,89</sup>, whereas none of the fish species were assumed to have constraints on the utilization of the prey groups that were actually present.

**DNA metabarcoding in gut contents.** DNA was extracted from 143 gut content samples of anchovy and sardine analysed previously under the microscope (Table 6), using the NZY Tissue gDNA Isolation kit as per manufacturers protocol (NZYTech, Lisbon, Portugal). Prior to DNA extraction, vials were shaken by hand to homogenise the gut contents of each individual sampled fish. An extraction blank was included in every DNA extraction round and treated as a regular sample to check for cross-contamination. DNA was resuspended in a final volume of 100 µL.

*Zooplankton and diatom characterization in diet.* For library preparation of zooplankton, a 313 bp-fragment of the Cytochrome Oxidase subunit I (COI) region was amplified using the primers COIntF (5' GGWACWGGW TGAACWGTWTAYCCYCC 3') and COI dgHCO2198 (5' TAAACTTCAGGGTGACCAAARAAYCA 3')<sup>90</sup>, to which the Illumina sequencing primer sequences were attached to their 5' ends. In order to limit the amplification of predator DNA, specific blocking primers were designed following Leray et al.<sup>91</sup> using Geneious 11.1.5 (www.geneious.com), based on COI sequences of *Sardina* sp. and *Engraulis* sp. The blocking primers *dgHCO2198_engraulis* (with 5'-3' sequence AAGAATCAGAATAGGTGTTGATAAAGAATC-C3) and *dgHCO2198_sardina* (with 5'-3' sequence AAGAATCAGAATAGGTGCTGATACAGAATG-C3) were used during PCR in order to prevent amplification from anchovy and sardine, respectively. A C3 CPG spacer was added to the 3' end of each blocking primer to prevent elongation.

Instead, for library preparation of diatoms, a fragment of around 312 bp of the rbcL chloroplast gene was amplified using the primers Diat_rbcL 708F_2 (5' AGG TGA AGT TAA AGG TTC WTA YTT AAA 3') and R3_1 (5' CCT TCT AAT TTA CCW ACW ACT G 3')<sup>92</sup>, to which the Illumina sequencing primer sequences were attached to their 5' ends.

The libraries were constructed by PCR amplification of a region of the COI gene for zooplankton and a region of the rbcL gene for diatoms. PCRs were prepared in a volume of 25 µL, containing 1 µL of template DNA, 0.5 µM of the primers, 10 µM of the blocking primers, 6.25 µL of Supreme NZYTaq 2× Green Master Mix (NZYTech, Lisbon, Portugal), and ultrapure water up to 25 µL. The thermal cycling profile included an initial denaturation at 95 °C for 5 min, followed by 35 cycles of 95 °C for 30 s, 53 °C for 30 s, 72 °C for 30 s, and a final extension step at 72 °C for 10 min. The same protocol but using only 5 cycles and 60 °C as the annealing temperature was used in a second PCR round to add index sequences for multiplexing different libraries in the same sequencing pool
(see Fig. 1 in Vierna et al.⁹³). The product sizes of both PCRs were checked by electrophoresis and visualised on 2% agarose gels. Then, COI and rbcL libraries were purified using the Mag-Bind RXNPure Plus magnetic beads (Omega Biotek, Norcross, USA), quantified using the Qubit dsDNA HS Assay (Thermo Fisher Scientific, Waltham, USA), and pooled equimolarly. Finally, the pool was sequenced in a MiSeq PE300 run (Illumina, San Diego, USA) at the facilities of Parque Científico de Madrid Foundation (Madrid, Spain).

The quality of the Illumina paired-end raw files was checked using the software FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc). Merging of paired-end reads was performed with FLASH2⁹⁴. The mismatch resolution in the overlapping region (<30 bp -base pairs-) was accomplished by keeping the base with the higher quality score. The CUTADAPT software 1.3⁹⁵ was used to remove sequences that did not contain the PCR primers (allowing up to 2 mismatches) and sequences that ended up between 293 and 333 bp for the COI region and between 250 and 375 bp for the rbcL region. Sequences were quality-filtered (Phred score > 20) and labelled using the script multiple split_libraries.py implemented in Qiime v1.9.1⁹⁶.

The resulting FASTA file was processed using the VSEARCH bioinformatic tool⁹⁷. Sequences were dereplicated (-derep full length), clustered with the SWARM 2.0 algorithm⁹⁸ with a *d* value of 13, and sorted (-sortbysize). Then, chimeras were de novo detected and removed using the UCHIME algorithm⁹⁹ implemented in VSEARCH. Finally, sequences were assigned to Operational Taxonomic Units (OTUs) (-usearch global).

In the case of the zooplankton taxonomic identification, a custom reference database was created using an in-house developed script to process the information from the BOLD Public Data Portal (accessed on July 2019). All the COI sequences for the following metazoan groups were retrieved: Actinopterygii, Appendicularia, Asciidiacea, Branchiopoda, Ciliphora, Chaetognatha, Cnidaria, Hexanauplia, Malacostraca, Mollusca, Polychaeta and Thaliacea. The taxonomic identification was performed by querying the clustered centroids against our BOLD reference database using the script *assign_taxonomy.py* implemented in Qiime and the UCLUST algorithm¹⁰⁰, with a sequence similarity threshold of 95%.

In the case of the diatoms taxonomic identification, it was performed by querying the sequences against the R-Sys reference database⁹². To do so, we used the naïve Bayesian method implemented in RDP¹⁰¹, with a confidence threshold of 80%.

OTUs occurring at a frequency below 0.005% in the whole dataset were removed¹⁰². In DNA metabarcoding studies it has been observed that a low percentage of the reads of a library can be assigned to another library. This phenomenon, referred to as *mistagging* is the result of the misassignment of the indices during library preparation, sequencing, and/or demultiplexing steps¹⁰³,¹⁰⁴. In order to correct for this phenomenon, the OTUs occurring at a frequency below 0.01% in each sample were removed.

OTUs that were not assigned to diatom taxa were removed from the final OTU table. In the case of zooplankton, OTUs without taxonomic assignment were queried directly against BOLD, using a similarity match of ≥95%. The OTUs that got a taxonomic assignment were included in the final OTU table. OTUs assigned to *Sardina* sp. and *Engraulis* sp. were removed from samples belonging to *Sardina pilchardus* and *Engraulis encrasicolus*, respectively. In the case of zooplankton, samples that ended up with <20 sequences were also removed from the final OTU table.

Samples in the OTU table were sorted by sampling site, ontogenetic stage, and area. Then, the resulting OTU tables were converted into a Biological Observation Matrix files (.biom) for the posterior analysis of the data. DNA metabarcoding analyses were carried out by AllGenetics & Biology S.L. (www.allgenetics.eu).

Presence-absence information of prey (zooplankton) OTUs was used for diet characterization as well as to calculate the species richness, the beta-diversity for anchovy and for sardine (considering each area separately) and niche overlap, according to the same equations applied to gut content characterization data (see previous section).

The diatoms analysis was considered for diet characterization as well as for prey richness approach. In order to ease the interpretation of graphical representation, diatom groups within the 70th percentile of occurrence frequencies (i.e. <6.25) were merged as 'Other diatom groups'. Diatoms without specific taxonomic information were presented as 'Diatom remains'.

All the raw reads of prey species obtained from gut contents were deposited in the BioProject database (GenBank - NCBI¹⁰⁵) as Bioproject PRJNA653773.

**Stable isotope analyses.** Stable isotope analyses of δ¹³C and δ¹⁵N were performed in muscle tissue. A small portion of the dorsal muscle without skin was extracted from each fish sample. Then, tissue samples were oven-dried at 60°C for 72 h. Dried samples were pulverized and 0.80 – 0.85 mg of muscle powder per fish was packed into tin capsules.

Isotopic analyses were performed at the Laboratory of Stable Isotopes of University of A Coruña, Galicia, Spain (Servicio de Analisis Instrumental (SAI)) through an elemental analyser (Carlo Erba CHNSO 1108) coupled to an isotopic ratio mass spectrometer (Finnigan Matt Delta Plus). The isotopic values are reported as δ¹³C (‰) and δ¹⁵N (‰) relative to the Vienna Pee Dee Belemnite or atmospheric nitrogen respectively¹⁰⁶. USGS40 and L-alanine from the International Atomic Energy Agency were used, as well as internal acetanilide standards. The accuracy (± SE) of the standards replicates and samples for the two isotopes is < 0.1 and < 0.3% respectively. When the C:N ratio was greater than 3.5, it indicates that lipids are present in the sample¹⁰⁷. Therefore, in these cases a correction was applied to the values of δ¹³C following the methodology of Post et al.¹⁰⁷.

To provide insight into species' trophic niche widths and to estimate the degree of isotopic niche overlap between the species, standard ellipse areas corrected for the sample size (SEA<sub>C</sub>; i.e. the area containing 40% of the data) as well as the Bayesian standard ellipse areas (SEA<sub>B</sub>⁹⁷) were calculated for each species, ontogenetic stage (i.e. juveniles vs adults) and area (i.e. GSA07, GSA06-North, GSA06-South). This metric represents a measure of the core isotopic niche and higher values of SEA<sub>C</sub> and SEA<sub>B</sub> represent a broader trophic niche width⁵⁵. Overlap
between SEA<sub>B</sub> were calculated as percentage of the area of species A overlapped with B and vice versa. SEA<sub>B</sub> was computed using 10,000 posterior draws<sup>37</sup>.

**Generalized additive models.** Generalized Additive Models (GAMs<sup>108</sup>) were implemented to test the influences of latitude, area (GSA07, GSA06-North, GSA06-South), depth (m), and total fish length (cm) on the SFD, Shannon diversity index, δ<sup>13</sup>C and δ<sup>15</sup>N of anchovy and sardine. Explanatory variables were tested for collinearity, correlation, outliers, and missing data, before using them in GAMs, following the procedure of Zuur et al.<sup>109</sup>. In particular, correlation among variables was checked by performing a Pearson's correlation test. Collinearity was tested by computing the generalized variance-inflation factors (GVIF), which are the corrected VIF values by the number of degrees of freedom of a predictor variable<sup>110</sup>. In all cases variables had a correlation lower than 0.70 and a GVIF lower than 3 and consequently were used together in models (Figure S08).

GAMs are often used for their ability to deal with non-linear and non-monotonic relationships between the response variable and the explanatory variables<sup>108,111</sup>. Separated GAMs for each one of the response variable and species were performed. When response variable was not normally distributed a logarithmic transformation was applied; when the normality was not achieved after transformation, a Gamma distribution was used. Accordingly, for both, anchovy and sardine, the SFD was not normally distributed (Shapiro-Wilk normality tests; W = 0.90, *p* value < 0.001 and W = 0.83, *p* value < 0.001, respectively) and therefore a Gamma distribution with a log link was used. On the contrary the Shannon index was normally distributed (Shapiro-Wilk normality tests; W = 0.95, *p* value < 0.001 and W = 0.97, *p* value = 0.06 for anchovy and sardine, respectively) and a Gaussian distribution with the identity link was implemented. For δ<sup>13</sup>C a Gaussian distribution with a log-link was applied (Shapiro-Wilk normality tests; W = 0.98, *p* value = 0.40 and W = 0.98, *p* value = 0.10 for sardine and anchovy, respectively), as well as for δ<sup>15</sup>N (Shapiro-Wilk normality tests; W = 0.96, *p* value < 0.001 and W = 0.95, *p* value < 0.001 for sardine and anchovy, respectively).

Each GAM was fitted using thin plate regression splines for non-linear covariate and restricting the number of knots at 4 to avoid additional over-fitting. GAMs were performed with all possible combinations of terms. Variables were selected with forward and backward stepwise procedures based on three different criteria such as Akaike Information Criterion (AIC) and deviance explained. The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, high D (deviance explained) values, and significant predictors.

**Software.** R software v.3.6.1<sup>112</sup> was used for all analyses and graphical representations (except Fig. 1, see corresponding legend). For diet composition figures the package 'ggplot2' v.3.2.1<sup>113</sup> was used. Diet overlap analyses in diet characterization data were performed using 'EcoSimR' package<sup>114</sup>. GAMs were performed using 'mgcv' package<sup>115</sup>, testing autocorrelation among variables and collinearity using 'corrplot' package<sup>116</sup> and 'corvif' function, respectively. Isotopic SEAs and their overlap were calculated using SIBER package<sup>37</sup>. OTU table-based matrices were directly imported into R<sup>112</sup>.

**Ethical statement.** Our study did not involve any endangered or protected species. No experimentation with animals was performed. No other ethical issues applied to the present research project. Special permissions or rules or sacrificing fish, from Institutional Animal Care and Use Committee (IACUC) or equivalent animal ethics committees, are at present non-existing in Spain for scientific fish sampling. Normally, the process of trawling and handling until biological sampling would lead to high mortality of the fish. Hence, fish were collected without unnecessary suffering.

Received: 27 April 2020; Accepted: 5 October 2020

Published online: 19 October 2020

**References**

1. FAO. *The State of Mediterranean and Black Sea Fisheries* (FAO, Rome, 2018).

2. Coll, M., Albo-Puigserver, M., Navarro, J., Palomera, I. & Dambacher, J. M. Who is to blame? Plausible pressures on small pelagic fish population changes in the northwestern Mediterranean Sea. *Mar. Ecol. Prog. Ser.* **617–618**, 277–294 (2019).

3. Palomera, I. et al. Small pelagic fish in the NW Mediterranean Sea: an ecological review. *Prog. Oceanogr.* **74**, 377–396 (2007).

4. Coll, M., Palomera, I., Tudela, S. & Dowd, M. Food-web dynamics in the South Catalan Sea ecosystem (NW Mediterranean) for 1978–2003. *Ecol. Modell.* **217**, 95–116 (2008).

5. Cardona, L., Martínez-Iñigo, L., Mateo, R. & González-Solis, J. The role of sardine as prey for pelagic predators in the western Mediterranean Sea assessed using stable isotopes and fatty acids. *Mar. Ecol. Prog. Ser.* **531**, 1–14 (2015).

6. Quattrocchi, F. Modelling the relationships of medium and long-term variations of the Anchovy and Sardine catches in the Catalan Sea (NW Mediterranean) with the environmental drivers. (Polytechnique University of Catalunya & Instituto de Ciencias del Mar (ICM), CSIC, 2017).

7. Van Beveren, E. et al. The fisheries history of small pelagics in the Northern Mediterranean. *ICES J. Mar. Sci.* **73**, 1474–1484 (2016).

8. SAC-WGSASP. *Technical Report of the Working Group on Stock Assessment of Small Pelagic Species (WGSASP)*. (2018).

9. Coll, M. & Bellido, J. M. Evaluation of the population status and specific management alternatives for the small pelagic fish stocks in the Northwestern Mediterranean Sea (SPELMED). (2018).

10. Pennino, M. G. et al. Current and future influence of environmental factors on small pelagic fish distributions in the Northwest-ern Mediterranean Sea. *Front. Mar. Sci.* https://doi.org/10.3389/fmars.2020.00622 (2020).

11. Brosset, P. et al. Spatio-temporal patterns and environmental controls of small pelagic fish body condition from contrasted Mediterranean areas. *Prog. Oceanogr.* **151**, 149–162 (2017).

12. Saraux, C. et al. Small pelagic fish dynamics: a review of mechanisms in the Gulf of Lions. *Deep. Res. Part II* (159), 52–61 (2018).
13. Costalago, D. & Palomera, I. Feeding of European pilchard (Sardina pilchardus) in the northwestern Mediterranean: from late larvae to adults. *Sci. Mar.* **78**, 41–54 (2014).

14. Le Bourg, B. et al. Trophic niche overlap of sprat and commercial small pelagic teleosts in the Gulf of Lions (NW Mediterranean Sea). *J. Sea Res.* **103**, 138–146 (2015).

15. Nikolioudakis, N., Palomera, I., Machias, A. & Somarakis, S. Diel feeding intensity and daily ration of the sardine Sardina pilchardus. *Mar. Ecol. Prog. Ser.* **437**, 215–228 (2011).

16. Nikolioudakis, N., Isari, S., Pitta, P. & Somarakis, S. Diet of sardine Sardina pilchardus: an 'end-to-end' field study. *Mar. Ecol. Prog. Ser.* **453**, 173–188 (2012).

17. Nikolioudakis, N., Isari, S. & Somarakis, S. Trophodynamics of anchovy in a non-upwelling system: direct comparison with sardine. *Mar. Ecol. Prog. Ser.* **500**, 215–229 (2014).

18. Tudela, S. & Palomera, I. Diel feeding intensity and daily ration in the anchovy Engraulis encrasicolus in the northwest Mediterranean Sea during the spawning period. *Mar. Ecol. Prog. Ser.* **129**, 55–61 (1995).

19. Tudela, S. & Palomera, I. A Trophic ecology of the European anchovy Engraulis encrasicolus in the Catalan Sea (northwest Mediterranean). *Mar. Ecol. Prog. Ser.* **160**, 121–134 (1997).

20. Brosset, P. et al. Linking small pelagic dietary shifts with ecosystem changes in the Gulf of Lions. *Mar. Ecol. Prog. Ser.* **554**, 157–171 (2016).

21. Albo-Puigserver, M. Ecological and functional role of small and medium pelagic fish in the northwestern Mediterranean Sea. (Polytechnique University of Catalunya & Instituto de Ciencias del Mar (ICM), CSIC, 2019).

22. Amundsen, P. & Sánchez-Hernández, J. Feeding studies take guts: critical review and recommendations of methods for stomach contents analysis in fish. *J. Fish Biol.* **95**, 1364–1373 (2019).

23. Jakubavičiūtė, E., Bergström, U., Eklöf, J. S., Haenel, Q. & Bourlat, S. J. DNA metabarcoding reveals diverse diet of the three-spined stickleback in a coastal ecosystem. *PLoS ONE* **12**, e0186929 (2017).

24. Nielsen, J. M., Clare, E. L., Hayden, B., Brett, M. T. & Kratina, P. Diet tracing in ecology: Method comparison and selection. *Methods Ecol. Evol.* **9**, 278–291 (2018).

25. Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents. *Front. Zool.* **10**, 34 (2013).

26. Berry, O. et al. Comparison of morphological and DNA metabarcoding analyses of diets in exploited marine fishes. *Mar. Ecol. Prog. Ser.* **540**, 167–181 (2015).

27. Casey, J. M. et al. Reconstructing hyperdiverse food webs: Gut content metabarcoding as a tool to disentangle trophic interactions on coral reefs. *Methods Ecol. Evol.* **10**, 1157–1170 (2019).

28. Thomas, A. C., Deagle, B. E., Eveson, J. P., Harsch, C. H. & Trites, A. W. Quantitative DNA metabarcoding: improved estimates of species proportional biomass using correction factors derived from control material. *Mol. Ecol. Resour.* **16**, 714–726 (2016).

29. Baker, R., Buckland, A. & Sheaves, M. Fish gut content analysis: Robust measures of diet composition. *Fish Fish* **15**, 170–177 (2014).

30. Pacella, S. R. et al. Incorporation of diet information derived from Bayesian stable isotope mixing models into mass-balanced marine ecosystem models: A case study from the Marennes-Oléron Estuary France. *Ecol. Modell.* **267**, 127–137 (2013).

31. Parnell, A. C. et al. Bayesian stable isotope mixing models. *Environmentmetrics* **24**, 387–399 (2013).

32. Stock, B. C. et al. Analyzing mixing systems using a new generation of Bayesian tracer mixing models. *PeerJ* **6**, e5096 (2018).

33. Carreon-Martinez, L., Johnson, T. B., Ludsin, S. A. & Heath, D. D. Utilization of stomach content DNA to determine diet diversity in piscivorous fishes. *J. Fish Biol.* **78**, 1170–1182 (2011).

34. Hunter, J. R. & Kimbrell, C. A. Egg cannibalism in the Northern anchovy Engraulis mordax. *Fish. Bull.* **78**, 811–816 (1980).

35. Vander Zanden, M. J., Clayton, M. K., Moody, E. K., Solomon, C. T. & Weidel, B. C. Stable isotope turnover and half-life in animal tissues: a literature synthesis. *PLoS ONE* **10**, e0116182 (2015).

36. Bearhop, S., Adams, C. E., Waldron, S., Fuller, R. A. & MacLeod, H. Determining trophic niche width: a novel approach using stable isotope analysis. *J. Anim. Ecol.* **73**, 1007–1012 (2004).

37. Jackson, A. L., Inger, R., Parnell, A. C. & Bearhop, S. Comparing isotopic niche widths among and within communities: SIBER—Stable Isotope Bayesian Ellipses in R. *J. Anim. Ecol.* **80**, 595–602 (2011).

38. Cannas, R. et al. Report on bioinformatic analyses of the GBS data and report of the population genetic analyses. *Evaluation of the population status and specific management alternatives for the small pelagic fish stocks in the Northwestern Mediterranean Sea* (SPELMED). (2018).

39. Koleff, P., Gaston, K. J. & Lennon, J. J. Measuring beta diversity for presence-absence data. *J. Anim. Ecol.* **72**, 367–382 (2003).

40. Southwood, T. R. E. & Henderson, P. A. *Ecological Methods* (Backwell Science, USA, 2000).

41. Pianka, E. R. Niche Overlap and Diffuse Competition. *Proc. Natl. Acad. Sci. U. S. A.* **71**, 2141–2145 (1974).

42. Bachiller, E. & Irigoien, X. Allometric relations and consequences for feeding in small pelagic fish in the Bay of Biscay. *ICES J. Mar. Sci.* **70**, 232–243 (2013).

43. Bachiller, E. & Irigoien, X. Trophodynamics and diet overlap of small pelagic fish species in the Bay of Biscay. *Mar. Ecol. Prog. Ser.* **534**, 179–198 (2015).

44. Raab, K. et al. Anchovy Engraulis encrasicolus diet in the North and Baltic Seas. *J. Sea Res.* **65**, 131–140 (2011).

45. Karachle, P. K. & Stergiou, K. I. Feeding and ecomorphology of three clupeoids in the N Aegean Sea. *Mediterr. Mar. Sci.* **15**, 9 (2013).

46. Costalago, D., Garrido, S. & Palomera, I. Comparison of the feeding apparatus and diet of European sardines Sardina pilchardus of Atlantic and Mediterranean waters: ecological implications. *J. Fish Biol.* **86**, 1348–1362 (2015).

47. Costalago, D., Navarro, J., Álvarez-Calleja, I. & Palomera, I. Ontogenetic and seasonal changes in the feeding habits and trophic levels of two small pelagic fish species. *Mar. Ecol. Prog. Ser.* **460**, 169–181 (2012).

48. Garrido, S. et al. Trophic ecology of pelagic fish species off the Iberian coast: diet overlap, cannibalism and intraguild predation. *Mar. Ecol. Prog. Ser.* **539**, 271–285 (2015).

49. Bacha, M. & Amara, R. Spatial, temporal and ontogenetic variation in diet of anchovy (Engraulis encrasicolus) on the Algerian coast (SW Mediterranean). *Estuar. Coast. Shelf Sci.* **85**, 257–264 (2009).

50. Brosset, P. et al. Body reserves mediate trade-offs between life-history traits: New insights from small pelagic fish reproduction. *R. Soc. Open Sci.* **3**, 160202 (2016).

51. Bachiller, E., Cotano, U., Boyra, G. & Irigoien, X. Spatial distribution of the stomach weights of juvenile anchovy (*Engraulis encrasicolus* L.) in the Bay of Biscay. *ICES J. Mar. Sci.* **70**, 362–378 (2013).

52. Ventero, A., Iglesias, M. & Villamor, B. Anchovy (Engraulis encrasicolus) otoliths reveal growth differences between two areas of the Spanish Mediterranean Sea. *Sci. Mar.* **81**, 327 (2017).

53. Costalago, D., Palomera, I. & Tirelli, V. Seasonal comparison of the diets of juvenile European anchovy *Engraulis encrasicolus* and sardine *Sardina pilchardus* in the Gulf of Lions. *J. Sea Res.* **89**, 64–72 (2014).

54. Borne, D., Tirelli, V., Brandt, S. B., Umani, S. F. & Arneri, E. Diet of Engraulis encrasicolus in the northern Adriatic Sea (Mediterranean): ontogenetic changes and feeding selectivity. *Mar. Ecol. Prog. Ser.* **392**, 193–209 (2009).

55. Layman, C. A. et al. Applying stable isotopes to examine food-web structure: an overview of analytical tools. *Biol. Rev.* **87**, 545–562 (2012).
56. Albo-Puigserver, M., Navarro, J., Coll, M., Layman, C. A. & Palomera, I. Trophic structure of pelagic species in the northwestern Mediterranean Sea. *J. Sea Res.* **117**, 27–35 (2016).

57. Petta, J. C. et al. Are you really what you eat? Stomach content analysis and stable isotope ratios do not uniformly estimate dietary niche characteristics in three marine predators. *Oecologia* **192**, 1111–1126 (2020).

58. Demestre, M. Growth and distribution of Solenocera membranacea (Risso, 1816) (Decapoda, Dendrobranchiata) in the north-western Mediterranean Sea. *Sci. Mar.* **57**, 161–166 (1993).

59. Tsagarakis, K., Giannoulaki, M., Somarakis, S. & Machias, A. Variability in positional, energetic and morphometric descriptors of European anchovy *Engraulis encrasicolus* schools related to patterns of diurnal vertical migration. *Mar. Ecol. Prog. Ser.* **446**, 243–258 (2012).

60. Coll, M., Shannon, L. J., Moloney, C. L., Palomera, I. & Tudela, S. Comparing trophic flows and fishing impacts of a NW Mediterranean ecosystem with coastal upwelling systems by means of standardized models and indicators. *Ecol. Modell.* **198**, 53–70 (2006).

61. Checkley, D. M., Asch, R. G. & Rykaczewski, R. R. Climate, Anchovy, and Sardine. *Ann. Rev. Mar. Sci.* **9**, 469–493 (2017).

62. Irigoien, X. & De Roos, A. The role of intraguild predation in the population dynamics of small pelagic fish. *Mar. Biol.* **158**, 1683–1690 (2011).

63. Bachiller, E., Cotano, U., Ibaibarriaga, L., Santos, M. & Irigoien, X. Intraguild predation between small pelagic fish in the Bay of Biscay: impact on anchovy (*Engraulis encrasicolus* L.) egg mortality. *Mar. Biol.* **162**, 1351–1369 (2015).

64. Purcell, J. E. Jellyfish and ctenophore blooms coincide with human proliferations and environmental perturbations. *Ann. Rev. Mar. Sci.* **4**, 209–235 (2012).

65. Naman, S. M. et al. Stable isotope-based trophic structure of pelagic fish and jellyfish across natural and anthropogenic landscape gradients in a fjord estuary. *Ecol. Evol.* **6**, 8159–8173 (2016).

66. Albo-Puigserver, M. et al. Trophic ecology of range-expanding round sardinella and resident sympatric species in the NW Mediterranean. *Mar. Ecol. Prog. Ser.* **620**, 139–154 (2019).

67. Morote, E., Olivar, M. P., Villate, F. & Uriarte, I. A comparison of anchovy (*Engraulis encrasicolus*) and sardine (*Sardina pilchardus*) larvae feeding in the Northwest Mediterranean: influence of prey availability and ontogeny. *ICES J. Mar. Sci.* **67**, 897–908 (2010).

68. Robinson, M. L., Gomez-Raya, L., Rauw, W. M. & Peacock, M. M. Fulton's body condition factor K correlates with survival time in a thermal challenge experiment in juvenile Lahontan cutthroat trout (Oncorhynchus clarki henshawi). *J. Therm. Biol.* **33**, 363–368 (2008).

69. Siokou-Frangou, I. et al. Plankton in the open Mediterranean Sea: a review. *Biogeosciences* **7**, 1543–1586 (2010).

70. Vila, M. & Masó, M. Phytoplankton functional groups and harmful algal species in anthropogenically impacted waters of the NW Mediterranean Sea. *Sci. Mar.* **69**, 31–45 (2005).

71. Percopo, I., Siano, R., Cerino, F., Sarno, D. & Zingone, A. Phytoplankton diversity during the spring bloom in the northwestern Mediterranean Sea. *Bot. Mar.* **54**, 243–267 (2011).

72. Devloo-Delva, F. et al. How does marker choice affect your diet analysis: comparing genetic markers and digestion levels for diet metabarcoding of tropical-reef piscivores. *Mar. Freshw. Res.* **70**, 8–18 (2019).

73. Forin-Wiart, M.-A. et al. Evaluating metabarcoding to analyse diet composition of species foraging in anthropogenic landscapes using Ion Torrent and Illumina sequencing. *Sci. Rep.* **8**, 17091 (2018).

74. Thuo, D. et al. Food from faeces: evaluating the efficacy of scat DNA metabarcoding in dietary analyses. *PLoS ONE* **14**, e0225805 (2019).

75. Coll, M., Pennino, M. G., Steenbeek, J., Sole, J. & Bellido, J. M. Predicting marine species distributions: complementarity of food-web and Bayesian hierarchical modelling approaches. *Ecol. Modell.* **405**, 86–101 (2019).

76. Feuilloley, G. et al. Concomitant changes in the environment and small pelagic fish community of the Gulf of Lions. *Prog. Oceanogr.* **186**, 102375 (2020).

77. Pennino, M. G. et al. Ingestion of microplastics and occurrence of parasite association in Mediterranean anchovy and sardine. *Mar. Pollut. Bull.* **158**, 111399 (2020).

78. Compa, M., Ventero, A., Iglesias, M. & Deudero, S. Ingestion of microplastics and natural fibres in Sardina pilchardus (Walbaum, 1792) and Engraulis encrasicolus (Linnaeus, 1758) along the Spanish Mediterranean coast. *Mar. Pollut. Bull.* **128**, 89–96 (2018).

79. Bertrand, J., Leonori, I., Première, P. Y. & Cosimi, G. The general specifications of the MEDITs surveys. *Sci. Mar.* **66**, 9–17 (2002).

80. Bertrand, J. A., De Sola, L. G., Papaconstantinou, C., Relini, G. & Souplet, A. The general specifications of the MEDITs surveys. *Sci. Mar.* **66**, 9–17 (2002).

81. Cole, M. et al. Isolation of microplastics in biota-rich seawater samples and marine organisms. *Sci. Rep.* **4**, 4528 (2014).

82. Hyslop, E. J. Stomach contents analysis-a review of methods and their application. *J. Fish Biol.* **17**, 411–429 (1980).

83. Bachiller, E., Skaret, G., Nottestad, L. & Slotte, A. Feeding ecology of Northeast Atlantic mackerel, Norwegian spring-spawning herring and blue whiting in the Norwegian Sea. *PLoS ONE* **11**, e0149238 (2016).

84. Somarakis, S. et al. Daily egg production of anchovy in European waters. *ICES J. Mar. Sci.* **61**, 944–958 (2004).

85. Palomera, I., Tejeiro, B. & Alemany, F. *Size at first maturity of the NW Mediterranean anchovy*. (2003).

86. Silva, A. et al. Temporal and geographic variability of sardine maturity at length in the northeastern Atlantic and the western Mediterranean. *ICES J. Mar. Sci.* **63**, 663–676 (2006).

87. Lawlor, L. R. Overlap, similarity, and competition coefficients. *Ecology* **61**, 245–251 (1980).

88. Barroeta, Z., Olivar, M. P. & Palomera, I. Energy density of zooplankton and fish larvae in the southern Catalan Sea (NW Mediterranean). *J. Sea Res.* **124**, 1–9 (2017).

89. Sabatés, A. Distribution pattern of larval fish populations in the Northwestern Mediterranean. *Mar. Ecol. Prog. Ser.* **59**, 75–82 (1990).

90. Wangensteen, O. S., Palacín, C., Guardiola, M. & Turon, X. DNA metabarcoding of littoral hard-bottom communities: high diversity and database gaps revealed by two molecular markers. *PeerJ* **6**, e4705 (2018).

91. Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents. *Front. Zool.* **10**, 34 (2013).

92. Vasselon, V., Rimet, F., Tapolczai, K. & Bouchez, A. Assessing ecological status with diatoms DNA metabarcoding: Scaling-up on a WFD monitoring network (Mayotte island, France). *Ecol. Indic.* **82**, 1–12 (2017).

93. Vierna, J., Doña, J., Vizcaino, A., Serrano, D. & Jovani, R. PCR cycles above routine numbers do not compromise high-throughput DNA barcoding results. *Genome* **60**, 868–873 (2017).

94. Magoc, T. & Salzberg, S. L. FLASH: fast length adjustment of short reads to improve genome assemblies. *Bioinformatics* **27**, 2957–2963 (2011).

95. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. *EMB net.journal* **17**, 10 (2011).

96. Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. *Nat. Methods* **7**, 335–336 (2010).

97. Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: a versatile open source tool for metagenomics. *PeerJ* **4**, e2584 (2016).

98. Mahé, F., Rognes, T., Quince, C., de Vargas, C. & Dunthorn, M. Swarm v2: highly-scalable and high-resolution amplicon clustering. *PeerJ* **3**, e1420 (2015).
99. Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C. & Knight, R. UCHIME improves sensitivity and speed of chimera detection. *Bioinformatics* **27**, 2194–2200 (2011).

100. Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. *Bioinformatics* **26**, 2460–2461 (2010).

101. Wang, Q., Garrity, G. M., Tiedje, J. M. & Cole, J. R. Naive bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. *Appl. Environ. Microbiol.* **73**, 5261–5267 (2007).

102. Bokulich, N. A. et al. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. *Nat. Methods* **10**, 57–59 (2013).

103. Esling, P., Lejzerowicz, F. & Pawlowski, J. Accurate multiplexing and filtering for high-throughput amplicon-sequencing. *Nucleic Acids Res.* **43**, 2513–2524 (2015).

104. Illumina, Inc [internet]. Effects of index Misassignment on multiplexing and downstream analysis (2017). Available from: https://www.illumina.com/. Accessed October 2020.

105. National Center for Biotechnology Information (NCBI) [internet]. Bethesda (MD): National Library of Medicine (US), National Center for Biotechnology Information (2018). Available from: https://www.ncbi.nlm.nih.gov/. Accessed August 2020.

106. Coplen, T. B. Guidelines and recommended terms for expression of stable-isotope-ratio and gas-ratio measurement results. *Rapid Commun. Mass Spectrom.* **25**, 2538–2560 (2011).

107. Post, D. M. et al. Getting to the fat of the matter: models, methods and assumptions for dealing with lipids in stable isotope analyses. *Oecologia* **152**, 179–189 (2007).

108. Hastie, T. & Tibshirani, R. Generalized Additive Models. *Stat. Sci.* **1**, 297–318 (1986).

109. Zuur, A. F., Ieno, E. N. & Elphick, C. S. A protocol for data exploration to avoid common statistical problems. *Methods Ecol. Evol.* **1**, 3–14 (2010).

110. Fox, J. & Weisberg, S. *An R Companion to Applied Regression* (Sage, Thousand Oaks, 2011).

111. Wood, S. N. Generalized additive models: an introduction with R. *J. Stat. Softw.* **16**, 2 (2006).

112. R Core Team. R: A Language and Environment for Statistical Computing. (2019).

113. Wickham, H. ggplot2: *elegant graphics for data analysis*. Springer (Springer, 2009).

114. Gotelli, N. J., Hart, E. M. & Ellison, A. M. EcoSimR: Null model analysis for ecological data. (2015). doi:https://doi.org/10.5281/zenodo.16522

115. Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. *J. R. Stat. Soc. Ser. B (Statistical Methodol).* **73**, 3–36 (2011).

116. Wei, T. et al. Package 'corrplot'. *Am. Stat.* **56**, 316–324 (2017).

117. QGIS Development Team. QGIS (Version 3.2.1-Bonn). (2018).

118. Hsieh, T. C., Ma, K. H., Chao, A. & McInerny, G. iNEXT: an R package for rarefaction and extrapolation of species diversity (ill numbers). *Methods Ecol. Evol.* **7**(12), 1451–1456 (2016).

## Acknowledgements

All gut content characterization analyses under the microscope were carried out in the laboratory of AZTI (Pasaia, Spain), under the supervision of U. Cotano (Marine Ecosystem Functioning Area, AZTI). This study was partially funded by the projects PELWEB (ES-PN-2017-CTM 2017-88939-R) and PELCAT (CAT 152CAT00013, TAIS ARP059/19/00005). MEDITS data collection has been co-funded by the EU through the European Maritime and Fisheries Fund (EMFF) within the National Program of collection, management and use of data in the fisheries sector and support for scientific advice regarding the Common Fisheries Policy. The authors express their gratitude to all the people that worked in the MEDITS surveys.

## Author contributions

M.C., M.A.P., J.G., J.M.B and E.B. contributed to the conception and design of the experiment. E.B., J.G., M.A.P. and N.M.M. conducted the experiments. E.B., M.C., M.A.P., J.G. and N.M.M. wrote the experimental part in the manuscript. E.B., M.A.P., J.G., M.G.P. and M.C. analysed the results. All authors reviewed and contributed to the writing of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-74602-y.

Correspondence and requests for materials should be addressed to E.B.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020