**Biologic Aspects of Juveniles of the Common Stingray,** *Dasyatis pastinaca* **(Linnaeus, 1758) (Elasmobranchii, Dasyatidae), from the Central Mediterranean Sea**

#### **Francesco Tiralongo 1,2,3,\* , Giuseppina Messina <sup>1</sup> and Bianca Maria Lombardo <sup>1</sup>**


Received: 24 February 2020; Accepted: 7 April 2020; Published: 10 April 2020

**Abstract:** Data on the biology of *Dasyatis pastinaca* are absent from the Ionian Sea and only a few studies were conducted in the Mediterranean Sea. Some biological and ecological aspects of *D. pastinaca* were investigated between November 2019 and February 2020 in the central Mediterranean Sea. In particular, we investigated several morphologic, population and ecological aspects of the species. The analysis of the stomach contents showed that *D. pastinaca* is a generalist carnivorous, mainly feeding on small crustaceans and polychaetes. The Levin's index value (B*i*) was 0.85. The sex ratio showed no significant differences from 1:1 ratio. Females were larger than males, but no statistical differences were found in disc width-weight and total length-disc width relationships between sexes. Most of the specimens caught were juveniles and inhabit shallow sandy bottoms.

**Keywords:** eastern Sicily; Batoidea; elasmobranchs; diet; coastal fishery

#### **1. Introduction**

Elasmobranchs are key top predators in most marine environments and play an essential role in regulating and structuring marine ecosystems [1,2]. On the other hand, due to their low fecundity and delayed age at maturity, elasmobranchs are highly vulnerable to fishing activity and are often affected by high by-catch rates [3,4]. All this is reflected in the dramatic decline of shark and ray populations and, as demonstrated by several researches conducted in the Mediterranean Sea over the last decades, most species of elasmobranchs have dramatically declined in number mainly due to overfishing and illegal fishing [5–9].

In Italian seas, stingrays (Dasyatidae) are represented by 4 species and 3 genera [10]: *Dasyatis* Rafinesque, 1810, *Pteroplatytrygon* Fowler, 1910 and *Taeniura* Müller & Henle, 1837. The genus *Dasyatis* comprises *Dasyatis pastinaca* (Linnaeus, 1758) and *Dasyatis centroura* (Mitchill, 1815). These two species are very similar but can be distinguished by the presence (*D. pastinaca*) or absence (*D. centroura*) of a dorsal keel behind the spine on the tail and by the presence of quite developed spines and tubercles on the dorsal surface of the body and tail in large specimens of *D. centroura* [11]. Furthermore, between all the Mediterranean species of the genus *Dasyatis*, *D. pastinaca* is the most abundant and widely distributed one [12]. *Pteroplatytrygon violacea* (Bonaparte, 1832) is a less common species than *D. pastinaca*, while *Taeniura grabata* (Geoffroy Saint-Hilaire, 1817) is very rare and was only occasionally caught in the Mediterranean, and only once in Italian waters [10,13]. Because of their low or no commercial value, stingrays are usually discarded [14–16]. However, due to their vulnerability to fishing activities, the conservation and correct management of these vulnerable K-selected species require special efforts to ensure their perpetuation.

The common stingray, *D. pastinaca*, is an Atlantic-Mediterranean species whose distribution extends from the southwestern Baltic Sea and throughout the Mediterranean Sea to Senegal [17]. It is listed as "data deficient" in the IUCN Red List [18]. The common stingray inhabits soft bottoms of coastal waters, down to a depth of about 200 m [19,20]. It is commonly caught with trammel nets and trawls and usually discarded [14–16]. The maximum recorded disc width (DW) is 64 cm, but common sizes range from 20 to 40 cm DW. The survival rate of this species, when caught with trammel net, is equal or close to 0 [16]. Indeed, because of their defensive behavior, fishermen are forced to kill the stingrays with the aim of avoiding serious injuries from their venomous spine on the tail. Another noteworthy fact of this species is the antimicrobial and anti-proliferative effects of its skin mucus [21].

Data on the biology and ecology of *D. pastinaca* in the Mediterranean Sea are scarce, and no data are available from the Ionian Sea. Therefore, the aim of this study is to provide new data on the biology and ecology of the common stingray, providing the first data from specimens from the Ionian Sea. In particular, we investigated size frequency distribution, sex ratio, disc width-weight relationships, total length-disc width relationships, diet composition and habitat of *D. pastinaca* from specimens caught as discards in trammel nets targeting cuttlefish, *Sepia o*ffi*cinalis* Linnaeus, 1758, in the central Mediterranean Sea.

#### **2. Materials and Methods**

A total of 120 specimens of *D. pastinaca* were collected along the coastline extending for about 21 km in the southeast coast of Sicily (Ionian Sea), from Avola to Marzamemi (Figure 1). The specimens were caught by professional fishermen at 5–30 m depth between the 10 November 2019 and 10 February 2020 with trammel nets targeting cuttlefish (*S. o*ffi*cinalis*) [16]. Trammel nets were deployed overnight (from 6 pm to 4 am) for about 10 h, on sandy and mixed bottoms (sand and rocks), close to *Posidonia oceanica* seagrass meadows. In order to better represent the population, specimens were randomly selected from different fishing vessels operating in the area. Data about bottom nature and depth of capture were collected through interviews to fishermen.

**Figure 1.** Study area (indicated in red) in the Ionian coast of Sicily (central Mediterranean Sea).

Each specimen was weighed and measured (disc width and total length). The sex was determined by the presence (male) or absence (female) of claspers. Weight and disc width measures were used for the disc width-weight relationships following the formula: W = aDWb, where W is the weight in grams (g), DW is the disc width in centimeters (cm), a is the intercept and b is the slope of the regression curve. Total length (TL) and disc width measures were used for the total length-disc width relationships following the formula: DW = yTL + x, where y is the slope and x is the intercept of the regression line. Disc-width frequency distributions were constructed for both sexes. A chi-square test was used to verify if there was a significant difference (α = 0.05) between the observed and the expected sex ratio (M:F, 1:1) of the whole sample. To test if the regressions of the weight (W) on disc width (DW) were significantly different (α = 0.05) for the two sexes, and to test if the regressions of the disc width on total length (TL) were significantly different (α = 0.05) for the two sexes, an analysis of covariance (ANCOVA) was employed. Five outliers (females > 2143 g) were excluded from the analyses.

The stomach was removed from each fish as soon as possible after landing, and its content analyzed. All the prey items in the stomachs were counted, washed in clean seawater and dried with blotter paper, identified under a microscope to the lowest taxonomic level possible and weighed to the nearest 0.01 g.

The frequency of occurrence (%F), percentage weight (%W), percentage abundance (%N) and the Index of Relative Importance (%IRI) were calculated for each prey category [22,23]. The vacuity index (percentage of empty stomachs) was also calculated.

According to the value of their percentage abundance (%N), prey were grouped into three categories [24]: dominant (N > 50%), secondary (10% < N < 50%) and accidental (N < 10%).

The feeding strategy of *D. pastinaca* was visually examined using a modified version of the Costello (1990) graph [25] by plotting the prey-specific biomass (P*i*) against their frequency of occurrence (%F) [26]:

$$\text{Pi} = \frac{\text{SWi}}{\text{SWt}\_{i}} \times 100\tag{1}$$

where P*i* is the prey-specific biomass of prey *i*, SW*<sup>i</sup>* the stomach content biomass of prey *i*, and SW*ti* the total stomach content biomass in those predators with prey *i* in the stomach.

Standardized Levin's index (B*i*) was used to evaluate the breadth of the diet [27]:

$$\mathbf{B} = \frac{1}{\sum p\_j^2} \tag{2}$$

$$\text{Bi} = \frac{\text{B} - 1}{\text{B}\_{\text{max}} - 1} \tag{3}$$

where *pj* is the relative frequency specimens in the jth prey item and B*max* is the total number of prey item categories found. B*i* is comprised between 0 and 1. The more the value of B*i* is close to 0, the narrower is the trophic niche of the species investigated. Conversely, the closer the value of B*i* is to 1, the wider the trophic niche. Hence, if B*i* < 0.40, the species is considered a "specialist feeder"; if 0.40 < B*i* < 0.60, the species is considered an "intermediate feeder"; if B*i* > 0.60 the species is considered a "generalist feeder" [28].

A cumulative prey curve [29] was computed with R Studio [30] in order to evaluate whether the number of analyzed stomachs was sufficient to describe the diet of the species. The estimated number of prey groups with the associated SD were plotted against the cumulative number of individuals whose stomach was examined.

#### **3. Results**

Out of the total 120 specimens sampled, 57 (47.5%) were males and 63 (52.5%) were females. The sex ratio did not deviate significantly from 1:1 (M:F, 0.9:1; *p*-value > 0.05). Females were on average larger than males, ranging from 31 to 76.4 cm in total length, from 16.8 to 45.5 cm in disc width and

from 144 to 3781 g in weight (Figure 2, Table 1). However, the ANCOVA test did not indicate significant differences (*p*-value > 0.05) in disc width-weight relationships between sexes and neither did total length-disc width relationships (*p*-value > 0.05). The disc width-weight relationships for both sexes and combined showed a positive allometry (Table 2, Figure 3). The parameters of the linear regression of the total length-disc width relationships were: DW = 0.641TL−4.2508 for males; DW = 0.641TL−4.122 for females; and DW = 0.641TL−4.1867 for combined sexes (Figure 4).

**Figure 2.** Sex composition by size classes of *Dasyatis pastinaca*.

**Table 1.** Sex distribution, total length, disc width and weight of *Dasyatis pastinaca* in the Ionian Sea (central Mediterranean Sea). The sex ratio (M:F, 0.9:1) was not significantly different from 1:1 (chi-square: *p* > 0.05).


**Table 2.** Disc width-weight relationships parameters of *D. pastinaca* in the Ionian Sea (central Mediterranean Sea); C.I. = 95% confidence interval.


**Figure 3.** Disc width-weight relationships of *D. pastinaca* for both sexes (**a**) and combined (**b**). A: in red for females (N = 63) and blue for males (N = 57).

**Figure 4.** Disc width-length relationships of *D. pastinaca* for both sexes (**a**) and combined (**b**). A: in red for females (N = 63) and blue for males (N = 57).

Considering the size at first maturity (L50 = 62.5 cm and 43 cm in DW for males and females, respectively) reported by a study conducted in the nearby Aegean Sea [31], a high percentage of specimens (virtually all males and the majority of females) were juveniles. All the specimens were caught on sandy bottom and within the whole bathymetric range investigated (5–30 m). A total of only 5 specimens (all females > 2143 g) were of relatively large size.

The cumulative prey curve showed that relatively few stomachs (about 50) are necessary for a reliable description of the diet of *D. pastinaca* (Figure 5). Out of the 120 specimens analyzed, only 67 had prey in their stomach (vacuity index of 55.83%). The stomach content analysis revealed that juveniles of *D. pastinaca* mainly feed on small benthic crustaceans and also, to a lesser extent, polychaetes (Table 3). Fish represented a negligible part of the diet. However, no preference for a particular prey was observed: the value of the standardized Levin's index (B*i*) was 0.85 (the digested category was included in the analysis), indicating a wide trophic niche. Most of the prey types were included into the category of "accidental prey" (N < 10%) (Table 3) and only a few (Amphipoda, Mysida and Lumbrineridae) were classified as "secondary prey" (10% < N < 50%). No dominant (N > 50%) prey types were found between the 13 identified items. Most of the prey were represented by small benthic invertebrates inhabiting sandy bottoms (Table 3). The values obtained for %IRI, %F and %W also indicated no clear dominance of any prey (Table 3). The same was clearly indicated by the plotted results of the main prey-specific biomass (P*i*) against the frequency of occurrence (%F), being in agreement with results of the other previous analysis (Figure 6).

**Figure 5.** Cumulative prey curve (in red) as a function of sample size for all stomachs analyzed of *D. pastinaca*. Standard deviation (SD) represented with blue lines.

**Table 3.** Diet composition of *D. pastinaca* from the Ionian Sea (central Mediterranean Sea). %F = percentage frequency of occurrence; %N = percentage in number; %W = percentage in biomass; IRI = Index of Relative Importance of prey items and its percentage (%IRI). Values > 10% are in bold.


**Figure 6.** Prey-specific biomass (P*i*) plotted against the frequency of occurrence (%F) of the main prey items for *D. pastinaca* from the Ionian Sea (central Mediterranean Sea).

#### **4. Discussion**

With the exception of a few large female specimens, all or most of the individuals were juveniles. The sex ratio (statistically not significantly different from 1:1) of the fish caught during the study period indicates that there is no sexual segregation in juveniles of *D. pastinaca*.

The disc width-weight relationships showed a positive allometric growth in both sexes. Population structure and growth parameters, such as length-weight relationships are valuable measurements in order to evaluate the health conditions of fish populations [32] and are useful for comparisons of life histories of species among regions [33]. Also, total length-disc width relationships are valuable measurements to investigate intra and interspecific variation among populations [34]. Females showed larger sizes than males, but no significant differences were found in disc width-weight and total length-disc width relationships between the two sexes. Males and females attained smaller sizes than those of specimens from the Aegean Sea and Cilician Basin [31,35]. On the contrary, they reached similar size to those from the Iskenderun Bay [36]. Size differences appears to be related to the sampling depth. Indeed, while our specimens and those from the Iskenderun Bay were caught in shallow waters (5–30 and 0–50 m, respectively), the specimens from the other two Mediterranean areas (Aegean Sea and Cilician Basin) were collected in deeper waters (>50 m, down to a depth of 500 m). In consideration of this, we can speculate that juveniles of *D. pastinaca* use the study area, consisting of sandy bottoms of shallow waters close to *Posidonia oceanica* meadows, as a feeding ground. A similar result was obtained for juveniles of the rough ray, *Raja radula* Delaroche, 1809 [37].

On the basis of the number of specimens collected during the study period—and especially from interviews with local fishermen—*D. pastinaca* can be considered common in shallow waters of the Ionian coast of Sicily (central Mediterranean Sea). We also point out that in the investigated area *D. pastinaca* represented about 19% of the total elasmobranch catches in trammel nets targeting *S. o*ffi*cinalis* [16]. In the southeastern Sicily, the species is sympatric with other common species of elasmobranchs: *R. radula* and *Torpedo torpedo* (Linnaeus, 1758). However, while fishermen caught both adults and juveniles of *T. torpedo* [38], in *R. radula*, as in the case of *D. pastinaca*, most of the specimens were juveniles [37]. Furthermore, considering the different feeding habits between *D. pastinaca* and *T. torpedo*, their trophic niche does not overlap, resulting in little or no direct competition between these two species. Indeed, while *D. pastinaca* mainly feeds on small benthic crustaceans, as also demonstrated in this study, *T. torpedo* is a generalist piscivore [38,39]. The different diets of the two species (*D. pastinaca* and *T. torpedo*) can be explained by a different predation strategy. Indeed, only *T. torpedo* is capable of producing electric discharges to stun prey before swallowing them, and therefore this ambush predator feeds on more mobile and larger prey such as fish. Conversely, between *D. pastinaca* and *R. radula* the diet overlap was consistent. Both species mainly feed on small benthic crustaceans, with *R. radula* being a more selective feeder [37]. Hence, the differences in the diets of these species and, consequently, the different trophic competition between them, could be the cause, at least in part, of the different abundances of the species in the area: ~19%, ~33% and ~47% of the total catch of elasmobranchs, respectively for *D. pastinaca*, *R. radula* and *T. torpedo* [16].

Crustaceans (total percentage in number > 75%), followed by polychaetes, were the dominant prey. Considering the IRI, the main prey among crustaceans were represented by mysid shrimps, *Liocarcinus* sp. and Caridea. Other studies from other Mediterranean areas also showed similar results. In particular, Yeldan et al. (2009) [35] showed that the diet of specimens from Cilician coasts (Turkey) was mainly composed by crustaceans. Another study from the Black Sea (Turkey) [40] showed a similar result and similar results were also reported from the northeastern Atlantic (Azores) and from the eastern Adriatic Sea (Croatia) [41,42]. Furthermore, although the crustacean species composition varies between our and the above mentioned studies, some decapods species seem to be quite common in the diet of *D. pastinaca* and were reported from two or more of the above mentioned studies (including our study): *Alpheus glaber* (Olivi, 1792), *Crangon crangon* (Linnaeus, 1758), *Liocarcinus depurator* (Linnaeus, 1758) and *Upogebia pusilla* (Petagna, 1792). In all cases, the high value of the Levin's index (0.85) and the plot of the prey-specific biomass (P*i*) indicated that *D. pastinaca* is a generalist feeder.

In the marine environment, sharks and batoids often occupy the highest trophic levels of the food chain [43]. The key ecological role of these species is however threatened by human activities such as illegal fishing, overfishing and habitat destruction [1], which have in some cases led to dramatic shift in relative abundance, sometimes resulting in a great increase in smaller, lower-trophic level species. However, despite the importance of understanding the feeding relationships and the energy transfer in marine ecosystems, little is known about the feeding ecology of most elasmobranch species, and this is particularly true for batoids, which have received considerably less attention than sharks [44,45]. The importance of stingray predation in regulating the abundance and composition of prey in shallow waters has been demonstrated by several authors [46,47]. Hence, there is a need to protect these vulnerable and important key species. It therefore appears of great importance in management and conservation strategies such as, for example, the establishment of marine reserve, to have a good understanding of the trophic differences and relationships among species. Hence, we also need to deepen the knowledge about the trophic ecology of these species [48].

In conclusion, our study provides new data on the biology and ecology of *D. pastinaca* in the central Mediterranean Sea (where the only available data were from the eastern Adriatic Sea), and first data from the Ionian Sea. *Dasyatis pastinaca* juveniles were found common in sandy bottoms of the investigated area. Juveniles are quite active feeders, and mainly feed on small benthic invertebrates, with a general preference for decapod crustaceans and polychaetes. However, there was no clear preference of any particular species. Additional targeted studies are needed to provide new useful data for the management and protection of this potentially vulnerable k-selected species. Still little is known about several biologic and ecological aspects of *D. pastinaca*, such as population structure, abundance, distribution and habitat selection and no data are available from the western part of the Mediterranean Sea. The experience and knowledge of fishery workers is a valid method for data collection [49,50]; a close collaboration between scientists and fishermen can lead to a better understating of the ecology and distribution of elasmobranchs.

Finally, particular attention is needed concerning the impacts of coastal fisheries on this and other coastal elasmobranch species. Considering the relatively high percentage of juveniles of the species of batoids caught in the area, a reduction of fishing pressure on these resources should be encouraged. A successful strategy could be the release of juvenile and adult specimens that are still alive and in good conditions in nearby selected areas immediately after landing.

**Author Contributions:** Conceptualization, F.T.; methodology, F.T. and G.M.; software, F.T.; validation, B.M.L.; formal analysis, F.T. and G.M.; investigation, F.T. and G.M.; data curation, F.T., G.M. and B.M.L.; writing—original draft preparation, F.T.; writing—review and editing, F.T., G.M and B.M.L.; supervision, B..M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We are grateful to the fishermen of Avola and Marzamemi for their support and willingness to freely share the information they had and for their help with the collection of the specimens.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **Review of Estimating Trophic Relationships by Quantitative Fatty Acid Signature Analysis**

**Junbo Zhang 1,2,3,4,\*, Chonglan Ren 1,**†**, Hu Zhang 5, Fang Yin 6, Shuo Zhang 1,\*, Rong Wan 1,4,\* and Daisuke Kitazawa <sup>7</sup>**


Received: 26 October 2020; Accepted: 14 December 2020; Published: 18 December 2020

**Abstract:** The dynamic predator–prey relations in the food web are vital for understanding the function and structure of ecosystems. Dietary estimation is a research hotspot of quantitative ecology, providing key insights into predator–prey relationships. One of the most promising approaches is quantitative fatty acid signature analysis (QFASA), which is the first generation of statistical tools to estimate the quantitative trophic predator–prey relationships by comparing the fatty acid (FA) signatures among predators and their prey. QFASA has been continuously widely applied, refined and extended since its introduction. This article reviewed the research progress of QFASA from development and application. QFASA reflects the long-term diet of predator, and provides the quantitative dietary composition of predator, but it is sensitive to the metabolism of predator. The calibration coefficients (CCs) and the FA subset are two crucial parameters to explain the metabolism of predators, but the incorrect construction or improper use of CCs and the FA subset may cause bias in dietary estimation. Further study and refinement of the QFASA approach is needed to identify recommendations for which CCs and subsets of FA work best for different taxa and systems.

**Keywords:** quantitative fatty acid signature analysis; aquatic food webs; dietary estimation

#### **1. Introduction**

To understand the structure and function of aquatic ecosystems and realize the trophic predator–prey relationships, it is necessary to obtain accurate dietary composition of predators. In aquatic ecosystems, the impracticality and limitations of direct observation method of feeding have prompted the development of indirect dietary estimation methods [1], such as the stomach content analysis [2], the stable isotopes (SI) analysis [3] and fatty acid (FA) analysis [4]. As one of the promising methods, FA analysis is based on the fact that specific FA can only be produced by certain primary producers (usually algae and bacteria), and higher trophic levels cannot synthesize by themselves; therefore, it can be used to track food sources. Fatty acids (FAs) have been widely applied in qualitative

studies for analyzing the trophic predator–prey relationships. Then, Iverson et al. [5] developed a new method named quantitative fatty acid signature analysis (QFASA), which was the first time FAs have been used to quantitatively estimate the diet of predators.

The QFASA model has been continuously refined and extended, though it has been extensively applied in diet estimation. This article reviews the development of QFASA in the dietary estimation of aquatic species, and discusses the application scope of QFASA based on its merits and demerits. The aim of this paper is to provide a reference when applying QFASA to the quantitative evaluation of trophic relationship in marine ecosystem.

#### **2. Fundamental Requirements of QFASA**

QFASA is related to the characteristics of FAs. The differences in FA biosynthesis among species make it possible to identify the original source of certain FA. Firstly, FAs in prey are transferred, largely unaltered or at least in a predictable manner, to lipid storage sites (i.e., adipose tissue) of their predators, and reflect the long-term diet. Secondly, FAs derived from food are called dietary FAs, some of which can only be produced by specific species (such as algae and bacteria), and are essential for the survival of organisms, also called essential FAs [6]. Thirdly, in contrast to primary producers, other species usually can only produce some simple FAs (e.g., 14:0, 16:0, and 18:0 saturated FA; 14:1n-5, 16:1n-7, and 18:1n-9 monounsaturated FA) [7].

To understand the QFASA model, it is necessary to understand the terms of FA signature and calibration coefficient (CC).

The first issue is FA signature. According to Iverson et al. [5], the quantitative distribution of all FAs obtained in samples (prey or predators) is called FA signature. Each prey species has a unique pattern of FA signature due to the difference of life history and dietary habits [8]. As these signatures are conserved through the food chain, they can act as indicators of dietary composition [8].

The second issue is CC, which is a coefficient to calibrate the lipid metabolism process in the fat tissue of the predator [5]. Due to the biosynthesis, mobilization, modification and deposition of FAs, the FA signature of the predator is not completely consistent with its prey. Calibration coefficients (CCs) were determined previously from feeding experiments of individuals on a given long-term diet, with the assumption that after such long-term feeding, the adipose tissue stored in the predator will be similar to the given diet [5,9–11]. Then it can be calculated as the discrepancy or ratio, between adipose tissue and diet levels of each FA [5], as follows:

$$\text{CC} = \text{Predator/Prey} \tag{1}$$

Fundamental requirements for the application of QFASA are as follows: (i) a predator tissue properly sampled, stored, extracted and analyzed for FA signature; (ii) an appropriate sampling analysis of FA signature of potential prey species for a predator; (iii) the lipid metabolism and deposition of predator estimated, usually through captive feeding experiments; (iv) a quantitative model to minimize the distance of FA signatures between a predator and its prey [5]. More detailed information about the initial QFASA model can be found in the study of Iverson et al. [5], and further discussion of the model can be found in the article by Budge et al. [12].

QFASA is the first generation of statistical tools to evaluate the quantitative trophic predator–prey relationships by comparing the FA signatures among predators and their prey [5]. Taking into account the metabolism of the predator, the technique aims to find the combination of prey, based on their FA signatures, which yields a combined FA signature similar to that of the predator [5]. QFASA is developed based on two main hypotheses: the representative prey signature database contains all potential prey taxa and CCs are known and accurate [13].

QFASA can be used for individuals as well as food webs. Figure 1 is revised from Magnone et al. [14], who applied QFASA to generate a food web model in the aquatic environment to find out the trophic

relationships among the species. The arrows point from all predators to their prey and the percentage represents the proportion of prey in the diet of predator.

**Figure 1.** Quantitative trophic predator–prey relationships. (Revised from Magnone et al. [14]).

#### *Ecological Modeling*

The dietary assessment of each predator is acquired from the weighted combination of its prey's FA signatures, and then the weighted coefficient with the lowest statistical distance from the predator's FA signature is selected. This is achieved by minimizing the statistical distance between the prey's FA signatures weighted combination and predator's FA signature. The flow chart of QFASA application is shown in Figure 2 and the ecological modeling is as follows [5]:

The QFASA model considering the vectors of prey (*x*) and predator (*y*), originally developed by Iverson et al. [5], combined with the Kullback–Liebler (*KL*) distance, can be expressed as:

$$\mathfrak{F} = \sum\_{k} p\_k \mathfrak{k}\_k \tag{2}$$

$$KL = \sum\_{j} (y\_j - \hat{y}\_j) \log \left(\frac{y\_j}{\hat{y}\_j}\right) \tag{3}$$

where *yˆ* is weighted sum of the FA of the prey, *x*ˆ*<sup>k</sup>* is a vector of the prey species *k*, *x*ˆ*kj* represents the mean of each FA *j* of the prey species *k*, *pk* is a weight coefficient corresponding to the evaluated proportion of the *k* in the diet of predator, *yj* represents the proportion of each FA *j* of the predator, and *y*ˆ*<sup>j</sup>* represents the estimated value of each FA *j* of the predator.

The aim of the model is to select the *pk* that makes the estimated value *yˆ* as near as possible to the true value *y*. Because the *pk* values are greater than or equal to zero and the sum of *pk* values is 1, the evaluation can be transformed into a constrained optimization problem [14,15],

$$\text{minKL} = \sum\_{j} (y\_j - \hat{y}\_j) \log \left( \frac{y\_j}{\hat{y}\_j} \right) \tag{4}$$

$$\sum\_{k} p\_k = 1, \ p\_k \ge 0 \tag{5}$$

**Figure 2.** The application flow of quantitative fatty acid signature analysis (QFASA).

#### **3. The Research Progress of QFASA**

#### *3.1. The Improvement of QFASA*

The improvement of the QFASA model mainly includes the testing of the method in new taxa, selection and refinement of prey library, calculation of CC, selection of the FA subsets, and optimization of the statistical model.

#### 3.1.1. The Testing of the Method in New Taxa

At first, QFASA was intended for estimating the dietary composition of marine mammals [5]. Subsequently, it has been used in multiple predators in various ecosystems, including marine mammals [16–19], seabirds [9–11,20], and fishes [14,15,21–23]. The development history of QFASA is shown in Figure 3, and each node (indicated by a circle) in the figure represents the first application to this species type.

In previous research, as shown in Figure 3, Iverson et al. [5] developed QFASA to estimate the dietary composition of *Halichoerus. grypus* and *Phoca. groenlandica*, and concluded that QFASA can be applicable to various predators and ecosystems; Iverson et al. [9] also applied QFASA to seabirds for the first time, and estimated the diets of *Uria. aalge* and red-legged *kittiwakes*; Wang et al. [10] estimated the dietary composition of threatened species (*Somateria. Wscheri* and *Polysticta. stelleri*) by QFASA, and inferred that QFASA can be used with other wild eiders or birds, and can estimate their diet composition at different life stages; Budge et al. [21] applied QFASA to marine fish for the first time and studied the diets of Atlantic salmon (*Salmo. salar*) (Figure 3); in addition, Magnone et al. [14] used

QFASA to generate a food web model in the aquatic environment to find out the trophic relationships among the species, and this was the first time that QFASA has been used to construct food web.


**Figure 3.** Development history of QFASA's dietary items.

In recent studies, Goetsch et al. [19] estimated the dietary composition of *Mirounga angustirostris* across five years; Conners et al. [20] studied the dietary composition of two albatross species through QFASA, which was beneficial to understand their fishery exploitation; Knox et al. [24] also evaluated the dietary composition of *Arctocephalus pusillus doriferus* by QFASA. It is obviously shown that the predator taxa of QFASA are dominated by marine species, mainly including marine mammals, seabirds, and fishes. In addition, as shown in Figure 3, there have been several studies focused on the diet of freshwater fishes in recent years [22,23,25,26], such as *Salvelinus namaycush* by Happel et al. [22].

#### 3.1.2. The Selection and Refinement of Prey Library

A pre-condition of QFASA is a prey library which contains the potential prey species of predator [5], and each prey species has a unique FA signature. QFASA uses the average value of each prey to evaluate its contribution to the signature of the predator. However, wild animals consume individuals of different species rather than homogeneous species [5]. Although it is likely to identify the species in the ecosystem by its FA signature, there may be considerable differences within species, such as the same species but different ages or sizes. For instance, studies have found that the lipid content of sprat changes drastically with age and size [27]. Moreover, Magnone et al. [14] treated the young and adult *Paralichthys orbignyanus* as two types, and found that the adult *P. orbignyanus* was the top predator in the ecosystem studied while the young *P. orbignyanus* was prey of other species, and the dietary composition was different from each other. Therefore, species subgroups can be introduced into the QFASA model to provide extra details of diet. It is recommended that the age and size of prey should be taken into account when applying the QFASA model, and species with drastic changes in lipid content could be regarded as two types.

Besides, it is necessary to ensure that the quantity of FAs used in model is no less than the amount of potential prey in the database [19], otherwise, the estimated value of diet will be non-unique [18,28]. Three solutions are used to solve this problem based on previous studies. Firstly, multiple prey types with FAs alike could be assembled as a single type before applying QFASA [18,29,30]. Secondly, a post hoc method, assessing the diet of each prey as a unique type, is then pooling the estimated values of those prey into their respective communities [17]. Thirdly, Goetsch et al. [19] proposed an ad hoc approach (Drop Core Prey Analysis) to recognize the prey that did not contribute much to the diet and eliminate them from the model. In general, the first method may cause the mean FAs of the mixed prey species to be dissimilar to the practical prey, and the second method cannot ensure if the estimated value of the mixed diet is unique [19]. In contrast, the third method could avoid the deficiencies of the two methods above.

#### 3.1.3. The Calculation of CCs

In addition to a detailed knowledge of the FA signatures of all potential prey, an understanding of the biochemistry and metabolism of FA within predators is essential to use this method accurately. At present, QFASA has two crucial parameters to explain the metabolism of predators: the set of CC for individual FA and the subset of FA used in modeling diets [5].

The first issue is CC, which accounts for the lipid metabolism of the predator. CCs were determined previously from captive feeding experiments. Although the CCs for a variety of birds, mammals, and fishes [11,15] have been determined, which is shown in Table 1, feeding experiments of most predators have not been conducted. In that condition, the developed CC of a surrogate species is normally used based on the evolutionary and ecological similarities. For instance, Thiemann et al. [16] used the CC of *Mustela vison* to replace the CC of polar bears; Meynier et al. [17] compared five sets of CC calculated by Iverson et al. [5] and Tollit et al. [31] then chose available CC that can be used for the application of QFASA on *Phocarctos hookeri*; Haynes et al. [30] used three CCs obtained by Iverson et al. [9] and Wang et al. [10] to replace the CC of *Gavia adamsii* by feeding trials.

Some studies have shown that diet estimation is very sensitive to the selection of CC, and the errors in CC would cause bias in dietary estimation [13,30]. Whereas, the factors such as feeding time and the type of food fed may affect the adaptability of the obtained CC in the study of wild predators [21,32] and the accuracy of CC in dietary estimation of wild animals cannot be verified. These issues have caused considerable criticism, and some recent studies emphasized and pointed out that using QFASA without CC is also sufficiently robust [15,21,22,32,33]. For instance, Budge et al. [21] used five calibration scenarios (four sets of CC were determined by captive feeding experiments, and the fifth calibration scenario was modeled without CC) to verify the application of QFASA in fish. They found a tendency to overestimate the dietary components used as prey for CC calculations, which indicates that applying QFASA without CCs to quantitatively estimate the dietary components

and the trophic relationships in aquatic food web might be sufficiently robust. Magnone et al. [14] generated a food web model in the aquatic environment without CCs, found the results were similar to the previous knowledge; Happel et al. [22] illustrated that the prey used to make CCs can affect QFASA outputs and that perhaps ignoring CCs was a better method. Moreover, in a latest research, Bromaghin et al. [33] proposed a model that can simultaneously estimate CC and dietary composition on the basis of FA signature samples of predator and its preys alone. In summary, the application of CC depends on the species and system being studied, and the diet estimation of fish without CC is sufficiently robust, and for other taxa which are difficult to conduct captive feeding trials, it is recommended to refer to the model proposed by Bromaghin et al. [33].


**Table 1.** Captive feeding trials to improve the precision of quantitative fatty acid signature analysis (QFASA).

#### 3.1.4. The Selection of FA Subsets

The second issue, as stated above, is related to the FA subset used in assessing diets. More than 70 FAs have been identified through the gas chromatography and a polar capillary column in marine lipids studies [12]. Nevertheless, due to the influences of the predator's metabolisms, only some FAs can provide useful information for dietary estimation [38]. For example, short-chain or medium-chain FAs (i.e., <14 carbons) could only be produced by biosynthesis. In contrast, n-6 or n-3 polyunsaturated FAs usually only come from the diet. Due to these differences in origins, some studies have explored the utility of FA subsets rather than the full FAs.

Traditionally, the dietary or extended dietary subsets of FA investigated by Iverson et al. [5] is routinely used [10,11,13,16,30]. The former FAs are strictly from the diet, while the latter included several FAs that can be biosynthesized by the predator or ingested from diet [5]. In general, the extended dietary subset performed better in the modeling exercises, likely due to the increased information provided [8]. Besides, additional criteria were also applied by some investigators [18,19,21]. For example, Magnone et al. [15] selected seven FA sets to estimate diet; Budge et al. [21] used four FA sets; Goetsch et al. [19] adopted a new FA subset, which only included FAs derive from the diet (*n* = 46) and removed three FAs (16:4n-3, 18:1n-11, 22:1n-11) that were used in other studies. In conclusion, the optimum subsets of FA may depend on the type of predators (e.g., mammal, bird, fish) and the ecosystem studied (e.g., estuarine ecosystems, fresh water, and soil ecosystem). It is a significant area for further research.

There are generally three approaches to deal with the FA subset selected [13]: the "traditional" approach is to rescale the sum of the proportions of the FA components to 1 [5,13,19]; the "un-scaled" approach maintains the original ratio of FAs in the diet unchanged; and the "augmented" method is to add an additional ratio while maintaining the original proportions of FAs in the subset to make the sum equal to 1.

Bromaghin et al. [13] used computer simulations to compare the effects of three scaling methods based upon *KL* [5], Aitchison [39] and Chi square [39] distance measures on the variance and bias of dietary estimation. They recommended that researchers should check the sum of the original proportions and analyze the degree of difference among preys, and then determine whether the scaling will cause deviation.

#### 3.1.5. The Optimization of the Statistical Model

QFASA has been extensively used since it was proposed and some distance measurement methods have been applied in previous studies (e.g., [5,19,39,40]), such as the *KL* [5], Aitchison [19,39], Chi square [39] statistical distance, although the *KL* measure originally recommended by Iverson et al. [5] has been frequently used [39].

Dietary estimation was usually carried out in two spaces. One was the prey space, which was transformed by dividing predator FA signature by CC based on Formula (1) [5]; and the other was the predator space, which was obtained by multiplying prey FA signature by CC based on Formula (1) [18,19]. Bromaghin et al. [41] investigated the influence of the choice of distance measures (such as *KL* or Aitchison) and optimization spaces (i.e., prey or predator space mentioned above) on dietary estimation in the application of the QFASA model, and revealed that the choice of evaluation method may have a significant impact on dietary estimation. Additionally, Bromaghin et al. [42] developed a new algorithm that can objectively determine the bootstrap sample size to generate a pseudo-predator signature with actual attributes, thus improving the efficacy of computer modeling to evaluate the performance of the QFASA estimator. Besides, Bromaghin et al. [40] used computer modeling to study the robustness of the Aitchison and *KL* distance measure, kept a record of the deviations in dietary estimation, and concluded that the former was more robust in terms of CC errors while the latter was more robust to the ingestion of preys unrepresentative in the prey profiles database. Furthermore, Litmanen et al. [43] estimated the performance of several algorithms and found that some Bayesian algorithms take longer time to calculate than QFASA, recommending the use of Chi square or *KL* statistical distance. In conclusion, the distance measure that performs the best may depend on the ecosystem being studied, and it is recommended to test multiple measures to assess differences, or evaluate the ecosystem under controlled conditions (e.g., feeding experiments).

In the application of the QFASA model, the statistical characteristics of predator dietary estimation is usually estimated through computer simulations [42], e.g., commercial software, Fortran programs [18], R package [13,19,30,32,44,45], a combination program of R and Fortran [6,41], and Matlab with its optimization toolbox [14,15,17,36]. In recent years, Bromaghin summarized a new R package named QFASAR, calculating the goodness-of-fit diagnosis, which may enhance the performance of the prey signature database [19,44].

#### *3.2. Research Status of QFASA*

QFASA has become a promising approach in dietary estimation since it was proposed and its applications are shown in Table 2. In recent years, Wang et al. [11] estimated the diets of threatened *S. fischeri*, and concluded that infertile eggs yolks can be used to evaluate the diet of breeding female *S. fischeri* for better understanding the source and timing of nutrients during reproduction; Magnone et al. [15] estimated the diet of *P. orbignyanus* by QFASA, determined CC and validated the model with a controlled experiment, which quantified the diet of lower vertebrates (*P. orbignyanus*) for the first time; Magnone et al. [14] used QFASA to generate a food web model in the aquatic environment to find out the trophic relationships among the species; Conners et al. [20] estimated two albatross species by QFASA, and adapted the QFASA model, which introduced a combination of FAs and fatty alcohols; Goetsch et al. [19] estimated the diet composition of *M. angustirostris*, and suggested that the seals mainly fed on mesopelagic fishes; Knox et al. [24] estimated the diet of *A. pusillus doriferus*, finding that elasmobranchs accounted for more in the diet of males than previous reports, and showed that prey composition varied among males; Happel et al. [25] used a controlled trial to test the influence of intraspecific difference in FA signatures of prey on the QFASA model, indicating that QFASA used for steelhead trout may not only be used for a specific lake, but also for other freshwater systems with alewife and round goby as the main food; Tao [26] estimated the diet composition of *Channa. argus* using the QFASA model, and then calculated the biomagnifcation factors (BMFs) of the alternative halogenated flame retardants (AHFRs) based upon the results.




**Table 2.** *Cont.*

#### **4. Conclusions and Suggestions**

FAs have been widely applied in qualitative research for evaluating trophic predator–prey relationships. However, QFASA is the first attempt to use FAs to quantitatively estimate the diet of predators, which provides an insight into the function and structure of dynamic ecosystems. The core of QFASA is to find the FA signatures combination of prey that is most similar to the FA signatures of predator to infer the diet of predator. One of the most important assumptions of the QFASA model is that the metabolism of the predator is already known. Currently, the effect of predators on FA metabolism and deposition is explained by CCs, which were determined previously from captive feeding trials and could be found by calculating FA levels found in the predator over FA levels in the food. Additionally, CCs have already been evaluated in several species of marine mammals, seabirds and fish (Table 1). The refinement of which subsets of FAs are more effective for QFASA accuracy remains a topic of study. QFASA has been used in a variety of taxa (e.g., marine mammals, seabirds and fish) rather than being limited to high-trophic vertebrates, and presumably it will be tested in new taxa. If one predator in a given ecosystem meets all of the above requirements for using QFASA, it can be modeled.

Nevertheless, QFASA has several limitations. Firstly, since it is a non-probabilistic model, it is hard to evaluate the effects of the different uncertainties related with dietary proportion estimation. The handling of confidence interval in QFASA can be found in Stewart et al. [51,52]. Ecological mechanisms cannot be built directly into the model owing to the absence of an explicit model [53]. Secondly, some species are not suitable for captive experiments; therefore, accurate CCs cannot be obtained. Thirdly, the model is not applicable to all species, because some taxa lack specific fat storage organs. In addition, some factors may cause bias when QFASA is applied to fish and invertebrates, since FAs may be modified during metabolism and transportation. Some invertebrates tend to have greater ability to biosynthesize and modify the FA than higher trophic organisms, and fish have a stronger ability to modify some exogenous FAs than mammals and birds [7]; many studies have indicated that FAs with a high proportion in the diets of fish are easily catabolized [21]; the growth and reproductive stage (e.g., the proportions of 20:5n-3 and 22:6n-3 in fish eggs seems to be fixed) may affect FA signature of tissues [21]; some freshwater taxa (e.g., Daphnia) can elongate and desaturate 18:3n-3 to 20:5n-3 [7].

Several techniques have been used to study predator–prey relationships, such as stomach content analysis, IS analysis and FA analysis [54]. Traditional stomach content analysis can directly reflect the diet composition of predator, but it usually only represents a snapshot of diet and tends to underestimate the proportion of soft-bodied prey [55]. Chemical markers such as SI and FAs can reflect the long-term diet of predator. SI analysis plays an important role in evaluating the trophic level of predators, but the resolution is limited by the number of SI to be measured (e.g., typically only 2–3 SI are measured) [56]. FAs seem to provide more information than SI, because many potential FAs can be measured. However, the proportion of the measured FAs is always limited to a sum of one. Therefore, some studies have developed mixed models that combine SI and FAs, such as the Bayesian fatty acid-based mixing model (Fatty Acid Source Tracking Algorithm in R, FASTAR) [57,58] intended for zooplankton and benthic macro-invertebrates, a mixed model combining FAs with SI (FastinR), which could improve the dietary estimation with the available fat content and conversion coefficients [53], and the Bayesian mixing model framework (mixSIAR) [59], which evaluates the relative contribution of food sources in the diet of predator. As a promising method in quantitative ecosystem, QFASA relies on the distance measurement rather than model-based formulation to evaluate the most likely diet proportions, which is contrary to (Bayesian) SI and FA mixing models.

As for the future studies of QFASA, firstly, it is suggested to use QFASA in combination with other approaches (such as the SI analysis [8] or stomach content analysis [16]) to obtain more comprehensive information. Secondly, some simplification procedures of FAs extraction are expected to widen the use of the QFASA model in marine ecology and biology, such as the simplified method to extract polyunsaturated FAs [60]. Thirdly, the application of QFASA in invertebrates and low-trophic vertebrates needs further study. More studies are expected to use the QFASA model to build food

webs [14,61]. For instance, Magnone et al. [14] generated a food web model in the aquatic environment to find out the trophic relationships among the species by QFASA. Further, captive feeding trials may provide a good way of testing if QFASA will work and its intricacies, but the predator signature produced through trials may not relate to wild species of different age, size, or lipid status even if they are fed the same prey in the wild. There is still much remaining to be known and tested with how FAs can be used in quantitative models. For instance, Goetz et al. [62] indicated that within the lake trout species different ectomorphs accumulate lipids in different ways which may mean finding CCs for plastic species may be difficult and suggests that over time CCs would need to be updated as genetics drift.

In conclusion, it has been clearly shown that FA has become an important trophic tracer in the research of carbon transfer, predator–prey relationships, food webs, along with the function and structure of dynamic ecosystem. Although QFASA has some limitations, it is still the promising approach because of its several potential advantages [12]. QFASA reflects the long-term (a period of weeks to months) diet compared with the stomach content analysis, allows more than three prey types to be estimated, and can avoid the problem of underdetermined systems common in analysis with SI markers. Results of QFASA may be useful in studies into the effects of maternal diets on offspring or the assessment of nutrient deficiencies in marine organisms [63,64].

**Author Contributions:** Conceptualization, J.Z.; methodology, J.Z. and C.R.; data analysis, C.R., H.Z., F.Y. and S.Z.; resources, H.Z. and S.Z.; writing—original draft preparation, J.Z. and C.R.; writing—review and editing, J.Z., C.R. and D.K.; supervision, R.W. and D.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is partly supported by the National Key Research and Development Program of China (No. 2019YFC0312104), the Young Orient Scholars Program of Shanghai (No. QD2017038), and the National Natural Science Foundation of China (No. 41807341).

**Acknowledgments:** The authors would like to thank Yang Chenxing for discussing the English grammar of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Aspects of Reproductive Biology of the European Hake (***Merluccius merluccius***) in the Northern and Central Adriatic Sea (GSA 17-Central Mediterranean Sea)**

**Michela Candelma 1, Luca Marisaldi 1, Daniela Bertotto <sup>2</sup> , Giuseppe Radaelli 2, Giorgia Gioacchini 1, Alberto Santojanni 3, Sabrina Colella <sup>3</sup> and Oliana Carnevali 1,\***


**Citation:** Candelma, M.; Marisaldi, L.; Bertotto, D.; Radaelli, G.; Gioacchini, G.; Santojanni, A.; Colella, S.; Carnevali, O. Aspects of Reproductive Biology of the European Hake (*Merluccius merluccius*) in the Northern and Central Adriatic Sea (GSA 17-Central Mediterranean Sea). *J. Mar. Sci. Eng.* **2021**, *9*, 389. https://doi.org/ 10.3390/jmse9040389

Academic Editor: Francesco Tiralongo

Received: 26 February 2021 Accepted: 31 March 2021 Published: 7 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** The study focused on the macroscopic, histological, and biometric analysis of European hake females in GSA 17 (Central-North Adriatic Sea). From 2013 to 2015, 976 females were collected and analyzed monthly. Though females in spawning conditions were found during the whole year, the trend of GSI showed a peak of the reproductive season from April to July in 2014 and 2015. HSI and *Kn* reached the highest values in September, after the spawning peaks. In 2013, the trend of these indices did not highlight an evident peak, probably due to an adverse event that occurred in the previous winter in the Adriatic shelf. The length at first maturity (L50) was estimated by macroscopic and histological approaches, resulting in 30.81 cm for the macroscopical length and 33.73 cm for the histological length; both values are higher than the current catching legal size. For the first time in this area, batch and relative fecundity were estimated. Relative fecundity was similar to the Mediterranean and the Atlantic stocks, whereas batch fecundity values were lower compared to other fishing grounds. Overall, the analysis of reproductive parameters plays a fundamental role in the sustainable management of this resource in an area as overfished as the Central-North Adriatic Sea.

**Keywords:** European hake; *Merluccius merluccius*; fecundity; somatic indices; Adriatic Sea; L50

#### **1. Introduction**

In the last several decades, the excessive fishing effort, together with the increase of pollution, poor fishing management, and impairment of marine ecosystems, caused the depletion of fish stocks worldwide. The Food and Agriculture Organization reported that 57.4% of fish stocks are fully exploited, 29.9% are overexploited, and only 12.7% are not fully exploited [1]. In the Mediterranean Sea, the long-lasting, intense fishing pressure applied on fish and invertebrate stocks has led to declining population biomasses [2], which has also been reflected in a reduction of catches for the majority of stocks [3]. Mullon et al. [4], by analyzing the FAO dataset of world fisheries catches for the period 1950–2000, detected that the major collapse occurred for demersal species. Colloca and coworkers [5] assessed the impacts of the fishing pressure in the period 2002–2014 in the Mediterranean Sea and determined that the Central-North Adriatic represents the highest catching area of demersal species in the Mediterranean area.

Among demersal species, the European hake (*Merluccius merluccius*, L. 1758) received attention because it represents one of the principal fishing targets in the Northeast Atlantic Ocean and the Mediterranean Sea. In the Mediterranean Sea, a total landing of 22,547 tons was recorded in 2011 [6,7]; in the Central-North Adriatic Sea, the European hake is one of the leading commercial species [8] and represents 77% of landings from Croatia [9]. To avoid the collapse of hake stock, fishing should be conducted more sustainably by maintaining a spawning stock biomass level that is suitable to guarantee the renewal of the species. In this regard, egg-producing females contribute more to embryonic development success. The acquisition of information on reproductive cycle of the female Adriatic stock, including length at first maturity, fecundity, and spawning cycle, is essential to quantify the reproductive potential of this population [10]. Because of the great commercial importance of the European hake and the critical status of its stocks, several studies were performed in different areas of the Mediterranean Sea related to stock assessment, spawning cycle and fecundity estimation [11–16], feeding habits [17,18], analysis of lipid reserves [19], nursery area [6,20], juvenile recruitment [21], the selectivity of fishing gear [22], fisheries management [2,23,24], and reproductive physiology [25,26]. However, none investigated the fecundity, the spawning cycle, and the proportion of mature individuals at length by a macroscopic and histological approach in the Central-North Adriatic Sea (FAO Geographical Sub-Area 17, according to GFCM division), where this species represents a valuable economic resource.

In this regard, a multiannual survey (2013–2015) was conducted to improve understanding of the reproductive dynamics of European hake in this area. In particular, the work focused on females since the female conditions limit the egg production and the progeny production more than males, and it evaluated the size at first maturity in females using the histological and macroscopic analysis of the gonads. The condition indices such as Le Cren's condition factor (*Kn*), the hepatosomatic index (HSI), and gonadosomatic index (GSI) were calculated and compared across different ovarian stages, and their changes analyzed in the three years. Furthermore, batch and relative fecundities were estimated together with the analysis of ovarian stages.

This study has never been carried out in this area (GSA 17) and completes the scenario of this important resource of the Italian seas, in order to provide knowledge for a sustainable management.

#### **2. Materials and Methods**

#### *2.1. Sampling*

A total of 976 females ranging from 13 to 64 cm in total length (TL) were collected aboard commercial bottom trawler fishing vessels in the Northern and Central Adriatic Sea (FAO Geographical Sub-Area 17, according to GFCM division) monthly from January 2013 to December 2015. Guidelines of the Data Collection Framework Regulation (EU Reg.199/2008) followed those established by the Community system for the conservation and sustainable exploitation of fisheries resources under the Common Fisheries Policy (CFP). The sampling procedures did not include any animal experimentation and animal ethics approval was therefore not necessary under the Italian legislation (D.L. 4 March 2014, n. 26, art. 2). For each specimen, the following parameters were recorded: total weight (TW), gutted weight (GW), total length (TL), sex, macroscopic maturity stage of gonads, and gonad and liver weights.

#### *2.2. Reproductive Seasonality and Fish Condition Indices*

To evaluate the temporal variation of maturity and condition stage of females, the gonadosomatic index (GSI), hepatosomatic index (HSI), and Le Cren's relative condition factor (*Kn*) were calculated. These indices were defined by the following equations:

GSI (%) = gonad weight/gutted weight body\*100

HSI (%) = liver weight/gutted weight body\*10

*Kn* = W/*a* TL*<sup>b</sup>*

where *a* and *b* are the regression parameters of the length–weight relationship, W is gutted weight, and TL is the total length. GSI was evaluated only for the mature females; 313 mature animals were found from developing to regenerating phase. From January to March of 2013, livers were not available for sampling problems. The spawning season was also investigated by analyzing the monthly frequency of ovarian maturity stages. Number of samples of mature females used for GSI, HSI and Kn calculation for each month in three years are reported in Table S1.

#### *2.3. Histological Analysis*

173 female specimens at different maturity stages were randomly chosen for histological analysis while choosing a group that was representative for all sizes to confirm the macroscopic classification of ovarian development. Ovarian samples were fixed at 4 ◦C overnight in formalin at 10% in phosphate-buffered saline (PBS, 0.1 M, pH 7.4), after which they were washed in PBS and stored in 70% ethanol at 4 ◦C until use. The samples were then dehydrated through a series of graded ethanol, cleared in xylene, and embedded in paraffin wax. Five μm thick consecutive sections were cut using a microtome RM2125 RTS (Leica Biosystems, Wetzlar, Germany) and stained with Meyers's hematoxylin and eosin. Images were acquired using a Zeiss Axio Imager M2 microscope (Zeiss, Oberkochen, Germany). Assessment of the reproductive status of samples was performed according to the method of Brown-Peterson and coworkers [27] but adapted for European hake following the method of Candelma and coworkers [25]. Five ovarian stages were defined: (i) immature; (ii) developing; (iii) spawning capable and actively spawning; (iv) regressing (postspawning); (v) regenerating (spent) (Table 1).


**Table 1.** Criteria used to determine the maturational status of European hake females.

The measurement of oocytes was only from oocytes that had been sectioned through the nucleus. The oocyte classification followed Candelma et al. [26] and Murua et al. [28]: oogonia (O); primary oocyte (PO); cortical alveoli (CA); lipid stage (LS); early vitellogenesis (vtg 1); middle vitellogenesis (vtg 2); late vitellogenesis (vtg 3); migrating nucleus (Mn); hydration (H); postovulatory follicles (POFs). Cohen's *k* coefficient [29] was applied to assess the agreement between the histological and macroscopic classification.

#### *2.4. Size at First Maturity (L50) and Fecundity*

The proportion of maturity at length (PL) was estimated using both histological and macroscopic data with the following logistic function:

$$\text{PL} = \frac{1}{1 + \exp(\alpha + \beta \* TL)}$$

where *α* (intercept) and *β* (slope) represent the estimated parameters and *TL* is the total length. All specimens were used for macroscopic investigation; 173 individuals were examined using histological analysis. Females from the developing stage onwards were considered mature. The length at which 50% of the females are mature was computed as L50 <sup>=</sup> −*α*/*β*. The two regression logistic curves of histological and macroscopic data were compared with the likelihood ratio test.

With regards to batch and relative fecundities, 28 females at the actively spawning stage representing the total number of females available were used and the gravimetric method, based on the relationship between ovary weight and oocyte density as described by Murua et al. [30], was applied. Because the distribution of hydrated oocytes does not statistically vary in the ovaries [31], only a lobe of the hydrated ovary was used for analysis. The hydrated oocytes from fresh subsamples were manually counted under a stereomicroscope (Optika, Italy). Batch fecundity was estimated as the average of the hydrated oocyte number integrated over three subsamples multiplied by ovary mass for each specimen with the hydrated oocyte [11]. The relative fecundity was calculated as a ratio of batch fecundity and gutted weight for each fish.

#### *2.5. Statistical Analysis*

Statistical differences of the monthly variation of indices for each year and of indices per ovarian stages in each year were determined using a Tukey's multiple test comparison. Both analyses were performed using Prism 6 (GraphPad Software, San Diego, CA, USA). For the determination of size at first maturity, the R statistical environment [32] and the packages "Fisheries Stock Assessment" (FSA) [33], ggplot2 [34], and rel [35] were used. The *p*-values < 0.05 were considered as significant. Results were expressed as the mean ± SEM.

#### **3. Results**

#### *3.1. Reproductive Seasonality and Fish Condition*

The percentages of European hake females in different maturity stages observed monthly in the three years indicated that specimens from developing to regressing phase were present throughout the year (Figure 1).

In Figure 2, the GSI trend of mature females during three years of sampling is represented. The GSI values were lowest in 2013 compared to the other two years. In 2013, the highest value was recorded in April with a significant difference only compared to October (Tukey's multiple comparison test; *p* < 0.05) (Figure 2a). Differently from 2013, in 2014 and 2015, the highest peak was observed in June and the lowest in September (Figure 2b,c). In particular, in 2014, the GSI in June was significantly higher than February, September, October, November, and December (Tukey's multiple comparison test; *p* < 0.05) (Figure 2b), and the peak in July was significantly higher than September, December, and February (Tukey's multiple comparison test; *p* < 0.05). In 2015, June was significantly different only compared to September (Tukey's multiple comparison test; *p* < 0.05) (Figure 2c).

**Figure 1.** Monthly percentages of European hake female maturity stages. Imm/reg, includes individuals in immature and regenerating stages; dev, developing; spaw, spawning capable; ac spaw, actively spawning, regr, regressing. For their histological characteristic, the immature and regenerating stages were grouped (See Section 3.2).

**Figure 2.** Monthly variation of European hake female gonadosomatic index (GSI) for three years. Statistical differences were described in the result section. In January and March of 2013, no mature females were found. In August of the three years, no samples (n.s.) were sampled. The letters indicate a statistically significant difference among groups (*p* < 0.05) determined by Tukey's multiple comparison test.

The HSI analysis distinguished immature females from mature females (considering from developing to regenerating phase) (Figure 3). The HSI trend varied over the three years. In immature females from 2013, the highest values were in summer months, but significantly different only with September (Tukey's multiple comparison test; *p* < 0.05) (Figure 3a). In 2014, the highest peak was reached in December, significantly different compared to March (Tukey's multiple comparison test; *p* < 0.05) (Figure 3b), whereas in 2015, the highest value was in September and significantly different only compared to March (Tukey's multiple comparison test; *p* < 0.05) (Figure 3c). The lowest value was in March in both 2014 and 2015. In 2013, as in immature females, the mature females were also characterized by an increase of HSI values from May to July (Figure 3d), whereas the lowest measurement was in April significantly diverse compared to July, September, October, November, and December (Tukey's multiple comparison test; *p* < 0.05). Different peaks were recorded during 2014, the highest in December, significantly higher than March– July (Tukey's multiple comparison test; *p* < 0.05) (Figure 3e). In 2015, the HSI surged in September, significantly different compared to all months (Tukey's multiple comparison test; *p* < 0.05), except to January and December (Figure 3f), and the lowest value was in May.

**Figure 3.** Monthly variation of European hake female hepatosomatic index (HSI) for three years in different ovarian stages. Graphics (**a**–**c**) indicate the females in the immature phase, while graphics (**d**–**f**) show females from developing to regenerating phase. From January to March of 2013, no sampling of the liver was conducted. In August of the three years no samples (n.s.) were sampled. The letters indicate a statistically significant difference among groups (*p* < 0.05) determined by Tukey's multiple comparison test.

The *Kn* analysis distinguished immature females from mature females (considering from developing to regenerating phase) (Figure 4). In immature individuals collected in 2013, the highest peak was in October and significantly different only compared to April (Tukey's multiple comparison test; *p* < 0.05) (Figure 4a), while in mature animals no significant variation was observed (Figure 4d). The trend of *Kn* varied mainly in 2014 and 2015, with the highest peak recorded in September for both immature and mature females (Figure 4b–f). In immature individuals captured in 2014, the lowest value was in March and significantly different compared to September and November (Tukey's multiple comparison test; *p* < 0.05) (Figure 4b), whereas mature females showed the lowest value in June with significant differences respect to September and December (Tukey's multiple comparison test; *p* < 0.05) (Figure 4e). In 2015, the *Kn* of immature females was significantly higher in September compared to the rest of the year except for May, June, July, and October (Tukey's multiple comparison test; *p* < 0.05) (Figure 4c). In mature females captured in 2015, the highest peak of September was significantly different compared to January, April, and May (Tukey's multiple comparison test; *p* < 0.05) (Figure 4f), whereas the value was lowest in December.

**Figure 4.** Monthly variation of European hake female Le Cren's condition factor (*Kn*) for three years in different ovarian stages. Graphics (**a**–**c**) indicate the females in the immature phase; graphics (**d**–**f**) show females from developing to regenerating phase. In January and March of 2013, mature individuals were not found. In August of the three years, no samples (n.s.) were sampled. The letters indicate a statistically significant difference among groups (*p* < 0.05) determined by Tukey's multiple comparison test.

The analysis of the condition indices per ovarian stage per three years was macroscopically performed (Figure 5a–c). For their histological characteristic, the immature and regenerating stages were grouped. Regarding the GSI, the significantly highest peak was evidenced in the spawning phase for all three years (Tukey's multiple comparison test; *p* < 0.05) (Figure 5a), the trend increased from imm/reg to the spawning stage for all three years but was significant only in 2013 (Tukey's multiple comparison test; *p* < 0.05). For HSI, the tendency of mean values was lowest in the imm/reg stage for all three years (Tukey's multiple comparison test; *p* < 0.05) (Figure 5b). The other mean values were not significantly different among them within the same year, except for 2014, in which the spawning phase showed a significant decrease (Tukey's multiple comparison test; *p* < 0.05) compared to developing and regressing stages. The variation of *Kn* was less visible compared to the other two indices in the period analyzed (Figure 5c). The only significant difference (Tukey's multiple comparison test; *p* < 0.05) was evidenced in the spawning capable stage compared to the imm/reg and developing stages in 2014 and only compared to imm/reg phase in 2015.

#### *3.2. Histological Analysis: Ovarian Classification and Patterns*

The female reproductive cycle of the European hake consists of five ovarian reproductive stages according to the histological classification (Table 1). Ovaries from immature individuals exhibited compact ovigerous lamellae with groups of oogonia and previtellogenic oocytes (diameter < 250 μm) and an ovarian wall usually thinner than that of the regenerating stage (Figure 6a). Entry into the developing phase is characterized by the appearance of early and the middle vitellogenic oocytes with a diameter ranged 250–550 μm (Figure 6b). The oocytes take up the yolk proteins. The yolk granules multiply and increase in size, forming a densely packed zone in the inner part of the cytoplasm. The oil droplets proliferate, becoming lipid globules and the zona radiata increases in width. The ovaries from spawning specimens were characterized by the presence of fully grown vitellogenic oocytes (vtg3) and/or migrating nucleus (Mn) stage as well as hydrated oocytes (H) (diameter 550–1150 μm) (Figure 6c–e). Nuclear migration and/or hydration distinguish the actively spawning subphase (Figure 6d). The presence of postovulatory follicles in the same ovary with vtg3, OM, and H was observed (Figure 6e). Extensive atresia and a reduced number of vitellogenic oocytes were considered markers of the regressing stage (Figure 6f). The females in the regenerating phase have spawned at least one time in life. Their ovaries are loose and have thickened wall and CA and LS oocytes are present among previtellogenic oocytes (diameter < 250 μm) (Figure 6g).

**Figure 5.** Values of the GSI (**a**), HSI (**b**) and *Kn* (**c**) for females of European hake per ovarian stages for three years. Imm/reg, includes individuals in immature and regenerating stages; dev, developing; spaw, spawning capable; regr, regressing. The letters indicate a statistically significant difference among groups (*p* < 0.05) determined by Tukey's multiple comparison test.

**Figure 6.** Tissue sections of European hake ovaries at different developmental stages. (**a**) Immature; (**b**) developing; (**c**) spawning capable; (**d**,**e**) actively spawning; (**f**) regressing; (**g**) regenerating. For a,b,g scale bar = 100 μm; for c,d,e,f scale bar = 200 μm. The abbreviations indicate: OW, ovarian wall; O, oogonia; PO, primary oocytes; CA, cortical alveoli; LS, lipid stage; vtg1, vitellogenesis 1; vtg2, vitellogenesis 2; vtg3, vitellogenesis 3; Mn, migrating nucleus; H, hydrated oocytes; POF (postovulatory follicle); Aα, alpha atretic oocyte; Aβ, beta atretic oocyte.

The histological investigation showed a 60.6% similarity with the macroscopic maturity stage classification of the females. Cohen's *k* was 0.45 (95% confidence interval: 0.35–0.55), which corresponds to a "Moderate" level of agreement [36]. Taking into consideration the macroscopic classification, the percentage of agreement between the histological and macroscopic analysis was different among phases. For the imm/reg phase was 91.8%, for the developing stage was 45.8%, for the spawning phase was 85.3%, and for regressing was 2.86%.

#### *3.3. Size at First Maturity (L50) and Fecundity*

L50 was estimated macroscopically and histologically to be 30.81 cm and 33.73 cm (Figure 7), respectively.

**Figure 7.** Estimated L50 and observed the proportion of mature individuals by 1 cm length classes. (**a**) Macroscopic value; (**b**) histological value. The shaded area indicates a 95% confidence interval.

The estimated parameters of the logistic regression were statistically significant (*p* < 0.05) and summarized in Table 2.


**Table 2.** Summary of the L50 estimates based on histological and macroscopic data.

The two curves of L50 were not significantly different (*p* > 0.05). The size at which all the females reproduced at least once was 39 cm. The shortest length at which the female specimens entered ovarian development (vitellogenesis) was 26 cm.

Batch fecundity estimation ranged between 7771 (TL 33 cm) and 137,256 (TL 32 cm) hydrated eggs. The shortest total length of a specimen with hydrated eggs in the ovary was 29 cm (batch fecundity of 45,004 eggs) and the longest was 42 cm (batch fecundity 70,519 eggs). The mean value of a relative fecundity was 205 eggs/g of body weight, the lowest value of relative fecundity was 34 eggs/g of body weight found in a female with TL of 33 cm, and the highest value of relative fecundity was 573 eggs/g of body weight derived from a female with TL of 32 cm.

#### **4. Discussion**

One of the targets of fisheries management is to conserve enough reproductive potential in a population to allow for sustainable exploitation [37]. To achieve this, the reproductive parameters, such as length at first maturity, fecundity, and spawning cycle determination, are usually applied as health indicators of the stock. The present study is the first attempt to highlight an exhaustive knowledge of the reproductive biology of *Merluccius merluccius* in the Northern and Central Adriatic Sea (GSA 17) to complete the scenario of knowledge in the Italian seas.

The easiest, cheapest, and most direct approach to investigate some aspects of reproductive biology is a macroscopic evaluation of the maturity stages of the gonads, but it is not always an accurate method [16,38,39]. This study considered the somatic indices to evaluate the reproductive state of the European hake. We focused on different indices, namely the gonadosomatic index (GSI), the hepatosomatic index (HSI), and Le Cren's condition factor (*Kn*) along three consecutive years. The GSI trend calculated on mature specimens sampled from the GSA17 indicated that the hake reproductive season reached a peak from March to July, even though a presence of spawning females was observed throughout the year, confirming that a protracted spawning period is typical for this species [11,40,41]. Unlike 2014 and 2015, in 2013, a lack of the GSI peak during the reproductive season was observed. The possible reason for this could be related to an exceptional event of dense water that occurred in the winter of 2012 in the Adriatic shelf [42]. Such a phenomenon was characterized by high-velocity currents that caused a cascading event leading to the transport of more than 50% of water volumes from the Northern Basin to the Southern area of the Adriatic Sea [42–44]. Directly or indirectly, unfavorable conditions for the successfulness of the hake reproduction could have emerged, as previously suggested by Recasens et al. [11], in which dense water formation caused intense cascading events shifting the hake reproduction in the Catalan Sea by several months. Thus, the oceanographic events that took place in the Adriatic Sea in 2012 could have negatively impacted the stock of the European hake in 2013 in the area. Regarding the GSI value in different ovarian stages, the

highest mean value was reached at the spawning stage, confirming its correlation with the sexual maturity of gonads.

In the European hake, lipids are stored mainly in the liver, confirming the critical role of this organ for energy storage in this species [19]. The analysis of HSI represents a good proxy of the energy reserve available in fish. In the Central-North Adriatic Sea, HSI showed different trends within the three years studied in both immature and mature females. In 2014 and 2015 for all ovarian phases, the highest recorded values of HSI after summer are probably due to the increase in the food availability for this demersal predator during this period. The energy stored during richer months foodwise is not followed by reproductive activity, allowing the European hake to survive during the winter period, as previously reported in Galician shelf stock [45]. In 2013, the HSI reached the greatest value of the year in the summer months, probably because the absence of the reproductive activity allowed an incessant accumulation of lipids. Based on ovarian developmental stages, the lowest level of HSI was obtained in the imm/reg stage, demonstrating that juvenile fish use lipid energy for growth and do not accumulate them in the liver [46].

Together with HSI, Le Cren's condition factor indicates the energy storage in fishery species. It is given by weight and length data and assumes that heavier fish of a given length are in a better condition [47]. In 2014 and 2015, the *Kn* showed the highest values in September both for immature and mature individuals, just like HSI. In 2013, the autumn peak shifted one month only for immature females, whereas no variation was shown in mature animals, so we speculated that like for HSI, the *Kn* also indicates an incessant accumulation of energy in fish during such a period. In the spawning capable phase, the low levels of *Kn* occurred in 2014 and 2015; only HSI in 2014 could confirm that European hake uses part of the energy accumulated in previous months for reproduction if it is necessary.

The assessment of maturity stages of the ovaries is a crucial step for the estimation of the L50, and errors might lead to subsequent biased estimation of the spawning stock biomass [38,39,48]. In this contest, an accurate assessment of sexual maturity is an essential component of effective fisheries management. For European hake, some sexual maturity classifications exist, but often in these studies, a different terminology was used for the oocyte stages or ovarian phases [12,16,28,49]. Furthermore, some authors proposed only an oocyte classification, whereas other authors suggested a macroscopic and histological scale, but without a proper calibration between visual (macroscopic) and histological (microscopic) staging. We used the oocyte stage classification proposed by Candelma et al. [26] and Murua et al. [28] and the standardized terminology for ovarian phases described by Brown-Peterson et al. [27] to obtain a clear and accurate scale on ovarian and oocyte stages validated by the similarity between macroscopic and histological approach. The "Moderate" level of agreement in the assignment of reproductive stages between histological and macroscopic investigation revealed that our macroscopic classification for European hake does not suffer from relatively high error rates. An exception was the regressing phase, whose percentage of agreement was weak. Overall, our scale appeared an accurate tool for the sexual maturity determination by the macroscopic method, notwithstanding some stages require (i.e., the regressing phase and the immature stages) histological analysis to avoid misclassification. Moreover, the histology is the most suitable method for the studies on reproductive biology and it is strictly recommended a regular calibration (over the years) between the two approaches. Histological analysis confirmed that *M. merluccius* is a batch spawner with an asynchronous type of ovarian development.

The L50 is a fundamental parameter of fish population dynamics studies and plays a crucial role in species management purposes. In the present study, the length at first maturity was estimated as 30.81 cm by macroscopic evaluation of gonads maturity and 33.73 cm by histological analysis. The L50, determined in this work, is comparable to the value of 30.5 cm obtained in a previous study of Zupanovic and Jardas [50] in the same area, while Alegria Hernandez and Jukic [51] estimated an L50 equals to 31.3 cm in total length, but both only by macroscopic inspection. In the Southern Adriatic, the size of the first maturity calculated using the macroscopic method was 31.95 cm [16]. In comparison with previous studies on length at first maturity in the Mediterranean Sea, the values estimated were similar to the measure reported in the present study. On the Egyptian coast, Al-Absawy [14] found an L50 of 32.5 cm, while in the Algerian coast [13], it was 30.6 cm. The first maturity size is generally affected by variations of several environmental factors that characterize the different areas, such as the abundance and distribution of local populations, competition for space, and availability of nourishment. Furthermore, the overexploitation of this resource could explain the differences in this parameter that reduces the spawning biomass in different areas [12,16]. The L50 of 33.73 cm resulting from histological analysis evidenced that the macroscopic approach underestimates this value, but the two curves were not significantly different, indicating that the histological examination is not strictly necessary.

The L50 found in the present study contrasts the decision regulated by Annex III of the Council Regulation (EC) N◦ 1967/2006, in which the minimum legally allowed fishing size is 20 cm, implying the catch of sexually immature individuals. This excessive fishing effort could lead to a likely decrease in the stock in the area in upcoming years. Accordingly, as reported by Anderson et al. [52], the reduction of the age of a stock and average body size can reduce the ability of exploited species to survive to annual environmental variation.

In the present study, for the first time, batch and relative fecundities in the Northern-Central Adriatic Sea were analyzed. The mean value of a relative fecundity reported in this study (205 eggs/g of body weight) seemed comparable with the previous estimates reported in Catalan waters (204.29 eggs/g of body weight), the North Tyrrhenian (202.35 eggs/g of body weight) [11], and in the Eastern-Central Atlantic (228.33 eggs/g of body weight) [12], yet it was lower than the value estimated in the Central Tyrrhenian and Southern Adriatic Sea (281.6 eggs/g of body weight GSA10 and 262.2 eggs/g of body weight GSA18) [16]. The batch fecundity estimation ranged between 7771 and 137,256 eggs. It was lower compared to the Southern Adriatic Sea and Central Tyrrhenian, probably due to the major dimensions of reproductive animals in these areas compared to Northern and Central Adriatic Sea [16]. In fact, as reported by Recasens et al. [11], El Habouz et al. [12], Carbonara et al. [16] and Korta et al. [53], the relationship between fecundity estimations and hake size influences the total of eggs number per female that increases with size, weight, and gonad weight.

In conclusion, the scenario depicted in the present study evidenced that smaller mature females released fewer eggs per spawn than bigger mature females. As reported by Working Group on Stock Assessment of Demersal Species (WGSAD) [54], the European hake stock of Adriatic Sea results to be overfished and the overexploitation together with a low number of spawned eggs, and unfavorable environmental events could lead to a consequent decrease of the fish resource with the collapse of the stock. To avoid the collapse of European hake, especially in the Central-North Adriatic Sea, continuous monitoring programs, estimation of reproductive potential, and the stock assessment are essential. The easiest, cheapest, and most direct tool to investigate the reproductive stage of European hake, necessary for the determination of reproductive potential, could be represented by the maturity scale provided in this study. Furthermore, the reproductive cycle, the analysis of somatic indices, the L50 estimations, and the estimation of fecundity represent valuable information for a scientific decision-making process to establish suitable management measures aimed at tackling the continuing stock decline of European hake in an overfished area as Central-North Adriatic Sea. Finally, the result highlighted in the present study suggested increasing the minimum legal catch size as a tool to preserve this resource in this area.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/jmse9040389/s1.

**Author Contributions:** M.C. contributed to the acquisition of data, analysis of data, and drafting of the manuscript. S.C. was involved in contributions to conception and design of the manuscript, acquisition of data, and analysis of data and in revising it critically for important intellectual content. L.M. contributed to the analysis and interpretation of data. A.S. contributed to the conception and design of the project and to revising it critically for important intellectual content. D.B., G.R., and G.G. contributed to the analysis of data and to revising it critically for important intellectual content. O.C. conceived the experimental design, contributed to the analysis of data, and revised it critically for important intellectual content. All authors gave the final approval of the version to be published. Each author has participated in the work to take public responsibility for appropriate portions of the content and agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the Italian Ministry of Agricultural, Food, and Forestry Policies (MiPAAF) and European Union (Italian National Programs 2011–2013 and 2014–2016, in the ambit of Data Collection Framework).

**Institutional Review Board Statement:** The sampling procedures did not include any animal experimentation and animal ethics approval was therefore not necessary, under the Italian legislation (D.L. 4 March 2014, n. 26, art. 2).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request due to restrictions, e.g., privacy or ethical. The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors wish to thank the Filippo Domenichetti and Camilla Croci of Nation-al Research Council (CNR) Institute of Biological Resources and Marine Biotechnologies (IRBIM) Ancona, Italy, and Captain Giordano and crew of the "Orizzonte" vessel for their support in sampling.

**Conflicts of Interest:** Sabrina Colella, Alberto Santojanni, and Oliana Carnevali contributed equally to this work. The authors declare that there is no conflict of interest.

#### **References**


*Article*

## **Bioaccumulation of Metals**/**Metalloids and Histological and Immunohistochemical Changes in the Tissue of the European Hake,** *Merluccius merluccius* **(Linnaeus, 1758) (Pisces: Gadiformes: Merlucciidae), for Environmental Pollution Assessment**

**Antonio Salvaggio <sup>1</sup> , Roberta Pecoraro <sup>2</sup> , Chiara Copat <sup>3</sup> , Margherita Ferrante <sup>3</sup> , Alfina Grasso 3, Elena Maria Scalisi 2, Sara Ignoto 2, Vincenza Serena Bonaccorsi 2, Giuseppina Messina 2, Bianca Maria Lombardo 2, Francesco Tiralongo 2,**† **and Maria Violetta Brundo 2,\*,**†


Received: 5 August 2020; Accepted: 11 September 2020; Published: 15 September 2020

**Abstract:** Pollution and other types of environmental stress do not spare marine environments, especially those affected by high industrial pressure. Fish, especially coastal species, are used for monitoring the marine environment because they are particularly efficient as bioindicators thanks to their ability to bioaccumulate and biomagnify along the trophic chain. The aim of this research is to evaluate the bioaccumulation and the indirect bioindication ability of the European Hake, *Merluccius merluccius* (Linnaeus, 1758), one of the most important commercial fish species of the Mediterranean Sea. Morphological and histological alterations of the main target organs, such as liver and gills, have been investigated and the results showed a steatosis in the hepatic tissue. The accumulation of heavy metals has been analyzed by inductively coupled plasma mass spectrometry and for several metals it was showed a different concentration in the two sexes. Moreover, the expression of metallothioneins 1 and Heat Shock Protein 70 has been assessed by immunohistochemistry and did not show high level of expression. We underline the importance of contamination evaluation in commercial fish species and the utilization of the ichthyofauna as bioindicator of environmental quality.

**Keywords:** commercial fish species; Mediterranean Sea; environmental health; heavy metals; biomarkers

#### **1. Introduction**

The significant industrial and technological development recorded in recent years has brought undoubted benefits both in terms of quality increase and life expectancy. Despite such benefits,

this rapid progress has greatly and negatively increased the impact on the environment and on human health.

An economy based on encouraging the consumption of goods and services promotes a progressive increase in anthropogenic pressure on the environment, due to the increasing loads of substances of all kinds released into the air, water and soil. Exposure to these contaminants can cause harmful effects on the health of organisms and humans.

One of the most dramatic aspects that humans will have to face himself in the immediate future is represented by waters pollution, determined by the growing release of heavy metals which, being unable to be degraded or destroyed, tend to accumulate in the environment and in organisms with consequent and substantial risks both for the various living forms and for humans who consume foods with a content of significantly high metals [1].

The effects of heavy metal pollution on the environment, on the conservation of biodiversity and on human health are now well known [1]. Immediate actions are therefore necessary both aimed at preventing the continuation of uncontrolled discharges and at imposing compliance with current regulations in terms of protection of the environment and public health. These interventions should be preceded by an environmental biomonitoring program that allows rapid identification of the risk areas.

In biomonitoring programs, the choice of species to be used is particularly important. In recent decades among the new generation bioindicators, fish have been taking on a growing interest in assessing environmental quality in various aquatic ecosystems, for their ability of bioaccumulation and biomagnification along the trophic chain [2,3]. It is also necessary to remember that fish are a fundamental component in human nutrition, representing a significant source of protein, polyunsaturated fatty acids and micronutrients. However, through the consumption of fish products, humans are thus exposed to various contaminants.

Heavy metals are elements with a characteristic shiny appearance and a good electric conductivity. In chemical reactions, heavy metals frequently occur as cations [4]. However, non-heavy elements, such as aluminum (density 1.5), can be dangerous under particular conditions. Indeed, it becomes a powerful toxic if placed in acidic waters [4].

From an environmental point of view, the term heavy metals are commonly referred to all metals or non-metals which represent a danger to health and to the environment. They are substances present in the environment as constituent elements of the same, as well as introduced following industrial emissions. They can undergo biological and chemical transformations, which entail their accumulation in the environment and in organisms, both vegetable and animal, thus manifesting their deleterious action. Although currently exposure to anthropogenic sources is of prevailing toxicological importance, exposure to natural sources has proven to be fundamental for the development, in living organisms, of detoxification mechanisms, elimination and use aimed at reducing the danger of metals. These mechanisms allow some animal species to withstand high metal concentrations, which vice versa can be toxic to others, without suffering any damage [5]. Some metals are required by organisms in limited quantities; in particular, Zn, Cu, Fe, and Mg, even if present at low concentrations, perform a series of fundamental activities for the cell behaving as essential micronutrients and participating in numerous biochemical processes responsible for cell growth and life [6].

Zinc, for example, is involved in the replication, transcription and translation processes [7,8], acts as a cofactor for over 200 metalloenzymes [9] and performs regulatory functions, as in the case of modulation of synaptic transmission [10]. Copper, in low concentrations, is essential for breathing, for defense against free radicals and for the synthesis and release of neurotransmitters [11]. Other metals, such as Cd, Cr, Al, Hg and Pb are not normally present in the cells not even in traces and therefore, even at low concentrations, they can be very toxic. These metals can cause delays in embryonic development [12,13], in growth [14–17], as well as a long series of pathologies, including cancer [18,19]. Fish species, for example, hardly eliminate the absorbed mercury and the metal halving times vary from 6 months for mussels to 2 years for pike. The accumulation in fish is greater in muscle tissue than in adipose tissue and about 90–99% of the mercury present in fish is in the form of

methylmercury, an extremely toxic form of this metal [20]. Therefore, fish in risk assessment studies are useful bioindicators and can represent an early warning system of environmental damage, which can also be used for the assessment of potential risks to human health.

*Merluccius merluccius* (Linnaeus, 1758) is a demersal fish whose depth range normally extends from 70 to 400 m, although it can be found in shallower waters and up to about 1000 m depth. Its distribution range extends from the eastern Atlantic (from Norway and Iceland to Mauritania) to the Mediterranean Sea and southern part of the Black Sea. Unlike the small specimens (<14 cm TL), which feed mainly on euphausiids and mysids, the larger specimens (>32 cm TL) of the European Hake, *M. merluccius*, are ichthyophagous [21]. Furthermore, this species has been successfully used as bioindicator [22]. Hence, considering the commercial importance of this predator, its position in the trophic web and its ability to be used as a good bioindicator of marine water pollution, our choice to use *M. merluccius* in the current study.

In this perspective, our work aims to evaluate the impact of metals/metalloids contamination in the commercial fish species *M. merluccius*, commonly known as "European Hake", sampled in the Ionian Sea. Previous research in this area of study report a total metal load higher in pelagic fish than demersal and benthic ones [23,24]. Although fish species resulted stressed by environmental conditions with a certain degree of oxidative stress in liver tissue [23], from a chemical point of view, the analyzed species were healthy for human consumption and the human risk to develop chronic systemic and carcinogenic effects due to their consumption was low [24]. Nevertheless, studies focused on the biomonitoring of the Sicilian ionic coast highlighted a significant metal load in areas subject to heavy anthropogenic pressure, particularly the Augusta coastal water, where it is hosted the largest industrial complex of Sicily [25,26].

The present study has been based on analysis of the gastrointestinal content to evaluate the presence of metals/metalloids; quantitative and qualitative analysis of metals/metalloids in the tissues; analysis of the toxicological effects of metals/metalloids through biomarkers detection (Heat Shock Proteins 70 and Metallothioneis 1) and histological analysis.

#### **2. Materials and Methods**

For this study, 20 specimens (10 males and 10 females) of *M. merluccius* fished in Food and Agriculture Organization of the United Nations (FAO) area 37 (Marzamemi, southeastern Sicily, Ionian Sea) were analyzed. Fish size ranged from 35 to 45 cm total length (TL), and were captured with longlines by local fishermen and transported fresh to the laboratory for analyses.

#### *2.1. Metals and Metalloids Analysis*

The method for extraction and quantification of metals and metalloids is described by Copat et al. [27]. Briefly, aliquots of 0.5 g of muscle, liver, gills and gastrointestinal content of fish were acid digested with 6 mL of 65% nitric acid (HNO3) (Carlo Erba) and 2 mL of 30% peroxide hydrogen (H2O2-Carlo Erba) in a microwave system (Ethos Touch Control, Milestone S.r.l., Italy). Analytical determination of arsenic (As), cadmium (Cd), cobalt (Co), chromium (Cr), copper (Cu), lead (Pb), mercury (Hg), manganese (Mn), nickel (Ni), vanadium (V), selenium (Se), antimony (Sb) and zinc (Zn) was performed with an ICP-MS Elan-DRC-e (Perkin–Elmer, United States). Blanks, standard and samples were prepared with the same reagents. A multi-elements certified reference solution ICP Standard (Merck) was used for the instrument calibration. Processed blanks were used to calculate the method detection limits (MDL) based on the following equation:

$$\text{MDL} = \text{One-tailed student's } t \text{-test (} p = 0.99\% \text{); } \text{df} = \text{n } -1 ) \times \text{Sr}$$

MDL (mg/kg ww) estimated for each trace elements are the following: As 0.013, Cd 0.002, Co 0.008, Cr 0.003, Cu 0.005, Pb 0.001, Hg 0.0025, Mn 0.005, Ni 0.007, V 0.025, Se 0.03, Sb 0.020 and Zn 0.109.

The quality control was performed with laboratory-fortified matrix (LFM) processed at each batch of digestion, obtaining a recovery ranges from 91.5 to 110% of the nominal concentration.

Statistical analysis was performed with the software SPSS (version 20.0, Inc., IBM, Armonk, NY, USA: IBM Corp.). Results below the MDL were elaborated as MDL/2. The normal distribution was verified using the Kolmogorov–Smirnov test. Since the low number of samples, the Mann-Whitney non-parametric test was used to compare median concentrations between tissues.

#### *2.2. Histological Analysis*

Histological analysis was performed according to our standard laboratory procedures [28]. Liver and gills were fixed in 4% formaldehyde (Bio-Optica, Milano, Italy) in PBS buffered to 0.1 M, pH 7.4 (Sigma Life Science) at room temperature for 48 h and processed with Tissue Processing Center TPC 15 Duo (MEDITE®, Burgdorf. Germany). The sections were stained with Haematoxylin-Eosin (HE) (Bio-Optica) and observed under optical microscope (Leica DM750, Monument, CO, USA) equipped with a digital camera (Leica DFC500, Monument, CO, USA).

#### *2.3. Immunohistochemical Analysis*

The immunohistochemical protocol was performed on sections to detect mouse monoclonal anti-HSP70 (Gene Tex, 1:1000) and to detect mouse polyclonal anti-MT1 (Abcam, 1:1000); secondary antibody used is FIT-conjugated goat anti-mouse IgG (Sigma-Aldrich, 1:1000). Analysis were performed according to our standard laboratory procedures [29–31]. Slides after mounted with mounting medium containing DAPI (Vectashield, Vector Laboratories, Burlingame, United States), were observed with NIKON ECLIPSE Ci fluorescence microscope and the images taken with the NIKON DS-Qi2 camera.

#### **3. Results and Discussion**

The chemical analysis showed that the bioaccumulation of metals in gills was higher in males for Cu (*p* < 0.001), Hg (*p* < 0.01), Se (*p* < 0.001) and Zn (*p* < 0.001), and in females for Cd (*p* < 0.001), Co (*p* < 0.001), Cr (*p* < 0.01), Pb (*p* < 0.001) and V (*p* < 0.05) (Table 1, Figure 1). Nevertheless, the concentrations of all the metals examined were low in both sexes, with higher concentrations of the essential metals Zn and Mn. The analysis of the gastrointestinal content showed a predominance content of Co (*p* < 0.05) and Mn (*p* < 0.01) in males versus an higher concentration of As (*p* < 0.01), Pb (*p* < 0.001), Ni (*p* < 0.01), V (*p* < 0.01) and Se (*p* < 0.001) in females (Table 1, Figure 2). In muscle tissue of both sex, it was observed a lower content of all metals versus the concentrations find in gill and gastrointestinal content. In females, it was found a higher concentration of Cr (*p* < 0.01), Cu (*p* < 0.05), Mn (*p* < 0.001), Ni (*p* < 0.001) and Zn (*p* < 0.001) than males. In males, only the bioaccumulation of Pb (*p* < 0.001) was found higher than females (Table 1, Figure 3).

The results obtained agreed with Bosch et al. [32] which describes a differential bioaccumulation. Particular attention is given to some heavy metals: cadmium, mercury and lead.

The authors hypothesize that the concentrations of heavy metals depend on factors influencing the absorption of metals, such as the simultaneous presence of other xenobiotics, the geographical distribution and the specific biological factors of the species. Overall, the metal concentrations found by us are comparable to those of other studies present in literature [33–35]. These concentrations of metals present in the edible portion of the fish (i.e., muscle) is acceptable as it falls within the limits of international legislation and the specimens analyzed seem to be safe for human consumption. Lead, cadmium and mercury were recorded with values well below the limits imposed by law (0.30 mg/kg for lead, 0.05 mg/kg for cadmium, 0.50 mg/kg for mercury).

*J. Mar. Sci. Eng.* **2020**, *8*, 0712



**Figure 1.** Mean values of metals and metalloids concentrations in gills (mg/kg ww). Significant differences between sex: a, *p* < 0.001; b, *p* < 0.01, c: *p* < 0.05.

**Figure 2.** Mean values of metals and metalloids concentrations in gastrointestinal content (mg/kg ww). Significant differences between sex: b, *p* < 0.01, c: *p* < 0.05.

As far as histological investigations are concerned, morphological anomalies of the gill lamellae were not found in the analyzed specimens (Figure 4A,B), but marked liver steatosis was evident in the liver sections (Figure 4C,D). The immunohistochemical investigation revealed a weak expression of MT1 (Figure 5A,C) and HSP70 (Figure 5B,D) in both the target organs (gills and liver) of both sexes.

**Figure 3.** Mean values of metals and metalloids concentrations in muscle (mg/kg ww). Significant differences between sex: a, *p* < 0.001; b, *p* < 0.01, c: *p* < 0.05.

**Figure 4.** Sections stained with Hematoxylin-Eosin. (**A**,**B**); gills sections (**C**,**D**); male liver sections, is evident a diffuse steatosis. Scale bar A,C and D: 200 μm; B 100 μm.

**Figure 5.** MT 1 expression in male gills (**A**) and liver (**C**). HSP70 expression in male gills (**B**) and in liver (**D**). Using specific antibodies anti-MT1 and anti-HSP70 it has been showed a weak expression of these proteins in cells (green). Scale bar A, C and D: 200 μm; B: 50 μm.

The histological analysis, despite being non-specific about pathogenic noxae, offers the possibility of highlighting any alterations in the gill and liver tissues. The use of these tissue is an efficient biomarker in a preliminary screening for fish and water quality. *Merluccius merluccius*, being a voracious predator, can be considered a useful species for monitoring the environmental situation over time through bioaccumulation and biomagnification studies. Heavy metals contamination in fishery products remains an almost inevitable condition due to the environmental situation facing most of the national fishing areas. Therefore, as already pointed out by European Food Safety Authority (EFSA) [36] in order to reduce this risk to an acceptable level, it is necessary to bridge the lack of information on the conscious consumption of the various fish species and, more generally, it is necessary that consumers can receive a concrete information on the possible risks associated with the consumption of certain species.

In recent decades, with the rapid evolution of molecular, biochemical and pathological technologies, biomarkers have found a very wide application. The toxicity of a contaminant on an organism is usually expressed at a biochemical and molecular level and, because of this, at cellular and tissue levels. As a response to the action of a toxic agent, the body develops adaptive responses that tend to bring the system back into a balanced state. However, when the homeostatic mechanism is not sufficient to balance the action of the toxic agent, the negative effect is manifested at the cellular, tissue and eventually organ level. It can therefore be said that the different responses, homeostatic and otherwise, generated by an organism towards a toxic agent, represent potential markers that can be used in ecotoxicological investigations [37,38].

It can be said that the gills, because of the role of primary importance in interfacing the fish with the surrounding aquatic environment and for their demonstrated reactivity and sensitivity, are optimal candidates for the role of versatile biomarkers in the biomonitoring of aquatic ecosystems. They lend themselves to a study in the field of environmental protection, also aimed at protecting the human being as a user, direct and indirect, of the resources that natural ecosystems can offer. However, also the investigation of ecological indices represents a good method for the monitoring of the environmental quality [39] and these types of studies should be associated to the bioaccumulation studies for a more complete and wide point of view.

#### **4. Conclusions**

The diversity of cultures bordering the Mediterranean and the level of anthropization present make interventions aimed at environmental protection difficult. The same international fishing legislation of the Mediterranean Sea is applied with different weights, according to the different social and economic realities on which it acts. To protect the Mediterranean environment, which currently has different pressures at the expense of coastal and marine habitats, it is increasingly necessary, in addition to implementing and applying existing environmental laws, to also adopt integrated approaches based on knowledge of ecosystems in order to better understand their organization and functioning. The result of these negative impacts is an ecosystem in which the survival of fish species is dependent on human intervention. This study analyzes some of the complex variables that currently affect the fish populations of the Mediterranean Sea, intended both as a fish resource of great commercial value, and as a symbol of an ecosystem that now lives below the limits of sustainability. In particular, the contamination level recorded by us in *M. merluccius* in the current study suggested that this important commercial species can be considered sure for human consumption, although we recorded a marked liver steatosis that can be the result of stressing environmental conditions. Further investigation in this direction may help to better understand contamination dynamics and effects on fish population and their potential risks to human health.\

**Author Contributions:** All authors have made substantial contributions to the conception and the design of the work; have approved the submitted version have agree to be personally accountable for the author's own contributions and for ensuring that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and documented in the literature. Methodology, A.S., V.S.B., C.C., A.G.; visualization, S.I., G.M., M.F., B.M.L.; writing, review and editing, M.V.B., F.T., R.P., E.M.S., A.S.; conceptualization M.V.B. and F.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the Department of Biological, Geological and Environmental Science and by the Department of Medical Sciences, Surgical and Advanced Technologies "G.F. Ingrassia"—Hygiene and Public Health, University of Catania.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18