Ontogenetic, Spatial and Inter-Annual Variability in the Diet of European Hake Merluccius merluccius Linnaeus, 1758, in the North Aegean Sea

Abstract

This study contributes to filling knowledge gaps regarding recent information on the diet of the European hake, Merluccius merluccius Linnaeus, 1758, in the Greek seas, particularly its ontogenetic shifts and its spatiotemporal variability. The trophic preferences of M. merluccius were investigated in the North Aegean Sea during the summers of 2019 to 2023 with visual stomach content analysis to assess its composition, diversity, and variability across body-size classes, years, and subareas. The identified prey are functionally diverse and in many cases also primary targets for local fisheries. The ontogenetic trophic niche of hake was characterized by two distinct shifts, delineated by 10 and 50 cm body-size thresholds. Cephalopods were a prevalent dietary component for large hake individuals. The intermediate body-size classes demonstrated greater trophic niche breadth concerning prey diversity and absolute prey-size ranges. A feeding strategy characterized by the specialization of individuals was also revealed. The effect of the temporal and spatial context on the hake diet was occasionally correlated with spatiotemporal variations in the sizes of their populations. In conclusion, ontogenetic and spatiotemporal variability in the hake diet were found to be significant and should be considered in the data collection design and analyses of the trophic interactions of the species.

1. Introduction

Merluccius merluccius Linnaeus, 1758 (European hake, hereafter hake) is an important demersal and benthopelagic predator of the continental shelf and upper-slope ecosystems of the Mediterranean [1]. It is also one of the main fisheries-targeted species in its distribution range [2]. Hake has been reported to be the most overexploited stock across the Mediterranean Sea [3], where its landings have decreased from 52,394 tonnes in 1994 to 17,824 tonnes in 2021 [4]. Moreover, undersized individuals still constitute a large part of hake catches and are either discarded or landed [5], despite the endorsement of legislation regulating minimum conservation reference size and minimum trawl codend mesh sizes for hake [6,7]. However, hake stock status in the Mediterranean has been assessed to be slowly improving at a regional scale, demonstrating a stable average F/FMSY overexploitation ratio [4]. In Greece, hake is among the main fisheries-targeted demersal species [8,9]. Its abundance is higher at depths of 100–200 m, while it is fished down to 550 m [10], principally in the North and Western Aegean Sea, mainly with trawlers, as well as longliners and netters [10]. Hake production decreased in the Aegean Sea from 4867 tonnes in 1994 to 2744 tonnes in 2021 [11] and the status of hake stock in the Aegean Sea has been characterized as bad (ratio of observed biomass to the biomass that would provide maximum sustainable yield, B/BMSY = 0.52) and overexploited (ratio of observed fishing mortality to the fishing mortality that would produce maximum sustainable yield, F/FMSY = 1.57) [12].

Hake, an active predator of mobile and patch/school-forming prey [2,10,13], exhibits dynamic feeding habits. It engages, particularly young individuals, in diel vertical movements to pursue migrating prey, thus linking the pelagic and benthic domains [10]. Research has also indicated significant ontogenetic variation in hake’s diet, alongside seasonal, bathymetric, and geographic variability [10,14]. Numerous studies have explored the trophic preferences of hake in the Mediterranean, investigating both ontogenetic diet changes and spatial variability [15,16,17,18,19]. However, there remains a notable paucity of published time-series data on hake diet variability in the Mediterranean, with few exceptions [20]. Furthermore, recent research using DNA metabarcoding [21,22] or stable isotope analysis [23] has shed light on the intricate trophic relationships of Mediterranean hake, offering insights beyond what visual stomach contents analysis is capable of revealing. While the diet of hake has been extensively studied in the Greek seas [10,24,25,26,27], relevant literature has become scarce since the 2000s. Moreover, data on the temporal variability of its diet in the Greek seas are not available. Similarly, besides qualitative insights [28] and data on the diet of young hake individuals [29], detailed quantitative data on ontogenetic shifts in the diet of hake in the area are lacking in the peer-reviewed literature.

The ecosystem-based approach to fisheries management (EBFM) [30,31] has received support from the EU through initiatives such as the Common Fisheries Policy (CFP) [32], the Marine Strategy Framework Directive (MSFD) [33], and the Marine Spatial Planning Directive (MSPD) [34]. Accordingly, EU member states have been gathering data on trophic interactions of priority species within the multiannual EU program for the collection, management, and utilization of fisheries sector data (EU MAP). This EU MAP activity aims to inform assessments of the impacts of fisheries on biological marine resources, e.g., stock assessments and ecosystem models, with data on key components of the food web regarding their trophic interactions and natural mortalities. The collection of stomach samples was integrated into EU MAP as a pilot study from 2019 until 2021, and has become a regular component since 2022. Given its significance as both a primary target in Mediterranean fisheries and a key component in the marine food webs of the region, acting as both predator and competitor for species targeted by fisheries [2], hake has been selected as the first target species for stomach sampling in the Mediterranean under EU MAP.

This study quantifies and evaluates the ontogenetic, spatial, and inter-annual variability of the trophic preferences of hake in the North Aegean Sea, using a spatiotemporal resolution that has not been offered by previous studies in the area. The aim of the study is to provide comprehensive and up-to-date information that could enhance our understanding of the trophic interactions of hake in the marine food webs in the area. The assessment was based on the visual analysis of the contents of hake stomachs collected in the North Aegean Sea during the summers of 2019 to 2023 within the framework of EU MAP.

2. Materials and Methods

2.1. Study Area

The study area comprised the North Aegean Sea (Figure 1), which is spatially divided into five subareas according to the Greek multiannual program: Northeastern Aegean Sea (CHIO-MIT), Evia (EVIA), Thermaikos Gulf and adjacent seas (THERM), Thracian Sea (THR-LIM), and Thessaly and Sporades Islands (VOL-SPOR). Commercial fisheries in the North Aegean Sea are characterized by high catches of small pelagic fish predominantly captured by purse seiners. Bottom trawls are employed to target hake, mullets, shrimps, and cephalopods, while multi-gear small-scale coastal fisheries target a variety of pelagic and demersal species [35]. In 2022, the total quantity of small- and medium-scale fisheries landings in the North and Northeastern Aegean Sea reached 35,636 tonnes, constituting 55% of the overall Greek fisheries landings. The highest landings took place in the Thracian Sea (16,299 tonnes) and the area of Thermaikos Gulf and adjacent seas (14,121 tonnes), while landings in the Northeastern Aegean Sea amounted to 3630 tonnes [36]. Hake landings by trawlers in the North Aegean Sea during 2019–2022 ranged from 1318 to 1370 tonnes, corresponding to 59–66% of the hake landings by trawlers in Greece. Total hake production in the North Aegean Sea in 2022 was 2135 tonnes, 51% of the total hake production in Greece [36].

2.2. Field Work

The hake samples were collected in the International Bottom Trawl Survey in the Mediterranean (MEDITS) samplings conducted by the Fisheries Research Institute in the North Aegean Sea during the summers of 2019 to 2023 (June–July) in the framework of EU MAP (Figure 1). Fishing was carried out with an experimental bottom trawl with a codend net mesh opening of 40 mm, and a standard cruising speed of 3 knots. The duration of the hauls was 30 min for stations shallower than 200 m, and 60 min at deeper stations. All hake samples were collected during daylight hours, from late in the morning to early in the afternoon.

Sampling and lab analysis of the stomachs were implemented using the protocols and guidelines originally recommended in deliverable D3.3 of the MARE/2014/19 “Strengthening regional cooperation in the area of fisheries data collection in the Mediterranean and Black Sea” project [37] and further elaborated in the report of the RCG-Med&BS-organized “Workshop on Sampling, Processing, and Analyzing the Stomach Contents (WKSTCON)” that was held at Palma de Mallorca, Spain, in 2018 [38], and in deliverable D4.1 of the MARE/2016/22 “Strengthening regional cooperation in the area of fisheries data collection in the Mediterranean and Black Sea (STREAM)” project [39].

Total length (mm) and total weight (g) of hake individuals were measured onboard. In each annual survey, hake individuals were sampled for every one of six 100 mm (TL) body-size classes: 0 (<100 mm), 1 (100–199 mm), 2 (200–299 mm), 3 (300–399 mm), 4 (400–499 mm), and 5 (>500 mm). Sample collection aimed to reach an annual total of at least 100 stomachs. The fish were dissected, and non-everted stomachs were removed and preserved in a freezer at −20 °C until their analysis in the lab.

2.3. Lab Work

In the lab, the stomachs were partitioned between those that contained prey items (“full stomachs”), and those that did not (“empty stomachs”). The stomachs that contained only mucus that was over 0.2% of body weight were considered full [40] and were included as such in the calculation of the vacuity and repletion indices. The stomachs were weighed and then dissected for the removal of prey items, which were identified to the lowest possible taxonomic level, counted, and weighed to the nearest 0.01 g. Taxonomic identification of the prey items was carried out macroscopically, or with the aid of a stereomicroscope for prey items of small body sizes. Furthermore, the identification of fish species was aided by microscopic examination of the otoliths.

2.4. Data Analysis

To evaluate whether the number of stomachs was sufficient for a valid characterization of hake diet, prey accumulation curves were computed [41,42].

Feeding intensity was assessed with the vacuity index [43] and the repletion index [44]. The vacuity index (VI) was calculated as the percentage of the analyzed stomachs that were empty, whereas the repletion index (RI) was calculated as the ratio of stomach contents weight to body weight, expressed as a per-thousand value. The χ² test of independence was used to test for the null hypothesis that the proportions of full and empty stomachs were independent of body-size class, year, or subarea. Permutational multivariate analysis of variance (PERMANOVA) [45] was used to test for statistically significant differences in stomach repletion between body-size classes, years, and subareas. PERMANOVA was run on a Euclidean similarity matrix of the hake individuals computed on the RI index values. All factors were treated as fixed in the analysis design. The sum of squares type selected was Type III (partial), and the permutation method involved a permutation of residuals under a reduced model for the three-way layout or unrestricted permutation of raw data for the one-way layout. The number of permutations was set at 9999.

PERMANOVA was also used to test for statistically significant differences in the composition of the diet of hake, both by numbers and by weight, between years, subareas, and body-size classes. Three-way layout PERMANOVA was run on a Bray–Curtis similarity matrix of the hake stomachs computed on the prey numbers or weights values. All factors were treated as fixed in the analysis design. The sum of squares type selected was Type III (partial), and the permutation method involved a permutation of residuals under a reduced model. The number of permutations was set at 9999.

The proportions of prey taxa in the stomach contents were calculated with the frequency of occurrence, and the numerical and gravimetric methods [46]:

The frequency of occurrence of a prey taxon in the stomach contents was calculated as:

$$\%F_i=\left(S_i/S_F\right) \times 100$$

(1)

where S_i is the number of stomachs that contained prey i, and S_F is the number of full stomachs.

The percentage of the total number of prey individuals in the stomach contents belonging to a prey taxon was calculated as:

$$\%N_i=(n_i/{\displaystyle \sum }_{i=1}^Q{n}_i)\times 100$$

(2)

where n_i is the total number of prey taxon i individuals and ${\displaystyle \sum }_{i=1}^Q{n}_i$ is the total number of individuals of all prey taxa (Q) in all full stomachs.

The percentage of the total weight of prey in the stomach contents belonging to a prey taxon was calculated as:

$$\%W_i=(w_i/{\displaystyle \sum }_{i=1}^Q{w}_i)\times 100$$

(3)

where w_i is the total wet weight of prey taxon i and ${\displaystyle \sum }_{i=1}^Q{w}_i$ is the total wet weight of all prey taxa (Q) in all full stomachs.

The relative importance of each prey taxon in the diet of hake was estimated with the calculation of the % form of the Index of Relative Importance (IRI) [47]. This index combines the information provided by the %F, %N, and %W indices in a single value. A modified version of IRI [48] was used, which is calculated as:

$$IR{I}_i=\left(N_i+W_i\right)\times F_i$$

(4)

The calculation of the % form of IRI was carried out as follows [49]:

$$\%IR{I}_i=(IR{I}_i/{{\displaystyle \sum }}_{i=1}^{Q}IR{I}_i)\times 100$$

(5)

The trophic preferences of hake, including the variations between body-size classes, years, and subareas, was assessed primarily based on the %IRI index values. The data for the %F, %N, and %W indices are provided in Supplementary Table S1 and illustrated in Supplementary Figures S1–S3.

The trophic links between hake size classes and their prey were represented in a prey–predator network graph, aiming to depict the trophic preferences of size classes while distinguishing between shared and size-class specific prey. The network graph allowed the differentiation of prey according to their out degrees, which correspond to the number of hake size classes that used them. It also depicts the strength of the trophic links via variations in edge width. The graph was laid out with the “force-atlas” algorithm.

The richness and diversity of prey taxa in stomach contents were estimated aiming to quantify the width of the trophic niche width of hake and assess its variability between body-size classes, years, and subareas. The indices used measure diversity in units of effective numbers of species; i.e., the number of equally abundant species that would be needed to give the same value of a diversity measure [50,51]. Moreover, diversity accumulation curves [52] for ⁰D (taxon richness), ¹D (the exponential of Shannon entropy), and ²D (the reciprocal of Simpson concentration) were computed with the iNEXT R package [53] to rigorously compare diversities by taking into account the effect of sample size. Differences in diversity were considered significant (at a level of p = 0.05) when the 95% CI bands of the diversity accumulation curves did not overlap [52]. Confidence intervals were computed with bootstrap resampling and 100 bootstrap replications. The methodology of this analysis is presented in detail in Appendix A.

Intraspecific variation in resource use has been argued to have significant ecological effects and be important in niche theory [54,55,56]. Total niche width (or degree of specialization) of a population can be broken down into two components [57]: the individual niche width (within-individual component), and the variance in resource use between individuals (between-individuals component). The relative importance of the components of hake total niche width, as a facet of the feeding strategy of the species, was assessed with the modified Costello graphical analysis [58,59]. In this method, species feeding strategy is graphically depicted using a 2-D representation, where the prey-specific abundance of prey taxon i was plotted against its frequency of occurrence in the stomachs with food contents. The prey-specific abundance of prey taxon i was calculated with the following formula [59]:

$$P_i=\frac{\sum S_i}{\sum S_t}\times 100$$

(6)

where P_i is the prey-specific abundance of prey taxon, S_i is the stomach content in numbers corresponding to prey i, and S_t the total stomach content only in stomachs that contained prey i. This method assesses simultaneously prey importance, feeding strategy, and inter- and intra-individual components of the population trophic niche width. This information is obtained by the visual inspection of the distribution of the points representing prey in the produced plots across the bottom left-top right diagonal (rare vs. dominant prey groups), top-bottom axis (specialization vs. generalization in the diet), and top left-bottom right axis (specialization at the individual vs. at the population level).

Hake diet overlap between body-size classes, years, and subareas was quantified with the Renkonen (or Schoener) proportional overlap index [60]:

$$C_{xy}=1-0.5\times (\sum \left|p_{xi}-p_{yi}\right|)$$

(7)

where p_xi and p_yi are the proportions of prey i (in terms of numerical abundance) in the diet of the species in the size class, year, and area x and y, respectively. C_xy ranges from 1 (same prey items in the same proportions) to 0 (no common prey items). This index is a measure of the actual area of overlap of resource utilization curves and is not sensitive to how resource states are divided [60].

All analyses were performed using the R Statistical Software (v4.3.2) [61] except the prey–predator network graph, which was plotted using the Gephi graph visualization and manipulation software [62], and the PERMANOVA tests, which were carried out in PERMANOVA 1.0.8 [63].

3. Results

3.1. Sampling Effort

A total of 769 hake stomachs were collected during the five annual MEDITS samplings that were conducted in the study area from 2019 to 2023 (Table 1). The number of stomach samples varied between years, ranging from 95 to 232, and between subareas, ranging from 90 to 303. Most stomach samples were collected from hake individuals of body-size classes 1 (225), 2 (261), and 3 (203). Sample sizes for classes 0, 4, and 5 were 26, 40, and 14 stomachs, respectively. The overall length range of hake individuals that had stomachs with prey in their contents was 72–572 mm (TL).

3.2. Feeding Intensity

Out of the 769 collected stomachs, 522 were categorized as full and 247 as empty. Thirty-six stomachs contained only mucus that exceeded 0.2% of body weight. These stomachs were regarded as full in the calculations of the VI and RI indices but were not considered in further analyses. VI values ranged between body-size classes from 25 to 35.71%, between years from 17.11 to 37.89%, and between subareas from 24.11 to 46.41% (Table 1). However, the χ² test of independence failed to reject the respective null hypotheses that the percentages of full and empty stomachs were independent of body-size classes, years, or subareas.

No significant effect of body size, year, or subarea on RI was revealed by PERMANOVA using a three-way layout. However, one-way PERMANOVA revealed that mean RI significantly differed among body-size classes (Pseudo-F = 6.95, p(perm) < 0.001), with mean RI values ranging from 21.05 to 77.10‰ (Table 1). Pairwise PERMANOVA tests indicated that mean RI was significantly higher in class 5 individuals compared to class 0 (t = 2.18, p < 0.05), 2 (t = 3.23, p < 0.01), 3 (t = 4.74, p < 0.001), or 4 individuals (t = 3.33, p < 0.01). The lowest mean RI value was observed in class 3, individuals, which was significantly lower than in class 1 (t = 4.43, p < 0.001), 2 (t = 3.52, p < 0.001), or 5 individuals (t = 4.74, p < 0.001). One-way PERMANOVA did not reveal any significant differences in mean RI between years or between subareas.

3.3. Sufficiency of the Number of Examined Fish Stomachs

The prey accumulation curves analysis results suggest that the sampling effort was not entirely sufficient for any of the body-size classes or for all stomachs pooled (Figure 2).

3.4. Diet Composition

The prey items found in the stomach contents analysis were classified into 58 taxa (Table 2). Among these, 40 taxa were identified at the species level, 6 at the genus level, and 7 at the family level, while the remaining 5 taxa were classified at higher taxonomic levels.

Among the identified prey taxa, 28 were benthic, demersal, or bathydemersal, while 21 were pelagic or bathypelagic. Additionally, four taxa were benthopelagic. The bathymetric ranges of the prey taxa varied, ranging from shallow-water continental shelf species (e.g., Sepiola sp., Alpheus glaber, Solenocera crassicornis, Lesueurigobius friesii, Sardina pilchardus, Serranus hepatus, Spicara flexuosum) to bathypelagic species (e.g., Pasiphaea multidentata, Sergestidae, Micromesistius poutassou, Nettastoma melanura, Myctophidae, Stomiidae), and bathydemersal species (Argentina sphyraena, Chlorophthalmus agassizi).

The identified prey taxa were also categorized into 10 functional groups based on their taxonomy, ecological roles, and preferred habitat (Table 2). These groups comprised octopuses and cuttlefish, squids, benthic decapods, bathypelagic decapods, benthopelagic decapods, non-decapod crustaceans, demersal fish, bathydemersal fish, bathypelagic fish, and benthopelagic fish.

The factors that had a significant effect on the composition by weight of the diet of hake were year (Pseudo-F = 1.6022, p(perm) = 0.0249), body-size class (Pseudo-F = 1.8251, p(perm) = 0.0007), and the interactions year × subarea (Pseudo-F = 1.3094, p(perm) = 0.006), year × body-size class (Pseudo-F = 1.2923, p(perm) = 0.0075), subarea × body-size class (Pseudo-F = 1.2124, p(perm) = 0.0004), and year × subarea × body-size class (Pseudo-F = 1.3875, p(perm) = 0.0006). Regarding diet composition by numbers, the factors that had a significant effect were body-size class (Pseudo-F = 1.6532, p(perm) = 0.0078), and the interactions year × subarea (Pseudo-F = 1.3572, p(perm) = 0.0066), year × body-size class (Pseudo-F = 1.2591, p(perm) = 0.0262), and year × subarea x body-size class (Pseudo-F = 1.3523, p(perm) = 0.0004).

3.5. Importance of Prey Taxa in Hake Diet

Crustacean taxa were important prey for body-size class 0 individuals (%IRI = 77.64). Their importance was reduced for classes 1–4 individuals (%IRI = 5.57–19.29), whereas they were not present in the diet of the individuals of class 5. The %IRI of crustacean prey across body-size classes was 11.7. Crustaceans other than Decapoda were important prey essentially only for class 0 individuals (%IRI = 30.91). Notably, Parapenaeus longirostris emerged as the most important decapod prey for hake, with a maximum %IRI of 13.45 corresponding to the class 4 individuals. Its %IRI value across body-size classes was 2.72 (Table 2, Figure 3).

The relative importance of Cephalopoda in the diet of hake peaked in body-size class 5 individuals (%IRI = 36.44), while cephalopods were also present in the diet of class 2–4 individuals (%IRI = 0.04–0.76). The %IRI of cephalopods across body-size classes was 0.91. Illex coindetii stood out as the most important cephalopod prey for hake, registering an %IRI of 34.45 in class 5 individuals and 1.28 across body-size classes (Table 2, Figure 3).

Teleostei was by far the most important prey group in the diet of hake, with an %IRI of 87.4. Teleostei were more important prey for the individuals of the intermediate size classes, where their %IRI ranged from 79.95 (class 3) to 93.68 (class 4). Their importance was comparatively less for class 5 (%IRI = 63.56) and especially for class 1 individuals (%IRI = 22.36) (Table 2, Figure 3).

Engraulis encrasicolus stood out for its importance in the diet of hake, with an overall %IRI of 32.18. Its %IRI for the different size classes ranged from 8.6 (class 5) to 54.68 (class 1), while it was absent only from the diet of class 0 individuals. Unidentified Myctophidae were also important in the diet of individuals of the intermediate size classes, with their %IRI ranging from 0.32 (class 4) to 14.71 (class 2), while their %IRI across body-size class was 9.09. Other fish species that were among the main prey for hake included S. flexuosum, with an %IRI of 6.08 for class 3 individuals, M. poutassou, which was important for class 4 (%IRI = 14.86) and class 5 individuals (%IRI = 25.84), and S. pilchardus, which was of particular importance for class 4 individuals (%IRI = 30.38). Lastly, Trachurus trachurus and C. agassizi were among the main prey for class 5 individuals (%IRI = 14.05 and 8.26, respectively) (Table 2, Figure 3).

Across the years, the following persistent pattern regarding the relative importance of prey taxonomic groups was observed (Figure 3): Teleostei were the prey of the highest relative importance, the importance of Crustacea was intermediate, and the importance of Cephalopoda was the lowest. However, the relative importance of Teleostei and Crustacea varied between years, with the relative importance of Crustacea being higher in 2019, 2020, and 2023 (%IRI = 25.58, 24.37, and 31.87, respectively). Engraulis encrasicolus was the most important prey among fish species across the entire period (%IRI = 14.27–43.3), and Myctophidae was the second most important (%IRI = 6.24–14.58). The relative importance of S. flexuosum in 2021 was also notable (%IRI = 15.73).

A similar pattern in the relative importance of prey taxonomic groups was also observed across subareas (Figure 3): Teleostei and Crustacea were the most important prey groups in all subareas. Cephalopoda, represented mainly by I. coindetii, were less important prey for hake in THERM and THR-LIM (%IRI = 0.53 and 0.06, respectively) than in the other subareas (%IRI = 2.47–3.75). Engraulis encrasicolus was particularly important in THERM, THR-LIM, and VOL-SPOR (%IRI = 32.24, 54.85, and 52.82, respectively). However, its importance in CHIO-MIT and EVIA was greatly reduced (%IRI = 1.41 and 1.04, respectively), whereas the relative importance of Myctophidae there was the highest (%IRI = 20.87 and 18.56, respectively). The relative importance of C. agassizi and M. poutassou was also remarkable in CHI-MIT (%IRI = 5.23 and 2.61, respectively) and EVIA (%IRI = 1.23 and 4.55, respectively). Parapenaeus longirostris was relatively important in CHI-MIT, THR-LIM, and EVIA (%IRI = 2.43, 4.16, and 1.56, respectively). In the latter area, the relative importance of Plesionika shrimps was particularly notable (%IRI = 9.43).

In a force-based network graph layout algorithm like the “force-atlas” employed in the present study, linked nodes attract each other within the graph, while non-linked nodes repel each other. Indeed, prey with low out-degrees, i.e., linked to one to few hake body-size classes, were plotted at more peripheral locations in the network, whereas prey with high out-degrees, i.e., that were shared by several body-size classes, were plotted more centrally in the graph (Figure 4). Similarly, body-size classes 1–3 were clustered together because of their comparable trophic preferences. The strongest links were those between the prey taxa E. encrasicolus, unidentified Teleostei, unidentified Decapoda, and unidentified Crustacea, other than decapods and hake body-size classes 1, 2, and 3. All body-size classes, except class 0, exhibited connections with both exclusive (prey out-degree = 1) and shared prey taxa (prey out-degree > 1).

3.6. Trophic Niche Breadth

Prey numerical abundance in the stomach contents varied from 1 to 30, with most occurrences involving single individuals. Prey taxa with more than a few individuals in a stomach belonged to crustacean species other than decapods (up to 30 individuals), Decapoda (up to 15 individuals), and the Teleostei E. encrasicolus, S. pilchardus, and Myctophidae (up to 8 individuals). The richness of prey taxa in the stomachs ranged from 1 to 3 per stomach, with most stomachs containing a single prey taxon.

A statistically significant rank ordering of body-size classes in terms of prey diversity of order 0 was 2, 3 > 1, 4 (Figure 5, Supplementary Table S2). Also, prey diversity of order 0 in class 0 was significantly lower than all other classes except class 5, while the confidence intervals for class 5 curve overlapped also with those for class 4. In terms of prey diversity of orders 1 and 2, the rank ordering was 3, 4 > 1, 2 > 0, while the confidence intervals for class 5 overlapped with those for classes 2–4.

Regarding interannual differences in prey diversity, a statistically significant rank ordering of years in terms of prey diversity of all orders was 2020, 2021, 2022, and 2023 > 2019 (Figure 5, Supplementary Table S2).

A statistically significant rank ordering of subareas in terms of prey diversity of order 0 was not possible at all due to overlapping confidence intervals (Figure 5, Supplementary Table S2). For prey diversity of order 1, the only statistically significant difference was EVIA > VOL-SPOR.

Most of the prey taxa are depicted in the modified Costello graph situated in the upper-left section, having low frequencies of occurrence, but high prey-specific abundances. This pattern was observed in the plot for the pooled stomachs and in the plots for individual body-size classes, but to a lesser degree in class 5 (Figure 6).

3.7. Diet Overlap

Diet overlap, as quantified with the percentage overlap index (Table 3), was found to be higher between the intermediate body-size classes 1–4 (C_xy range = 0.45–0.69, mean = 0.56) and comparatively low between the extreme classes 0 and 5 and the intermediate classes 1–4 (C_xy range = 0–0.41, mean = 0.18). The index values exhibited a mean of 0.52 across years, with a range spanning from 0.36 to 0.63. The lowest diet overlap was observed between 2019 and the other years. The mean of the Renkonen index for diet overlap between subareas was 0.54, with values ranging from 0.40 to 0.70.

4. Discussion

The trophic preferences of M. merluccius were investigated in the North Aegean Sea during the summers of 2019 to 2023 with visual stomach contents analysis to assess its composition, diversity, and variability across body-size classes, years, and subareas. The principal prey were pelagic and benthopelagic organisms that are abundant in the area, and, in many cases, also primary targets for local fisheries. The identified prey are functionally diverse and inhabit a broad array of depths and habitats, underscoring hake’s importance as a pivotal component of the marine food web in the area. The ontogenetic trophic niche of hake was characterized by two distinct shifts: from a diet predominantly comprising crustaceans in individuals smaller than 10 cm, to a diet primarily including teleosts at intermediate body sizes, and ultimately transitioning to a diet based on teleosts and cephalopods in individuals larger than 50 cm. Cephalopods were a prevalent dietary component for large hake individuals, a finding that, to our knowledge, has not been reported in the literature. In comparison to the extreme body-size classes, the intermediate classes demonstrated greater trophic niche breadth concerning prey diversity and absolute prey-size ranges. Furthermore, a feeding strategy characterized by specialization of individuals was revealed, implying a dominant between-individuals component of the total trophic niche width. The effect of the temporal and spatial context on the hake diet was less pronounced than the ontogenetic changes and involved variability in the relative importance of prey groups and taxa that was occasionally correlated with spatiotemporal variations in the sizes of their populations.

4.1. Ontogenetic Shifts in the Diet of Hake

Numerous species undergo ontogenetic shifts in their food and habitat utilization, driven by variations in their vital rates, resource acquisition capacities, and predation risk relative to body size. Shifts in the species’ ontogenetic niche are thought to optimize their growth rates or surplus energy [64], while they are also considered to influence community and food web structure and functioning [65,66]. Ontogenetic dietary shifts are widespread among fishes [64,66], and several fish species have been used as models in studies of trophic ontogeny across all types of aquatic habitats [67,68,69,70]. The number of dietary shifts observed in fishes during their life histories varies considerably among species [64,66]. Similarly to hake in the study area, two dietary shifts have also been observed in Perca fluviatilis, in that case from planktivorous in juveniles to feeding on benthic invertebrates at intermediate sizes, and finally, when individuals become large enough, to a diet based mainly on fish [71]. Furthermore, profound multiple shifts in the diet occurring over the life cycle are frequently seen in piscivorous fish species [66,67]. The number of ontogenetic dietary shifts has also been found to vary between populations of the same fish species. For example, two dietary shifts were observed in the life history of Micropterus salmoides in Japan [72], whereas several shifts in the diet of the same species were observed in Spain [73].

In the present study, an abrupt shift from a diet based on crustaceans to one primarily based on teleosts was identified to occur in hake at a size threshold of 10 cm (TL). This shift was apparent under all approaches used to measure prey importance. Similar ontogenetic diet shifts have been regularly reported for hake in the Mediterranean and beyond, although the body sizes at which the shifts occurred often differed between studies. For example, in the area of the Ionian Sea, fish were found to gradually become the dominant prey of hake in terms of weight in individuals with sizes between 150 and 200 mm (TL) and in terms of frequency of occurrence in individuals between 200 and 300 mm long. In terms of numbers, crustaceans remained the dominant prey group across all body-size classes [25]. On the other hand, in the Tyrrhenian Sea, hake changed its diet from crustacean-based to one based on decapods and teleosts at a size of 10.5 cm (TL) [19], which coincides with the threshold identified in the present study. Beyond the Mediterranean, hake was found to be completely piscivorous in the Bay of Biscay and the Celtic Sea in individuals over 25 cm long, while individuals <20 cm in length fed primarily on crustaceans, mainly euphausiids [13,74].

A further notable change in the diet of hake that was observed in the present study concerned the significant increase in the relative importance and percentage, both by numbers and weight, of cephalopods, primarily I. coindetii, among hake individuals larger than 50 cm. Although cephalopods were part of the diet of hake also in body-size classes 2–4, their importance in those classes was comparatively low. Body-size class 5 individuals that fed on cephalopods were collected at sampling stations over the continental slope, with a mean depth of 410 m. Similarly, the body-size classes 2–4 individuals that preyed on cephalopods were collected at sampling stations with a mean depth of 271 m. Large hake individuals are capable of capturing prey at the size of cephalopods, which they may prefer as a food source because of their higher energetic values and digestibility compared to fish [75]. Additionally, cephalopods are more readily available in deep waters, where large hake are predominantly distributed [76]. For example, I. coindetii, an important cephalopod prey for hake in the study area, is known to be more abundant in deeper bathymetric zones [77].

Cephalopods generally comprise important food resources for fish, and generalist predators like various species of hake have occasionally been shown to consume large quantities of them [78]. Indeed, small cephalopods have been found to be important prey for small hake individuals in the Central Mediterranean, thought to reflect the local abundances of these species [79], while hake fed almost exclusively on the locally abundant cephalopods in Pomo Pit, in the Adriatic [22]. However, no cephalopods were found in the stomach contents of large hake individuals that exceeded 50 cm in length collected in a previous study in Greece [28], while the presence of cephalopods in the diet of hake collected in Saronikos Gulf was considered incidental [24]. The significance of cephalopods in hake diet was minimal also in the Egyptian Mediterranean waters [18], the Tyrrhenian Sea [19], and the Bay of Biscay, whereas they were not found among hake’s prey in the Celtic Sea [13]. Nevertheless, the transition in large individuals to a diet almost equally reliant on teleosts and cephalopods as was observed in the present study, to our knowledge, has not been previously documented for European hake, although it has been observed in the congeneric Merluccius paradoxus in the Southeastern Atlantic [75].

Diet overlap was low between the extreme and the intermediate body-size classes, while diet overlap between the latter was much larger. The findings of the present study support the delineation of three distinct phases in the ontogenetic trophic niche of hake in the North Aegean Sea, delimited by 10 and 50 cm (TL) body-size thresholds. However, it is important to note that these phases do not appear to align with specific life stages of the species. Specifically, the <10 cm phase corresponds to immature individuals less than 1 year old, as young-of-the-year lengths were reported to range in the Aegean Sea from 10 to 16 cm (TL) [10]. Additionally, in Thermaikos Bay, which is located in the North Aegean Sea, M. merluccius was observed to mature at lengths of 350 mm and 410 mm (TL) for males and females, respectively [10]. Consequently, the 10–50 cm phase corresponds to both immature and mature individuals, while the >50 cm phase includes only mature individuals.

Euphausiacea and/or Mysida are relatively small-sized crustaceans that have been shown to be important prey for hake in several studies across its distribution range [10,13,19]. It can thus be safely assumed that the unidentified small-sized crustaceans that were important for the diet of body-size class 0 individuals included mainly Euphausiacea and/or Mysida. Despite larger hake individuals primarily feeding upon larger, energetically richer teleosts, decapods, and cephalopods, small-sized crustaceans were also present in their diets, except in individuals larger than 50 cm. Feeding upon both small- and large-sized prey in the intermediate size classes indicates a widening of the range of absolute prey-size ratios with increasing predator size in hake, a finding that is in agreement with the results of previous studies for fish predators [80]. This pattern is considered to reflect improvements with increasing body size in predator behavioral and morphological capacities for capturing and ingesting larger prey, as well as increased encounter rates with smaller prey, which are also more susceptible to capture [80].

4.2. Spatial Variability in the Diet of Hake

While teleosts and crustaceans constitute the primary prey groups for hake across its range, their proportions and species composition in the diet of hake have often exhibited significant spatial variability, usually attributed to differences in their local abundances [10]. Indeed, the major hake prey in the study area, i.e., E. encrasicolus and Myctophidae, followed by P. longirostris and I. coindetii, are species whose populations are abundant in the North Aegean Sea ecosystem [35]. In particular, 75% of European anchovy production in Greece in 2022 was caught in the Northern Aegean [36]. Engraulis encrasicolus and Myctophidae were found to be significant prey for hake also in Western Greece [25], as well as in the Tyrrhenian Sea [19]. However, in Egyptian Mediterranean waters, Sardinella aurita and Caranx rhonchus were identified as more important prey for hake than E. encrasicolus [18]. Likewise, in the Celtic Sea, E. encrasicolus was less significant compared to T. trachurus or small hake [13].

Nevertheless, large-scale spatial variability in hake’s diet has not been attributed solely to differences in prey availability. In the Celtic Sea, variations in the proportions of main prey such as M. poutassou, T. trachurus, Scomber scombrus, and clupeoids across the different subareas were considered an indicator of hake “appetite” rather than reflecting spatial variability in prey abundances [74]. A more formal explanation is offered by the optimal foraging theory [81], which posits that predators strategically select prey that offer the highest energetic return relative to the effort expended in capturing and consuming them, thereby maximizing their fitness. Moreover, fish foraging behavior may also be influenced by past experience, the presence of competitors, or predation risk [64,82].

The influence of prey availability on the composition of the diet of hake in the study area was also evident when considering its variability between subareas. Indeed, the relative importance of E. encrasicolus as a prey for hake was diminished in CHIO-MIT and EVIA, where its biomass density was lower compared to the other subareas. Similarly, the increased relative importance of C. agassizi in CHI-MIT and EVIA, as well as of M. poutassou in EVIA, could be attributed to the higher biomass densities of these species in those areas. Likewise, P. longirostris was relatively more important as a prey for hake in subareas with a higher supply of this species, namely, CHI-MIT, THR-LIM, and EVIA. However, spatial variability in the relative importance of other prey in the diet of hake, like Cephalopoda, Myctophidae, and Plesionika spp., did not appear to be linked to the spatial variability of their biomass densities (Figure 7 and Supplementary Figure S7).

4.3. Interannual Variability in the Diet of Hake

The results of the PERMANOVA that was performed on the multivariate diet composition data revealed that the effect of the year was less important than the effect of body size, becoming significant only when diet composition was measured by weight. However, the effect of body-size class on hake diet composition was found to depend on the year, regardless of how prey proportions were calculated. This is also illustrated in the interannual variability in the %N and %W indices, as well as in the %IRI index, across body-size classes (Supplementary Figures S4–S6). Indeed, although the importance of Crustacea and Cephalopoda as prey for body-size classes 0 and 5, respectively, was evident across all years, the respective indices’ values differed considerably between years. However, this may be a biased result, as sample sizes were low for body-size classes 0 and 5. On the other hand, Teleostei was the most important prey group for all intermediate body-size classes across years, with few exceptions that occurred mostly when prey proportions were calculated in terms of numbers.

Studies of the trophic preferences of hake that are based on time-series data are scarce [20]. The results of the present study revealed that hake diet overlapped between years only to an intermediate extent, indicating that temporal factors had a less pronounced effect than body size/ontogenesis. The interannual variations in the diet of hake involved changes in the relative importance of Teleostei and Crustacea. At the species level, the prevalence of S. flexuosum in the diet of hake in 2021, rivaling the importance of hake’s main fish prey, E. encrasicolus, for that year, is particularly noteworthy. This finding can be attributed to the temporal biomass fluctuations of these species in the study area, as they experienced their respective minimum and maximum biomass densities for the study period in 2021 (Figure 7 and Supplementary Figure S7). The reduced abundance of E. encrasicolus in the North Aegean Sea in 2021 is possibly linked to the exceptional plankton bloom that occurred in the area during the spring and summer of that year. Furthermore, although P. longirostris was noted to be absent from hake stomach contents in previous analyses conducted in Greek seas, despite this species consistently constituting a significant proportion of Greek trawl catches [10], we identified P. longirostris as the primary decapod prey for hake in the study area during the study period across all body-size classes except class 5.

4.4. Cannibalism

Cannibalism was essentially not detected in the diet of hake in the study area during the sampling period (summer), as it was observed in only one individual. However, hake cannibalism has been observed in Greek seas [25]. The limited cannibalism observed in hake in the Celtic Sea [65] was attributed to spatial segregation of hake according to size. Cannibalism on young hake was also found to be important in the Mediterranean Iberian coasts [15], possibly related to the high abundance of hake recruits in the sampling period. However, it has been reported that cannibalism was absent in the populations of hake in the Tyrrhenian Sea [19]. The degree of cannibalism in hake has been associated with the presence of spatial overlap of predator and prey individuals and/or the availability of other prey [83]. In the study area, the MEDITS hauls during 2019–2023 that captured hake individuals of body-size class 0 also contained individuals of all other classes (Supplementary Figure S8). However, as hake individuals <10 cm long were collected mostly at sampling stations <350 m deep (Supplementary Figure S8) and hake body size was found to increase with depth (Supplementary Figure S9), the hypothesis of spatial segregation of hake in the study area according to size cannot be refuted.

Nevertheless, the availability of prey other than conspecific recruits has probably played an important role in the observed absence of cannibalism in hake in the study area. Less than one-third of the stomachs collected were empty, suggesting that, on average, food availability was not a significant limiting factor for hake in the North Aegean Sea during the study period. Furthermore, repletion was notably higher in individuals of body-size class 5 compared to most other classes. Similar substantial increases in food consumption among hake individuals longer than 30 cm have been reported in previous studies on the Greek seas [28,29]. Increased repletion in large hake individuals could be attributed to their capacity to prey upon fish and cephalopods of greater body sizes.

4.5. Trophic Niche Breadth

As sample sizes differed considerably between body-size classes, years, and subareas, the computation of diversity accumulation curves was considered necessary to estimate prey diversity. In terms of prey taxa richness (order 0 diversity), hake individuals of intermediate-size classes were found to have significantly more diverse diets. A similar pattern was observed when prey abundances in the stomach contents were also considered (orders 1 and 2 diversities). On the other hand, the lower prey diversity in size class 5 suggests enhanced dietary specialization in that class. As hake individuals > 50 cm long consumed only fish and cephalopods, it is assumed that ratio-based prey-size range in the particular body-size class was also reduced. This result is in agreement with previous studies [80] that found that ratio-based trophic-niche breadths in marine fish predators tended to narrow during ontogenesis for the largest predators (>50 cm long).

4.6. Feeding Strategy

Regarding the feeding strategy of hake and its variations between body-size classes, most prey taxa were depicted in the modified Costello graphs in the upper-left quadrant, i.e., were present in a few stomachs but with medium-to-high prey-specific abundances. This pattern, which was identified in both the plot for the pooled stomachs and the plots for the individual body-size classes, suggests that hake exhibited significant between-individuals variation (specialization of individuals) in its diet in the study area [59]. Most of the dietary niche width of hake was thus attributable to the between-individuals component [57]. However, prey-specific abundances of prey in body-size class 5 individuals were comparatively decreased, indicating reduced specialization of individuals, while the unsatisfactory taxonomic resolution of crustacean prey in body-size class 0 did not allow for a robust assessment of feeding strategy in the particular class.

4.7. Study Shortcomings and Future Directions

Seasonal variability in the diet of hake has been observed in the Greek seas and is considered to be related to the seasonal availability of its prey [10], with exceptions; for example, little seasonal changes in hake diet were found in the area of the Ionian Sea [25]. The reliance on samples collected only during summer is an important limitation of the present study, thereby not allowing the assessment of potential seasonal variability in hake diet in the study area. Additionally, the unsatisfactory taxonomic resolution attained for small-sized crustacean prey did not permit an unbiased estimation of trophic niche breadth and feeding strategy for individuals of body-size class 0. Nevertheless, limited taxonomic resolution in species diets is considered a characteristic shortcoming of visual analysis of stomach contents that modern DNA metabarcoding techniques are promising to overcome [84]. Furthermore, a sample size for body-size class 5 larger than that obtained in the present study would permit more robust conclusions to be made regarding the relative importance of cephalopods as prey for large hake in the study area. However, obtaining a large sample size for hake individuals larger than 50 cm is challenging, as they are typically collected from deeper waters, and their stomachs often become everted when the trawl net is lifted onboard. Lastly, although a considerable number of hake stomachs containing prey items were collected in the present study, the prey accumulation curves indicated that larger sample sizes would still reveal more prey taxa. This may be due to the high spatial resolution and multiannual temporal scale of the present study, while it also likely reflects the opportunistic nature of hake predatory behavior [74,85]. However, similar issues were also encountered in other visual studies of the trophic preferences of hake [21].

Important directions that future research could follow regarding the trophic preferences of hake comprise the investigation of seasonal variability and the evaluation of similarities and differences between different areas. Additionally, a comparative assessment of the diet of the species via alternative methods, e.g., with visual stomach content analysis and DNA barcoding, would provide a more comprehensive and accurate picture of its trophic interactions and its role as a component of the food web. Finally, as environmental variables have been shown by previous studies to influence the feeding and diet of hake [15], an investigation of the effects of the abiotic environment, for example, temperature, on the variability of the diet of hake in the study area would also be important. Indeed, variability in seawater temperature has been shown to affect the growth rates and the metabolic demands of predators, as well as the behavior of prey species, e.g., their movement and activity rates, which have an impact on encounter rates and probability of capture success [86,87]. In addition, spatiotemporal variability in temperature, for example, the sea temperature rise due to global climate change, has been shown to drive the distributions of many marine species, thus affecting predator–prey dynamics [88,89,90].

Assessments of food web structure and functioning offer a broad perspective on ecosystem dynamics by considering the trophic interactions between multiple species. Such approaches complement single-species assessments, contributing to a comprehensive understanding of ecosystems that is essential in EBFM [91,92]. Indeed, quantification of species trophic interactions and food web structure can inform the formulation of promising new harvesting strategies for highly exploited multispecies fisheries, such as balanced harvesting [93,94]. However, gathering comprehensive data on marine species diets is challenging, as omnivory is widespread in marine food webs [95], trophic preferences typically change during species ontogenesis [66], and spatiotemporal changes in prey availability often affect diet composition in opportunistic predators [96]. The findings of the present study underscore that DCF represents a unique framework to obtain a time series of ontogenetic stage-specific and spatially comprehensive data regarding the trophic interactions of key biological resources, following a systematic and comprehensive approach.

5. Conclusions

This study provides comprehensive recent information on the ontogenetic and spatiotemporal variability of the diet of European hake in the North Aegean Sea in summer, thereby contributing to filling important knowledge gaps regarding a major biological resource and key component of the marine food web in the Mediterranean. European hake was found to be a generalist predator in the study area, with its wide trophic niche resulting from a combination of individual specialization and a high variance in resource use between individuals. Ontogenesis is accompanied by two major changes in the trophic preferences of the species early and late in its life history, while spatiotemporal variability in its diet was occasionally correlated with spatiotemporal variations in the sizes of prey populations. European hake is a key node in the food web of the Mediterranean and an overexploited major fisheries target. Data on its trophic interactions, such as the information presented in this study, are critical for the ecosystem-based approach to fisheries management to ensure the sustainable use of its populations.